數(shù)值分析實驗報告

1、天津工業(yè)大學數(shù)值分析實驗報告電子與信息工程學院例題一用Gauss消去法來解方程組x1+x2+x3=612x1-3x2+3x3=15-18x1+3x2-x3=-15模型：高斯（Gauss）消去法高斯消去法的基本思想是：先逐次消去變量，將方程組化成同解的上三角形方程組，此過程成為消元過程。然后按方程相反順序求解上三角方程組，得到原方程組的解，此過程稱為回代過程。這種方法稱為Gauss消去法。將方程組改成如下形式a111x1+a121x2+······+a1n1xn=b11a211x1+a221x2+···&#

2、183;··+a2n1xn=b21: : : : : : : : an11x1+an21x2+······+ann1xn=bn1（2-1）消元過程：第一步：設a1110，記li1=ai1/a111(i=2,3,······，n)，方程中第i個方程減去第一個方程乘以li1，完成第一次消元。將上述方程式可以改寫為a111x1+a121x2+······+a1n1xn=b11a222x2+·&

3、#183;····+a2n2xn=b22 : : : : : : an22x2+······+ann2xn=bn2（2-2）其中aij(2)=aij(1)-li1a1j(1)，bi(2)=bi(1)-li1b11(i,j=2,3·········，n)。方程組（2-2）簡記為A2x=b2.第二步：設a2220，記li2=ai22/a222(i=2,3,····

4、3;·，n)，方程中第i個方程減去第一個方程乘以li1，完成第二次消元。第n-1步：同上，得如下方程組a111x1+a121x2+······+a1n1xn=b11a222x2+······+a2n2xn=b22 : : : : : ann2xn=bn2簡記A(n)x=b(n)回代過程：按變量的逆次序逐步回代得到方程組的解xn=bn/annnxk=(bkk-l=k+1naklkx1)/akkk (k = n-1,n-2···，1)算法：

5、1. 輸入A=(aij)，b=(b1·······bn)，n，置k=1。2. 若akk0,轉3；否則輸出失敗信息，停機。3. 對,置i=k+1···，n置 aik/akkb1對bi-aikbkaij4. 若k=n -1,轉5；否則k+1k，轉2,。5. 若ann=0，輸出失敗信息，停機；否則，置 bn/annbn。6. 對k = n-1,n-2···，1置bk-l=k+1naklbl)/akkkbk7. 輸出x=b,停機代碼：步驟一：編寫高斯消去法函數(shù)functio

6、n X =Fgauss(A ,b)zengguang =A b;n=length(b);ra=rank(A);rz=rank(zengguang);temp1=rz - ra;if temp1>0, disp('錯誤')returnendif ra=rzif ra=n; X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:n m=zengguang(k,p)/zengguang(p,p); zengguang(k,p:n+1)=zengguang(k,p:n+1)- m*zengguang(p,p:n+1);endend

7、b=zengguang(1:n,n+1); A=zengguang(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q);endelse disp('錯誤')endend步驟二：編寫方程組信息a=1 1 112 -3 3-18 3 -1;b=6 15 -15'x=Fgauss(a,b)步驟三：運行主函數(shù)得到結果：X1=1.000X2=2.000 X3=3.000結果分析：高斯消去法簡單易行，但其計算過程中要求akk0，因而適用范圍小，只能適用于從1到n-1階順

8、序主子式不為零的矩陣A，計算實踐表明高斯消去法數(shù)值穩(wěn)定性差，當出現(xiàn)小主元素時，會嚴重影響計算精度，甚至導致錯誤的結果。例題二用雅克比（Jacobi）迭代法求解下列方程組10x1-x2-2x3=72-x1+10x2-2x3=83-x1-x2+5x3=42模型：雅克比（Jacobi）迭代法因為方程組a11x1+a12x2+···+a1nxn=b1a21x1+a22x2+···+a2nxn=b2 : : : : : : : :an1x1+an2x2+···+annxn=bn（3-1）的系數(shù)矩陣A非奇異，不妨設a

9、ii0(i=1,2···,n)。將方程組（3-1）變形為x1= b12x2+···+b1nxn+g1x2=b21x1 +···+b2nxn+g2 : : : : : : : :xn=bn1x1+bn2x2+···+bnn-1xn-1+gn（3-2）其中bij=-aijaii，ij，i，j=1，2，···，n，gi=biaii(i=1，2，···，n)。若記B = 0 b12b13 ···

10、 b1n-1b1nb22 0 b23 ··· b2n-1b2n：bn1bn2bn3 ··· bnn-1 0g=g1g2 ··· gnT則方程組（3-2）可以簡記為x=Bx+g（3-3）選取初始向量x(0)代入迭代公式x(k+1)=Bx(k)+g(k=0,1,2，···)（3-4）產生向量序列x(k)，由上述計算過程所給出的迭代法稱為Jacobi迭代法，又叫簡單迭代法。式（3-4）為它的迭代公式，迭代矩陣為B。這是Jacobi迭代法的計算公式。如果用系數(shù)矩陣A來表示，則有B=I

