數(shù)值分析計算實習題

1、插值法1.下列數(shù)據(jù)點的插值X 0 1 4 9 16 25 36 49 64y012345678可以得到平方根函數(shù)的近似，在區(qū)間0, 64上作圖.(1)用這9個點作8次多項式插值Ls (x).(2)用三次樣條(第一邊界條件)程序求S(x).從得到結(jié)果看在0, 64上，哪個插值更精確；在區(qū)間0,1上，兩 種插值哪個更精確解：(1)拉格朗日插值多項式，求解程序如下syms x 1;xl=0 1 4 9 16 25 36 49 64; yl=0 1 2 3 4 5 6 7 8; n=length(xl);Ls=sym(0); for i=l:n l=sym(yl (i); for k=l:i-l1=1

2、*(x-xl(k)/(xl(i)-xl(k); endfor k=i+l:n 1=1*(x-xl(k)/(xl(i)-xl (k); endLs=Ls+l; end Ls=simplify (Ls) %為所求插值多項式Ls (x).輸出結(jié)果為Ls = -/*xA2+95549/72072*x-l/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6（2）三次樣條插值，程序如下xl=0 1 4 9 16 25 36 49 64;yl=0 1 2 3 4 5 6 7 8;x2二0:1:64;y3=s plin e (

3、xl, yl, x2);p=po Iyfit(x2,y3,3)1 %得到三次樣條擬合函數(shù)S=p(l)+p(2)*x+ p(3)*x-2+p(4)*xA3 % 得到 S(x)輸出結(jié)果為: S = /6464-2399/88*x+/1984*xA2+2656867/624*xA3 在區(qū)間0, 64上，分別對這兩種插值和標準函數(shù)作圖,Plot (x2, sqrt (x2), ' b'，x2, y2, 'x2, y3, ' y')藍色曲線為y=X函數(shù)曲線，紅色曲線為拉格朗日插值函數(shù)曲線，黃色曲線為三次樣條插值曲線10080604020-2010203040506

4、070可以看到藍色曲線與黃色曲線幾乎重合，因此在區(qū)間0, 64上三次樣條插值更精確。在0,1區(qū)間上由上圖看不出差別，不妨代入幾組數(shù)據(jù)進行比較，取 x4=0: :1準確值拉格 朗日插值三 次樣條插值x4=0:1;sqrt (x4)subs (Ls, ' x' , x4)spline (xl,yl,x4)運行結(jié)果為ans 二0ans 二0ans 二 0從這幾組數(shù)值上可以看出在0,1區(qū)間上，拉格朗日插值更精確。數(shù)據(jù)擬合和最佳平方逼近2.有實驗給出數(shù)據(jù)表試求3次、4次多項式的曲線擬合，再根據(jù)數(shù)據(jù)曲線形狀，求一個另外函數(shù)的擬合曲線，用圖示數(shù)據(jù)曲線及相 應的三種擬合曲線。解：(1)三次擬合

5、，程序如下sym x;x0=;y0= ; cc=polyfit (xO,yO,3);S3=cc(l)+cc Q*x+cc 3*xA2+cc (4) *x3 %三次擬合多項式xx=x0(1):x0(length(xO);yy=polyval (cc, xx); plot(xx,yy,'-，);hold on; plot (xO, yO, ' x') ; xlabel (' x*) ; ylabel('y');運行結(jié)果S3 二V624+437/28*x-0524945/28*x2+409/*x3圖像如下2. 51. 50. 50 L10. 10.20

6、.3WMM «0:10.40.50. 60.70.80.91(2)4次多項式擬合sym x;x0=;y0=;cc=polyfit (xO,yO,4);S3=cc(l)+cc Q*x+cc(3)*xA2+cc (4) *x-3+cc *XA4xx=x0(l) : :x0(le ngth(xO);yy=po lyval(cc, xx);plot (xx, yy,' r')；hold on;plot (xO, yO, ' x');xlabel (' x');ylabel (' y');運行結(jié)果S3 =00396215/624-

7、347/28*x+/28*xA2-/624*xA3+849/9007*xA4圖像如下2. 51. 50. 500. 10. 20.30.40. 50.60.70. 80.9(3)另一個擬合曲線，新建一個M-file程序如下：function C,L=lagra n( x, y) w=le ngth(x);n二wT;L二zeros (w, w);for k=l: n+1v=l;for j=l： n+1if r=jV=conv(V, p oly(x(j)/(x(k)-x(j) ; end endL(k, :)=V;endC=y*L在命令窗口中輸入以下的命令：X二口;y二口；cc=po lyfit

8、(x,y,4);xx=x(l) : :x(le ngth(x);yy=po lyval(cc,xx);P lot (xx, yy, T);hold on;plot (x, y, ' x');xlabel (' x');ylabel (' y');x=:y< ; %y中的值是根據(jù)上面兩種擬合曲線而得到的估計數(shù)據(jù)，不是真實數(shù)據(jù)C, L=lagra n(x, y);xx=0:;yy=po lyval(C, xx);hold on;plot (xx, yy, ' b'，x, y,'.')；圖像如下解線性方程組的直接解

9、法3 .線性方程組Ax=b的A及b為A=rlO 7 87 568 6 10 95 9 loj715r3223,b=331則解X= 1用MATLA內(nèi)部函數(shù)求detA及A的所有特征值和cond ( A) 2.若令A+ S A=rW 7 8. 17. 08 5. 048 5, 986 . 99 57 216 589 99. 98，求解(A+S A) ( x+ S x)=b,輸出向量xAx=b解得相對誤差I(lǐng)I SXIJ 2/2及A的相對誤差H S A2/ I|A|和II S刈2.從理論結(jié)果和實際計算兩方面分析線性方程組的關(guān)系.解：(1)程序如下clear;A=10 7 8 7;7 5 6 5;8 6

10、10 9;7 5 9 10;det (A)con d(A, 2)eig(A)輸出結(jié)果為ans 二(2)程序如下A=10 7 ; 6 5;8 9; 5 9 ;b=32 23 33 31'x0=l 111'x=Ab%擾動后方程組的解xl=x-x0Sx的值norm(xl,2)% Sx 的 2-范數(shù)運行結(jié)果為xl 二 ans 二程序如下A=10 7 ; 6 5;8 9; 5 9 ;A0=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10; b=32 23 33 31' x0=l 111'x=Ab; xl=x-x0; norm(xl, 2);Al=A-A

11、0 ;%5 A 的值norm(xl, 2)/norm(x0, 2)%II 5 x| j 2/| xj |2 的值norm (Al, 2)/norm (AO, 2)輸 *II 5 A | 2/| A | 2 的值出結(jié)果為ans =該方程可見A相對誤差只為， 而得到的結(jié)果x的相對誤差就達到了，組是病態(tài)的，A的條件數(shù)為遠遠大于1,當A只有很小的誤差就會給結(jié)果帶來很大的影響。非線性方程數(shù)值解法4 .求下列方程的實根(1) X 八 2-3x+2-e 八 x=0;(2) X 八 3+2x 八 2+10x-20=0.要求：(1)設(shè)計一種不動點迭代法，要使迭代序列收斂，然后再用斯特芬森加速迭代，計算到|x(k

12、)-x(k-1) |vlO/l(-8)為止。(2)用牛頓迭代，同樣計算到|x(k)-x(kT)|vlO>(-8) o輸出迭代初值及各次迭代值和迭代次數(shù)k,比較方法的優(yōu)劣。解：(1)先用畫圖的方法估計根的范圍ezp lot (' xA2-3*x+2-ex p (x)');grid on;2可以估計到方程的根在區(qū)間(0, 1);選取迭代初值為x0=；構(gòu)造不動點迭代公式 x(k+l) = ( x(k) 2+2-e x(k)/3;9 (x)= (2+2-e ) x)/3;當0VXV1時，07 (x)<l;9' (x)v'1;因此迭代公式收斂。程序如下 :fo

13、rmat long;f=i niin eC(x 2+2-ex p(x)/3, ) disp('x=')；x=fevai (f,;disp(x);Eps=lE-8;i=l；whiie 1xO=x;i=i+l; x=fevai (f,x);disp(x);if abs(x-xO)<Eps break;endend i, x運行結(jié)果為Iniine function: f(x)=(x 2+2-ex p(x)/31 二14斯特芬森加速法，程序如下format long;expf=i nii ne('x-(xA2+2-ex p(x)/3-x) A2/(xA2+2-ex p(x

14、)/3)A2+2-ex p(儀人2+2- (x)/3)/3-2*(xA2+2-ex p (x)/3+x)');disp(' x=');x=fevai(f,;disp (x);Eps=lE-8;i=l；whiie 1xO=x; i=i+1; x=fevai (f, x); disp(x);if abs (x-xO)<Eps break;endendi, x運行結(jié)果為x=牛頓迭代法，程序如下：format long;x=sym(，x');f=sym(，xA2-3*x+2-ex p(x)'); df=diff (f, x);FX=x-f/df;Fx=in

15、line (FX);disp(' x=');xl=;disp(xl);Eps=lE-8;i=0;while 1xO=xl;i=i+l; xl=feval(Fx, xl); disp(xl);if abs(xl-xO)<Epsbreak;endendi, xl運行結(jié)果如下：x= i =4 xl(2)先用畫圖的方法估計根的范圍根大約在區(qū)間(1, 2);選取初值xO=;構(gòu)造不動點迭代公式 x (k+l) = (2x (kF2-10x (k) +20F1/3;9 (x) =(-2x/< 2-10x+20F1/3;程序如下:format lo ng;f=i niin e(&g

16、t; (-2*x 2- 10*x+20) 1/3') dis p('x=');x=feval(f,;dis p(x)Ep s=IE-8;while 1xO=x;i=i+l; x=feval(f, x); dis p(x);if abs(x-xO)>1E1O break;endif abs(x-xO)<Eps break;end end i, x運行結(jié)果：Inline function:f(x) = (-2*xA2-10*x+20)Al/3 x=+002+005+0111 =6 x =+011迭代6次后x的值大得令人吃驚，表明構(gòu)造的式子并不收斂.也無法構(gòu)造出收

17、斂的不動點公式牛頓迭代法，程序如下：format long;x=sym(，x');f=sym(，x*3+2*x 2+10*x-20,);df=diff(f, x);FX=x-f/df;Fx=inline(FX);disp(，x=');xl=;disp(xl);Eps=lE-8;i=0;while 1x0=xl;i=i+l; xl=feval(Fx, xl); disp(xl);if abs(xl-xO)<Eps break;endendi, xl運行結(jié)果:X 二xl 二比較三種方法，牛頓法的收斂性比較好，相比不動點迭代法要構(gòu)造出收斂的公式比較難，牛頓法迭代次數(shù)也較少，收斂

18、速度快，只是對初值的要求很高，幾種方法各有利弊，具體采用哪種也需因題而異。常微分方程初值問題數(shù)值解法5.給定初值問題y' =-50y+50x N 2+2x, b W X W 1 ；y(O)=i/3；用經(jīng)典的四階R-K方法解該問題，步長分別取h=,，計算并打印x=(i=O, 1,，10)各點的值，與準確值 y(x)=l/3eA(-50x)+x )1 2 比 較。解：取步長h二，程序如下：憾典的四階R-K方法 clear;F=, -50*y+50*xA2+2*x,;a=0;b=l;h=; n= (b-a) /; X=a:b;Y=zeros (1, n+1);Y(l)=l/3; for i=l: n x=X(i); y=Y(i)； Kl=h*eval(F); x=x+h/2; 尸y+Kl/2;K2=h*eval(F);y=Y(i)+K2/2;K3=h*eval(F); x=X(i) +h; y=Y(i)+K3;K4=h*eval(F);Y(i+1) =Y(i) + (K1+2*K2+2*K3+K4) /6;end先準確值te mp=;f=dsolve (' Dy=-50*y+50*x,