1、數值分析實 驗 指 導 書濰坊學院數學與信息科學學院2012年04月文檔可自由編輯打印目 錄目錄I實驗一 插值與曲線擬合的最小二乘法1實驗二 數值積分4實驗三 解線性方程組的直接法9實驗四 解線性方程組的迭代法11實驗五 非線性方程的數值解法13實驗六 常微分方程數值解法17實驗一 插值與曲線擬合的最小二乘法一、實驗目的：1了解拉格朗日插值法、牛頓插值法、曲線擬合最小二乘法的基本原理和方法；2掌握拉格朗日插值多項式牛頓插值多項式的用法；3掌握最小二乘原理，會求擬合函數及超定方程組的最小二乘解。二、實驗內容：1用拉格朗日插值公式和牛頓插值公式確定函數值；2對函數f (x)進行拉格朗日插值和牛頓插

2、值；3利用Polyfit擬合冪函數，利用Polyfit擬合多項式。三、實驗過程：1給定函數四個點的數據如下:，試用插值公式確定函數在處的函數值。MATLAB程序如下：X=1.1,2.3,3.9,5.1; Y =3.877,4.726,4.651 ,2.117;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv (p1, p3), p4)/( X(2)- X(

3、1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),

4、運行后輸出的基函數l0，l1，l2和l3為l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*x輸入程序 L=poly2sym (P),x=2.101; Y = polyval(P,x)運行后輸出插值多項式和插值為L=-629/5376*x3+31433/53760*x2-631/2856*x+21353/562949953421312Y =4.5969輸入程序 L=poly2sym (P),x=4.234; Y = polyval(P,x)運行后輸出插值多項式和插值為L

5、=-629/5376*x3+31433/53760*x2-631/2856*x+21353/562949953421312Y =4.2244L=1645/90992*x3-2517512191700115/45496*x2+14477/6000*x+2007/1000Y = 3.4290輸入程序 syms M; x=2.101; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24運行后輸出誤差限為R3 =4351/363968*M2在區間上取結點數，等距間隔的節點為插值點，對于函數進行拉格朗日插值。MATLAB程序如下t=-5: 1:5;ft=(1+t.*

6、t).5;t1=-5:1:5;ft1=(1+t1.*t1).5;y1=Lagran(t1,ft1,t);plot(t,ft,b:,t,y1,g+);xlabel(x);ylabel(y);對一組數據做拉格朗日的M文件如下Lagranmfunction fi =Lagran(x,f,xi)fi=zeros(size(xi);npl=length(f);for i=1:nplz=ones(size(xi);for j=1: npl if i=j z=z.*(xi-x(j)/(x(i)-x(j); endendfi=fi+z*f(i);end3給出節點數據，作三階牛頓插值多項式，計算，并估計其誤差。

7、MATLAB程序如下syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,A,C,P=newdscg(X,Y,x,M)運行后輸出y = 22.3211R =65133/562949953421312*M（即R =2.3503*M）A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167C =0.9167 4.2500 -4.1667 1.0000P =11/12*x3+17/4*x2-25/6*x+1其中求牛頓插值多項式、差商、插值及其

8、誤差估計的MATLAB主程序如下：function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y;s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C

9、,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);4給定數據如下 試用冪函數擬合以上數據。Matlab程序如下x=0.15,0.4,0.6,1.01,1.5,2.2,2.4,2.7,2.9,3.5,3.8,4.4,4.6,5.1,6.6,7.6y=4.4964,5.1284,5.6931,6.2884,7.0989,7.5507,7.5106,8.0756,7.8708,8.2403,8.5303,8.7394,8.9981,9.1450,9.5070,

10、9.9115c=polyfit(x, y,1)5用以下數據求二次多項式的系數，并做出擬合曲線： Matlab程序如下x=0.1,0.4,0.5,0.7,0.7,0.9;y=0.61,0.92,0.99,1.52,1.47,2.03;cc=polyfit(x,y,2)xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy);hold onplot(x,y,x)axis(0,1,0,3)xlabel(x)ylabel(y)四、實驗要求：1學會Lagrange插值方法和牛頓插值方法并應用上述方法于實際問題；2將擬合的結果與拉格朗日插值的結果比較。3歸

11、納總結數值實驗結果，試定性地說明函數逼近各種方法的適用范圍，及實際應用中選擇方法應注意的問題。實驗二 數值積分一、實驗目的：1理解數值積分的概念，掌握各種數值積分方法，包括梯形公式、拋物線公式等，通過實際計算體會各種數值積分方法的精確度；2掌握復化梯形公式、復化Simpson公式及其截斷誤差的分析；3了解Romberg公式及Romberg積分法。二、實驗內容：1利用低階求積公式求積分的近似值；利用復化梯形公式和復化Simpson公式計算定積分；2利用復化梯形公式計算無窮限廣義積分；3利用Romberg積分法計算定積分。三、實驗過程：1用求截斷誤差公式的MATLAB主程序，求計算定積分d的近似值

12、的階牛頓科茨公式的截斷誤差公式。Matlab程序如下n=1, RNC1=ncE(n), n=2, RNC2=ncE(n), n=3, RNC3=ncE(n)n=4, RNC4=ncE(n), n=8, RNC8=ncE(n)運行后屏幕顯示結果如下n= 1 RNC1 =1/12*(b-a)3*fxn1n = 2 RNC2 =1/90*(1/2*b-1/2*a)5*fxn2n = 3 RNC3=3/80*(1/3*b-1/3*a)5*fxn1n = 4 RNC4=8/945*(1/4*b-1/4*a)7*fxn2n = 8 RNC8=353/6926923254988800*(1/8*b-1/8*

13、a)11*fxn2其中計算n階牛頓-科茨的公式的截斷誤差公式的MATLAB主程序function RNC=ncE(n)suk=1; p=n/2-fix(n/2);if p=0for k=1:n+2suk=suk*k;endsuk; syms t a b fxn2,su=t2; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC= (b-a)/n)(n+3)*fxn2*abs(y)/ suk;elsefor k=1:n+1suk=suk*k;endsuk; syms t a b fxn1,su=t; for u=1

14、:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC= (b-a)/n)(n+2)*fxn1*abs(y)/ suk;end2分別利用復化梯形公式和復化Simpson公式計算定積分 取，精確值為。復化梯形公式MATLAB程序如下clear;Iexact=4.006994;a=0;b=2;fprintf(n n I Errorn);n=1; for k=1:4n=2*n;h=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=sqrt(1+exp(x);I=travez_v(f,h);fprintf(%3.0f %10.5

15、f %10.5fn,n,I,Iexact-I);end其中復化梯形求積法的M文件如下trapez_v.mfunction I=trapez_v(f,h)I=h*sum(f)-(f(1)+f(length(f)/2;復化Simpson求積公式MATLAB程序如下：clear;Iexact=4.006994;a=0;b=2;fprintf(n n I Errorn);n=1;for k=1:4,n=2*nh=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=sqrt(1+exp(x); I=(h/3)*(f(1)+4*sum(f(2:2:n)+f(n+1); if n2 I=I+(h/3

16、)*2*sum(f(3:2:n); end fprintf(%3.0f %10.5f %10.5fn,n,I,Iexact-I);end3計算無窮限廣義積分，其精確值。用-10,10替換積分限，則有利用復化梯形公式，調用trapez_v.m文件，取n=10,20的MATLAB程序如下：I=1;a=-10;b=10;fprintf(n h n I1n);n=0;for k=1:2;n=n+20;h=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=exp(-x.2);I=trapez_v(f,h);I1=1/sqrt(pi)*I;fprintf(%3.1f %10.0f %10.8fn

17、,h,n,I1);end4取精度為，分別用和作為計算停止的條件，用Romberg程序計算d，取精度為,并與精確值比較。然后取精度為，用作為計算停止的條件，觀察用龍貝格求積公式計算的結果與精確值的絕對誤差是否滿足。MATLAB程序如下F=inline(1./(1+x); RT,R,wugu,h=romberg(F,0,1.5,1.e-8,13)syms x fi=int(1/(1+x),x,0,1.5); Fs=double(fi), wR=double(abs(fi-R), wR1= wR - wugu其中龍貝格積分的MATLAB程序function RT,R,wugu,h=romberg(f

18、un,a,b, wucha,m)n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4); RT(1,1)=h*(feval(fun,a)+feval(fun,b)/2;while(wuguwucha)&(km)|(k0,dispreturnendif RA=RB if RA=ndisp X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:n m= B(k,p)/ B(p,p); B(k,p

19、:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse dispendend四、實驗要求：1寫出高斯順序消去法、列主元消去法解線性方程組的算法，編寫程序上機調試出結果；2進一步加深對高斯消去法的理解。實驗四 解線性方程組的迭代法一、實驗目的：1熟悉解線性方程組的迭代法的有關理論和方法；2掌握用迭代法求解線性方程組的基本思想和計算步驟；3熟悉逐次超松弛迭代

20、法求解線性方程組的過程；4注意所用方法的收斂性及其收斂速度問題。二、實驗內容：用高斯-塞德爾迭代法、逐次超松弛迭代法求線性方程組的解。三、實驗過程：1用高斯-塞德爾迭代法解下列線性方程組,取初始值，要求當時，迭代終止。（a） （b）用高斯-塞德爾迭代定義解線性方程組Ax=b的MATLAB主程序如下：function X=gsdddy(A,b,X0,P,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); dD=det(D);if dD=0dispelsedispiD=inv(D-L); B2=iD*U;f2=iD*b;jX=Ab; X=X0

21、; n m=size(A);for k=1:max1X1= B2*X+f2; djwcX=norm(X1-X,P);xdwcX=djwcX/(norm(X,P)+eps);if (djwcXwucha)|(xdwcXwucha) return else k,X1,k=k+1;X=X1;endendif (djwcXwucha)|(xdwcXwucha) dispelsedispX=X;jX=jXendendX=X;D,U,L,jX=jX (a)的MATLAB程序如下：A=10 -1 -2;-1 10 -2;-1 -1 0.5; b=7.2;8.3;4.2; X0=0;0;0;X=gsdddy(

22、A,b,X0,inf, 0.001,100)(b)的MATLAB程序如下：A=3 4 -5 7;2 -8 3 -2;4 51 -13 16;7 -2 21 3;b=5;2;-1;21;X0=0;0;0;0;X=gsdddy(A,b,X0,inf,0.001,100)2用逐次超松弛迭代法解由一電路得到的方程組解的精度為0.001。逐次超松弛法的MATLAB程序如下：a(1,1)=1/2+1/4+1/3;a(1,2)=-1/4;a(1,3)=-1/3;a(2,1)=a(1,2);a(2,2)=1/4+1/3+1/5;a(2,3)=-1/5;a(3,1)=a(1,3);a(3,2)=a(2,3);a

23、(3,3)=1/3+1/3+1/5;y(1)=20/2;y(2)=0;y(3)=5/3;x=zeros(1,3);w=1.2;for it=1:50 error=0; for i=1:3 s=0;xb=x(i); for j=1:3 if i=j,s=s+a(i,j)*x(j);end end x(i)=w*(y(i)-s)/a(i,i)+(1-w)*x(i); error=error+abs(x(i)-xb); end fprintf(it.no.=%3.0f,error=%7.2en,it,error) if error/31)&(xdpiancha0.5)&(k3) disp retur

24、n; end if (piancha 0.001)&(xdpiancha3) disp return; endp=(i-1) piancha xdpiancha xk;用不同迭代法求解的MATLAB程序如下：function m=fun1 (x)m=(10-x*x)/2;k,piancha,xdpiancha,xk= diedai1(2,5)運行后輸出用迭代公式的結果k,piancha,xdpiancha,xk=1.000 1.000 0.33333333333333 3.0002.000 2.500 5.000 0.500 3.000 4.375 0.89743589743590 4.875

25、 4.000 11.75781250000000 1.751 -6.882812500000005.000 11.813 0.63167031671297 -18.68655395507813請用戶注意：此迭代序列發散k=5，piancha = 11.813，xdpiancha = 0.63167031671297，xk = -18.68655395507813由以上運行后輸出的迭代序列與根=2.316 624 790 355 40相差越來越大，即迭代序列發散，此迭代法就失敗.這時偏差piancha逐漸增大且偏差的相對誤差xdpiancha的值大于0.5。請重新輸入新的迭代公式function

26、 m=fun1 (x)m=10/(x+2);k,piancha,xdpiancha,xk= diedai1(2,5) 用迭代公式運行后輸出的結果k,piancha,xdpiancha,xk=1.000 0.500 0.200 2.500 2.000 0.27777777777778 0.125 2.222222222222223.000 0.14619883040936 0.173 2.36842105263158 4.000 0.555 0.116 2.289 5.000 0.128 0.115 2.33146067415730k=5,piancha =0.128,xdpiancha = 0

27、.115,xk = 2.33146067415730可見，偏差piancha和偏差的相對誤差xdpiancha的值逐漸變小，且第5次的迭代值xk =2.331 460 674 157 30與根=2.316 624 790 355 40接近，則迭代序列收斂，但收斂速度較慢，此迭代法較為成功。請重新輸入新的迭代公式function m=fun1 (x)m=x-(x*x+2*x-10)/(2*x+2);k,piancha,xdpiancha,xk= diedai1(2,5)用迭代公式運行后輸出的結果k,piancha,xdpiancha,xk=1.000 0.33333333333333 0.142

28、85714285714 2.333333333333332.000 0.667 0.432 2.31666666666667 3.000 0.690 0.822 2.31662479061977 4.000 0.437 0.412 2.316624790355405.000 0 0 2.31662479035540k = 5; piancha = 0; xdpiancha = 0; y = 2.31662479035540.可見，偏差piancha和偏差的相對誤差xdpiancha的值越來越小，且第5次的迭代值xk=2.316 624 790 355 40與根=2.316 624 790 35

29、5 40相差無幾，則迭代序列收斂，且收斂速度很快，此迭代法成功。2取初值，用牛頓法解非線性方程，要求精度為0.00001.牛頓法的MATLAB程序如下function x=Newt_n(f_name,x0)x0=10;x=x0,xb=x-999;n=0;del_x=0.01;while abs(x-xb)0.00001 n=n+1;xb=x; if n300,break;end y=feval(f_name,x);y_driv=(feval(f_name,x+del_x)-y)/del_x; x=xb-y/y_driv; fprintf(n=%3.0f,x=%12.5e,y=%12.5e,n,

30、x,y) fprintf(yd=%12.5en,y_driv)endfprintf(n final answer=%12.5en,x);其中調用的函數f_name(x)定義如下function y=f_name(x)y=x2-155;3用弦截法解由一物理問題抽象出的物理方程 弦截法的MATLAB程序如下function y=zlj(x0,x1)x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0);n=1;while(abs(x1-x0)1.0e-4)&(n=10000000) x0=x1;x1=x2; x2=x1-fc(x1)*(x1-x0)/(fc(x1)-fc(x0);

31、n=n+1; fprintf(n=%3.0fn,n-1) fprintf(fv=%12.5en,x0)endnx2四、實驗要求：掌握解非線性問題的二分法、牛頓法、弦截法，并會用牛頓法、弦截法解非線性方程。實驗六 常微分方程數值解法一、實驗目的：1熟悉求解常微分方程初值問題的有關方法和理論；2熟悉顯式歐拉公式、二龍格-庫塔公式和四階龍格-庫塔公式的使用；3會編制上述方法的計算程序；4通過對各種求解方法的計算實習，體會各種解法的功能、優缺點及適用場合，會選取適當的求解方法。二、實驗內容：1列出常微分方程并利用顯式歐拉公式求其數值解；2利用二龍格-庫塔公式和四階龍格-庫塔公式求微分方程的數值解。三、實驗過程：1某跳傘者在t=0時刻從飛機上跳出，假設初始時刻的垂直速度為0，且跳傘者垂直下落。已知空氣阻力為，其中c為常數，v為垂直速度，向下方向為正。(a) 寫出次跳傘者的速度滿足的微分方程；(b) 若此跳傘者的質量為M=70kg，且已知c=0.27kg/m，利用顯式歐拉公式畫出的速度曲線(取h=0.1s)。(a) 跳傘者速度