11、-D-1A，g=D-1b（3-5）其中D=diag（a11，a22，···，ann），D-1=diag(1a11，1a22，···，1ann)。算法：1. 輸入A=（aij），b=（b1，b2，···，bn），維數(shù)n，x0=(x10，x20，···，xn(0)，最大容許迭代次數(shù)N。2. 置k=13. 對 i=1，2，···，n xi=(bi-j=1jinaijxj0)/aij4. 若x-x(0)<，輸出x，停機；否則轉5。5. 若k<

12、N置k+1k，xixi0(i=1,2，···，n)，轉3；否則，輸出失敗信息，停機。程序：步驟一：編寫雅克比迭代法函數(shù)function x,k = Fjacobi( A,b,x0,tol )max1=300; D=diag(diag(A); L=-tril(A,-1); U=-triu(A,1); B=D(L+U); f=Db;x=B*x0+f;k=1; while norm(x-x0)>=tolx0=x;x=B*x0+f;k=k+1;if(k>=max1) disp('錯誤');returnendk ,x'end步驟二：編寫方

13、程組信息a=10 -1 -2;-1 10 -2;-1 -1 5;b=72 83 42'x0=0 0 0'x,k=Fjacobi(a,b,x0,1e-7)步驟三：運行主函數(shù)得到結果結果分析：由表所示方程組的精確解為x=(11 12 13 )T當?shù)螖?shù)增加時，迭代結果越來越接近精確解。在迭代到11次后數(shù)值不變k x1(k) x1(k) x1(k)2.0000 9.7100 10.7000 11.50003.0000 10.5700 11.5710 12.48204.0000 10.8535 11.8534 12.82825.0000 10.9510 11.9510 12.9414

14、6.0000 10.9834 11.9834 12.98047.0000 10.9944 11.9944 12.99338.0000 10.9981 11.9981 12.99789.0000 10.9994 11.9994 12.999210.0000 10.9998 11.9998 12.999711.0000 10.9999 11.9999 12.999912.0000 11.0000 12.0000 13.000013.0000 11.0000 12.0000 13.000014.0000 11.0000 12.0000 13.000015.0000 11.0000 12.0000 1

15、3.000016.0000 11.0000 12.0000 13.000017.0000 11.0000 12.0000 13.000018.0000 11.0000 12.0000 13.000019.0 11.0000 12.0000 13.0000例題三：已知函數(shù)表如下：-1 0 10.5 1 2利用線性插值和二次插值計算的近似值。模型：Lagrange插值函數(shù)在區(qū)間上個互異節(jié)點的函數(shù)值。求一個至多次多項式使其在給定點處與同值，即滿足插值條件。線性插值，即只有兩個節(jié)點的插值多項式。插值多項式設為且滿足插值條件 ，解此方程組得，所以兩點的一次插值多項式為（3-1）記，則（3-1）可寫成這是

16、用兩點的直線近似曲線故這種插值又稱線性插值。由此得到啟發(fā)，當節(jié)點增加到時，可以先構造次多項式，它們滿足，插值基函數(shù)為插值多項式為（3-2）式（3-2）稱為次Lagrange插值多項式，為了以后便于區(qū)別，常用代替以突出表示這是由Lagrange插值所得到的插值多項式，即 （3-3）特別地，=1時，兩點一次Lagrange插值（也叫線性插值）多項式為；=2時稱為二次插值，也叫拋物線插值，插值多項式為算法：（1） 輸入。（2） 對（3）（4） 輸出停機。程序：function yp=mlagr(x,y,xp)n=length(x);m=length(xp);yp=zeros(1,m);c1=ones

17、(n-1,1);c2=ones(1,m);for i=1:n xb=x(1:i-1,i+1:n); yp=yp+y(i)*prod(c1*xp-xb'*c2)./(x(i)-xb'*c2);end結果分析：在命令窗口輸入x=0 1;y=1 2; xx=0.3;yy1=mlagr(x,y,xx)輸出結果為yy1 = 1.300000000000000這是線性插值的求得的的近似值。在命令窗口輸入 x=-1 0 1; y=0.5 1 2; xx=0.3; yy2=mlagr(x,y,xx)輸出結果為yy2 = 1.247500000000000這是拋物線插值求得的的近似值。的準確解為

18、1.2311，可見拋物線插值比線性插值的誤差要小。這是因為拋物線插值選取了比線性插值更多的節(jié)點求解近似值。由插值余項可知，Lagrange插值法要提高插值多項式的函數(shù)逼近程度，要增加節(jié)點數(shù)，提高多項式的次數(shù)，但這樣往往不能達到預想的效果。為了克服Runge現(xiàn)象，我們常用分段線性插值來改進。例題四：分別用復化梯形公式和復化Simpson公式計算下列積分的近似值。模型：復化梯形公式和復化Simpson公式（1）復化梯形公式用分段線性插值函數(shù)近似被積函數(shù)，等于把積分區(qū)間分成若干小區(qū)間，在每個小區(qū)間上以梯形的面積近似曲邊梯形的面積。即用梯形公式求小區(qū)間上積分的近似值。如圖1圖1將積分區(qū)間等分，記，。在

19、每個小區(qū)間上用梯形公式求和，得整理得 （1）式（1）稱為復化梯形公式。如果，在小區(qū)間上，梯形公式的截斷誤差為因此因為在連續(xù)，有介值定理，存在使得從而有 （2）這就是復化梯形公式的截斷誤差。（2）復化Simpson公式如果用分段而插值函數(shù)近似被積函數(shù)，即在小區(qū)間上用Simpson公式計算積分近似值，就導出復化Simpson公式。將區(qū)間分成等分，分點為，在每個小區(qū)間上，用Simpson公式求積分，則有求和得整理后得到 （3）式（3）稱為復化Simpson公式。如果復化Simpson公式的截斷誤差為因為連續(xù)，故存在使得代入上式得 （4）式（4）表明，步長越小，截斷誤差越小。當時，用復化Simpson

20、公式所求的近似值收斂域積分值，而且算法具有數(shù)值穩(wěn)定性。算法：a.復化梯形公式（1） 輸入f(x)，a，b，區(qū)間等分數(shù)n（2） .（3）（4）（5） 停機。b.復化Simpson公式（1） 輸入f(x)，a，b，區(qū)間等分數(shù)n，n=2m（2） ,（3）（4） 停機。程序：a.復化梯形公式function s=mtrap(f,a,b,n)h=(b-a)/n;x=linspace(a,b,n+1);y=feval(f,x);s=0.5*h*(y(1)+2*sum(y(2:n)+y(n+1);endb.復化Simpson公式function s=msimp(f,a,b,n)h=(b-a)/n;x=lin

21、space(a,b,2*n+1);y=feval(f,x);s=(h/6)*(y(1)+2*sum(y(3:2:2*n-1)+4*sum(y(2:2:2*n)+y(n+1);end結果分析： 在MATLAB工作窗口輸入format longf=(x)x./(1+x.2);s=mtrap(f,0,1,8)p=msimp(f,0,1,8)輸出結果為s =0.345268948437223p =0.344490897538895s和p分別是復化梯形公式、復化Simpson公式對在的近似值。例題五確定方程x3-x+4=0的實根的分布情況，并用二分法求在開區(qū)間（-2，-1）內的實根的近似值，要求精度為0

22、.001.模型：二分法二分法也稱對分區(qū)間法。它的基本思想是：先確定方程f(x)=0含根的區(qū)間(a,b)，再把區(qū)間逐次二等分。設f(x)在區(qū)間a,b上連續(xù)，且有fafb<0則由連續(xù)函數(shù)介值定理，f(x)在(a,b)內必有零點，稱(a,b)為方程f(x)=0的有根區(qū)間。不妨設fa<0，fb>0，取x0=a+b2，若fx0=0,則x=x0就是方程f(x)=0的解。否則，若fx0<0，取a1=x0，b1=b；若fx0>0，取a1=a，b1=x0，則有 a1,b1a,b，b1-a1=b-a2且f(x)在a1,b1上連續(xù)，滿足fa1fb1<0。重復上述過程又得到區(qū)間a2

23、,b2，滿足a2,b2a1,b1，b2-a2=b1-a12，且fa2fb2<0。如此繼續(xù)下去，得到一個區(qū)間序列a,ba1,b1a2,b2an,bn滿足fanfbn<0，bn-an=b-a2n因而滿足an,bn(n=1,2,)均為方程f(x)=0的有根區(qū)間。由區(qū)間套定理，存在x*a,b，使得limnan=limnbn=x*且x=x*是方程f(x)=0的根。若取xn=an+bn2作為根x*的近似值，則誤差為 x*-xnbn-an2=b-a2n+1 （1-1）按上述方法求非線性方程f(x)=0的近似解稱為對分區(qū)間法，式（1-1）表明，只要f(x)連續(xù)，對分區(qū)間法總是收斂的。算法：（1）輸

24、入f(x)，a，b，誤差限，最大容許迭代次數(shù)N。（2）置n=0。（3）xn=a+b2。（4）若fxn=0或b-a2<，則輸出x*=xn，停機；否則轉5。（5）若n<N，轉6；否則輸出失敗信息，停機。（6）若fafxn<0，則xnb；否則，xna。（7）置n+1n，轉3。（c）程序：function y1=fun(x) y1=x3-x+4;function k,x,wuca,yx=erfen(a,b,abtol)a(1)=a;b(1)=b;ya=fun(a(1);yb=fun(b(1);if ya*yb>0,disp('注意：ya*yb>0，請重新調整區(qū)間端

25、點a和b'),returnendmax1=-1+ceil(log(b-a)-log(abtol)/log(2);for k=1:max1+1 ya=fun(a); yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2;k=k-1; k,a,b,x,wuca,ya,yb,yxif yx=0 a=x; b=x;elseif yb*yx>0 b=x; yb=yx；else a=x; ya=yx;endif b-a<abtolreturnendendk=max1;x;wuca;yx=fun(x);endk,x,wuca,yx=erfen(

26、-2,-1,0.001)結果分析：首先先確定方程x3-x+4=0的根的分布情況，在MATLAB工作窗口輸入程序x=-4:0.1:4;y=x.3-x+4;plot(x,y)grid,gtext('y=x.3-x+4')畫出函數(shù)f(x)=x3-x+4的圖像，如圖1-1所示，從圖像中可以看出，此曲線有兩個駐點±33都在x軸的上方，在（-2，-1）內曲線與x軸只有一個交點，則該方程在（-2，-1）內唯一的一個實根。圖1-1其次，再在區(qū)間（-2，-1）內求出實根的近似值，在MATLAB的工作窗口輸入程序：k,x,wuca,yx=erfen-2,-1,0.001運行后屏幕上顯示二

27、分法的計算結果，如表1-2所示：表1-2次數(shù)k左端點ak右端點bk中點xkbn-an2函數(shù)值fak函數(shù)值fbk函數(shù)值fxk0-2.0000-1.0000-1.50000.5000-2.00004.00002.12501-2.0000-1.5000-1.75000.2500-2.00002.12500.39062-2.0000-1.7500-1.87500.1250-2.0000.3906-0.71683-1.8750-1.7500-1.81250.0625-0.71680.3906-0.14184-1.8125-1.7500-1.78130.0313-0.14180.39060.12965-1

28、.8125-1.7813-1.79690.0156-0.14180.1296-0.00486-1.7969-1.7813-1.78910.0078-0.00480.12960.06277-1.7969-1.7891-1.79300.0039-0.00480.06270.02908-1.7969-1.7930-1.79490.0020-0.00480.02900.01219-1.7969-1.7949-1.79590.0010-0.00480.01210.0037經過計算，我們可以得出方程根的近似值為x=-1.7959，二分區(qū)間的次數(shù)為k=9，誤差為wuca=9.7656e-004?？梢姡址?/p>

29、的優(yōu)點是對函數(shù)的要求低，方法簡單、可靠、程序設計容易，可以容易估計計算次數(shù)，收斂速度恒定；缺點是不能求出方程的偶重根，收斂速度慢。例題六用迭代法求方程fx=x2+2x-10的一個正根。模型：迭代法迭代法的基本思想是利用某種迭代公式，使某個近似根逐步精確化，直到滿足精度要求的近似根為止。把給定的方程f(x)=0改寫成x=x（其中x稱為迭代函數(shù)），選擇合適的初始值x0，代入x算得的結果記為x1=x0，一般x1x0；把x1代入x算得的結果記為x2=x1,.若把xk代入x，一般有迭代公式 xk+1=xk(k=0,1,2,) (2-1)由此得到序列x0,x1,xk,.若迭代序列xk收斂到x*（及l(fā)imk

30、xk=x*），則當函數(shù)x連續(xù)時，由(2-1)得limkxk=limkxk=limkxkx*=x*，稱為的不動點，即x*為原方程f(x)=0的根。一般計算有限步，即可得到某種精度的近似根xk。這種求方程根的方法稱為簡單迭代法。算法：（1） 輸入x，初始值x0，誤差限，最大容許迭代次數(shù)N。（2） 置n=1。（3） 計算x=x0。（4） 若x-x0<或x-x0x<，則輸出x*=x，停機；否則轉5。（5） 若n<N，則置n+1n，xx0，轉3；否則輸出失敗信息，停機。程序：function y1=fun1(x)y1=g(x);functionk,piancha,xdpiancha,xk=diedai(x0,k)x(1)=x0;for i=1:kx(i+1)=fun1(x(i);piancha=abs(x(i+1)-x(i);xdpiancha=piancha/(abs(x(i+1)+eps);i=i+1;xk=x(i);(i-1)piancha xdpiancha xkendif(piancha>1)&&(xdpiancha>