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

1、.wd.wd.wd.數(shù)值分析實驗報告第二章實驗題目：分別用二分法、牛頓迭代法、割線法、史蒂芬森迭代法求方程fx=x2+1x-15=0的根x=1，觀察不同初始值下的收斂性，并給出結論。問題分析：題目有以下幾點要求：不同的迭代法計算根，并比較收斂性。選定不同的初始值，比較收斂性。實驗原理：各個迭代法簡述二分法：取有根區(qū)間a,b的重點x0，確定新的有根區(qū)間a1,b1的區(qū)間長度僅為a,b區(qū)間長度的一版。對壓縮了的有根區(qū)間a1,b1重復以上過程，又得到新的有根區(qū)間a2,b2，其區(qū)間長度為a1,b1的一半，如此反復，可得一系列有根區(qū)間，區(qū)間收斂到一個點即為根。牛頓迭代法：不動點迭代法的一種特例，具有局部二

2、次收斂的特性。迭代格式為xn+1=xn-fxnfxn , n=0,1,2, 割線法：是牛頓法的改進，具有超線性收斂的特性，收斂階為1.618. 迭代格式為xn+1=xn-fxnfxn-fxn-1xn-xn-1 , n=1,2, 史蒂芬森迭代法：采用不動點迭代進展預估校正。至少是平方收斂的。迭代格式為yn=xnzn=ynxn+1=xn-(yn-xn)2zn-2yn+xn 這里x可采用牛頓迭代法的迭代函數(shù)。實驗內容：寫出該問題的fx函數(shù)代碼如下：function py= f(x)syms k;y=(k2+1)*(k-1)5;yy=diff(y,k);py(1)=subs(y,k,x);py(2)=

3、subs(yy,k,x);end分別寫出各個迭代法的迭代函數(shù)代碼如下：二分法：function y=dichotomie(a,b,e)i=2;m(1)=a;while abs(a-b)e t=(a+b)/2; s1=f(a); s2=f(b); s3=f(t);if s1(1)*s3(1)=e s=f(x); t=x-s(1)/s(2); en=t-x; x=t; m(i)=t; i=i+1;endy=x,i+1,m;end牛頓割線法：function y=Secant(x1,x2,e)i=3;m(1)=x1,m(2)=x2;while abs(x2-x1)=e s1=f(x1); s2=f(

4、x2); t=x2-(x2-x1)*s2(1)/(s2(1)-s1(1); x1=x2; x2=t; m(i)=t; i=i+1;endy=x2,i+1,m;end史蒂芬森迭代法：Functionp=StephensonIterative(x,e)i=2;m(2)=x;en=2*e;while abs(en)=e y=fai(x); z=fai(y); t=x-(y-x)2/(z-2*y+x); en=t-x; x=t; m(i)=t; i=i+1;endp=x,i+1,m;end因為x經常被使用，故可以寫一個x函數(shù)。代碼如下：function y=fai(x) s=f(x); y=x-s(1

5、)/s(2);end可以繪制不同的圖形來比較不同迭代法的收斂性和不同初值下的收斂性。代碼如下：clear all;%一樣初始值，不同迭代法下的收斂x1=dichotomie(0,3,1e-10);x2=NewtonIterative(0,1e-10);x3=Secant(0,2,1e-10);x4=StephensonIterative(0,1e-10);x1(2),x2(2),x3(2),x4(2)figure,subplot(2,2,1),plot(x1(3:x1(2),title(二分法);subplot(2,2,2),plot(x2(3:x2(2),title(牛頓迭代法);subpl

6、ot(2,2,3),plot(x3(3:x3(2),title(牛頓割線法);subplot(2,2,4),plot(x4(3:x4(2),title(史蒂芬森迭代法);figure,subplot(2,2,1),plot(x1(4:x1(2)-1)-x1(1)./(x1(3:x1(2)-2)-x1(1),title(二分法);subplot(2,2,2),plot(x2(4:x2(2)-1)-x2(1)./(x2(3:x2(2)-2)-x2(1),title(牛頓迭代法);subplot(2,2,3),plot(x3(4:x3(2)-1)-x3(1)./(x3(3:x3(2)-2)-x3(1

7、),title(牛頓割線法);subplot(2,2,4),plot(x4(4:x4(2)-1)-x4(1)./(x4(3:x4(2)-2)-x4(1),title(史蒂芬森迭代法);%不同初始值，一樣迭代法下的收斂性x5=dichotomie(-1,1,1e-10);x6=dichotomie(-2,3,1e-10);x7=dichotomie(0,4,1e-10);x8=dichotomie(-4,4,1e-10);x9=NewtonIterative(-2,1e-10);x10=NewtonIterative(-4,1e-10);x11=NewtonIterative(4,1e-10);

8、x12=NewtonIterative(6,1e-10);figure,subplot(1,2,1),plot(1:x1(2)-2,x1(3:x1(2),1:x5(2)-2,x5(3:x5(2),1:x6(2)-2,x6(3:x6(2),1:x7(2)-2,x7(3:x7(2),1:x8(2)-2,x8(3:x8(2),title(二分法);subplot(1,2,2),plot(1:x2(2)-2,x2(3:x2(2),1:x9(2)-2,x9(3:x9(2),1:x10(2)-2,x10(3:x10(2),1:x11(2)-2,x11(3:x11(2),1:x12(2)-2,x12(3:x

9、12(2),title(牛頓迭代法);x13=Secant(-1,1,1e-10);x14=Secant(-4,5,1e-10);x15=Secant(0,7,1e-10);x16=Secant(-8,2,1e-10);x17=StephensonIterative(-1,1e-10);x18=StephensonIterative(-4,1e-10);x19=StephensonIterative(4,1e-10);x20=StephensonIterative(6,1e-10);figure,subplot(1,2,1),plot(1:x3(2)-2,x3(3:x3(2),1:x13(2)

10、-2,x13(3:x13(2),1:x14(2)-2,x14(3:x14(2),1:x15(2)-2,x15(3:x15(2),1:x16(2)-2,x16(3:x16(2),title(牛頓割線法);subplot(1,2,2),plot(1:x4(2)-2,x4(3:x4(2),1:x17(2)-2,x17(3:x17(2),1:x18(2)-2,x18(3:x18(2),1:x19(2)-2,x19(3:x19(2),1:x20(2)-2,x20(3:x20(2),title(史蒂芬森迭代法);實驗結果：各個迭代值分布不同迭代法下的得到的迭代值迭代值的情況如下：二分法牛頓迭代法牛頓割線法

11、史蒂芬森迭代法00001.50000000000.20000000002.00000000001.35555555560.75000000000.37049180320.33333333330.98161652831.12500000000.50764420760.38071968010.99994600030.93750000000.61461894470.49828334190.99999999951.03125000000.69738690980.57049963330.98437500000.76155380910.63938062441.00781250000.81154111860

12、.69427858790.99609375000.85067638570.74116926531.00195312500.88144821230.78027159970.99902343750.90572974000.8132927871當二分法的初始區(qū)間選為0,3，誤差限為110-10，牛頓迭代法初值選為0，誤差限為110-10，牛頓割線法初始點為0,2，誤差限為110-5，史蒂芬森迭代法初始點選為0，誤差限為110-10，迭代情況如以下列圖。迭代次數(shù)分別為38次，100次，140次，9次。故而，史蒂芬森迭代法速度最快，效果最好。收斂情況不同迭代法下迭代值得收斂情況二分法收斂效果較差，牛頓迭

13、代法和牛頓割線法相近，史蒂芬森迭代法收斂次數(shù)高于1，效果最好不同初值的收斂情況二分法，牛頓迭代法下不同初值的收斂情況牛頓割線法，史蒂芬森迭代法下不同初值的收斂情況二分法的五個初始區(qū)間分別為0,3,-1,1,-2,3,0,4,-4,4；牛頓迭代法的五個初始值分別為0,-2,-4,4,6；牛頓割線法的五個初始區(qū)間分別為0,2,-1,1,-4,5,0,7,-8,2；史蒂芬森迭代法的五個初始值分別為0,-1,-4,4,6；由圖可知，它們最終均到達收斂。收斂性分析及結論：二分法收斂較慢且不能求解崇根，但算法簡單；此處牛頓法具有了平方收斂；從迭代次數(shù)上看，牛頓割線法較牛頓法的多，所以收斂性較差，是超線性收

14、斂；史蒂芬森迭代法收斂效果最好。因為牛頓迭代法是局部的二次收斂，所以要注重初值的選取，本次實驗中選擇的初值均得到了收斂，效果比較好。牛頓割線法也應注意初值的選取。第三章實驗題目：區(qū)間-1,1作等距劃分：xk=-1+kh , h=2n , k=0,1,n以xk為結點對函數(shù)fx=15+x2進展插值逼近。分別取h=1,5,10,20,25用牛頓插值對fx進展逼近，并在同一坐標系下做出函數(shù)的圖形，進展比較。寫出插值函數(shù)對fx的逼近程度與節(jié)點個數(shù)的關系，并分析原因；試用三次樣條插值對fx進展逼近，在同一坐標下畫出圖形，觀察樣條插值函數(shù)對fx的逼近程度與節(jié)點個數(shù)的關系；整體插值有何局限性如何防止一組數(shù)據(jù)如

15、下，求其擬合曲線.數(shù)據(jù)表i012345678910 xi23478101114161819yi106.42108.2109.5110109.93110.49110.59110.6110.76111111.2求以上數(shù)據(jù)形如 yx=c0+c1x+c2x2 的擬合曲線及其平方誤差；求以上數(shù)據(jù)形如 yx=ae-bx 的擬合曲線及其平方誤差；通過觀察12的結果，寫出你對數(shù)據(jù)擬合的認識.問題分析：題目除上述要求之外還有以下幾點：明確整體插值和分段插值的不同。牛頓插值多項式屬于整體插值，三次樣條插值屬于低次分段插值。將結果在同一坐標下繪制出。但是為了方便分析節(jié)點個數(shù)對于插值效果的影響，也可以單獨繪制。第二題

16、中為了確定各個參數(shù)的大小，可以進展適當變換，轉化為線性，運用最小二乘法，得到擬合。實驗原理：牛頓插值多項式：對于給定的插值節(jié)點 ax0 x1xnb , 構造次數(shù)不超過 n 的插值多項式Nnx=a0+a1x-x0+a2x-x0 x-x1+anx-x0 x-x1x-xn-1使其滿足插值條件Pnxi=yi=fxi i=0,1,n這樣得到的插值多項式 Nnx 稱為 Newton 插值多項式。系數(shù) ai 為差商，可以通過構造差商表得到。三次樣條插值：三次樣條插值函數(shù) Sx 在每個小區(qū)間 xk-1,xk 上為三次多項式；在全進 a,b 上存在二階連續(xù)導數(shù)；其次符合插值條件。Matlab 中存在內置的三次樣

17、條插值函數(shù)，命令為 spline .實驗內容：第一題：牛頓插值函數(shù)的構造代碼如下：function f=Newton(x0,y0) %牛頓多項式插值函數(shù)syms x;SZ=size(x0,2);a(1)=y0(1);y(:,1)=y0;for j=2:SZ nx1=1;for i=1:SZ-j+1; nx2=nx1+j-1; y(i,j)=(y(i,j-1)-y(i+1,j-1)/(x0(nx1)-x0(nx2); nx1=nx1+1;endendf=y(1,1);for j=2:SZ ff=y(1,j);for i=1:j-1 ff=ff*(x-x0(i);end f=f+ff;endend

18、牛頓和三次樣條插值情況及比較：代碼如下：clear all;clc;syms x;fx=1/(5+x2);%牛頓多項式插值x0=-1:0.01:1;y0=subs(fx,x0);n1=1,5,10,20,25;n2=5,55,60,62,67;h1=2./n1;h2=2./n2;x1=-1:h1(1):1; y1=subs(fx,x1); f1=Newton(x1,y1); y01=subs(f1,x0); x2=-1:h1(2):1; y2=subs(fx,x2); f2=Newton(x2,y2); y02=subs(f2,x0); x3=-1:h1(3):1; y3=subs(fx,x3

19、); f3=Newton(x3,y3); y03=subs(f3,x0); x4=-1:h1(4):1; y4=subs(fx,x4); f4=Newton(x4,y4); y04=subs(f4,x0); x5=-1:h1(5):1; y5=subs(fx,x5); f5=Newton(x5,y5); y05=subs(f5,x0); figure,plot(x0,y0,x0,y01,x0,y02,x0,y03,x0,y04,x0,y05),title(所有結果);x6=-1:h2(1):1; y6=subs(fx,x6); f6=Newton(x6,y6); y06=subs(f6,x0)

20、; x7=-1:h2(2):1; y7=subs(fx,x7); f7=Newton(x7,y7); y07=subs(f7,x0); x8=-1:h2(3):1; y8=subs(fx,x8); f8=Newton(x8,y8); y08=subs(f8,x0); x9=-1:h2(4):1; y9=subs(fx,x9); f9=Newton(x9,y9); y09=subs(f9,x0); x10=-1:h2(5):1; y10=subs(fx,x10); f10=Newton(x10,y10); y010=subs(f10,x0); figure,plot(x0,y0,x0,y06,x

21、0,y07,x0,y08,x0,y09,x0,y010),title(龍格現(xiàn)象);%三次樣條插值spline自帶命令x0=-5:0.01:5;y0=subs(fx,x0);figure,subplot(2,3,1),plot(x0,y0),title(準確圖);y11=spline(x1,y1,x0);subplot(2,3,2),plot(x0,y11),title(n=1);y12=spline(x2,y2,x0);subplot(2,3,3),plot(x0,y12),title(n=5);y13=spline(x3,y3,x0);subplot(2,3,4),plot(x0,y13),

22、title(n=10);y14=spline(x4,y4,x0);subplot(2,3,5),plot(x0,y14),title(n=20);y15=spline(x5,y5,x0);subplot(2,3,6),plot(x0,y15),title(n=25);figure,plot(x0,y0,x0,y11,x0,y12,x0,y13,x0,y14,x0,y15),title(所有結果);第二題：代碼如下：clear all;clc;x=2 3 4 7 8 10 11 14 16 18 19;y=106.42 108.2 109.5 110 109.93 110.49 110.59 1

23、10.6 110.76 111 111.2;%Agenda 1A=(x.*x) x ones(11,1);A1=A*A;y1=A*y;c1=A1y1x1=2:0.1:19;y1=polyval(c1,x1);plot(x,y,k*,x1,y1,r),gtext(曲線一),hold ony01=polyval(c1,x);delta1_2=sum(y-y01).2)%Agenda 2xx=-1./x;yy=log(y);B=ones(11,1) xx;B1=B*B;yy1=B*yy;c2=B1yy1;syms s;a=exp(c2(1);b=c2(2);a bf=a*exp(-b/s);x2=x

24、1;y2=subs(f,x2);plot(x2,y2,b),gtext(曲線二);y02=subs(f,x);delta2_2=sum(y-y02).2)實驗結果：第一題：牛頓插值結果不同節(jié)點數(shù)的牛頓插值多項式綜合圖龍格現(xiàn)象由圖可知，三次樣條插值結果不同節(jié)點下的三次樣條插值結果不同節(jié)點數(shù)的三次樣條的綜合圖由圖可知，牛頓插值多項式在n較小的時候如5,10,20,25差值效果良好，當n變大時如60,62,65,67,100,200就出現(xiàn)了龍格現(xiàn)象，三次樣條在各個子區(qū)間內為三次多項式，擬合效果好.第二題第一問的多項式擬合得到的擬合曲線為yx=106.2927+0.6264x-0.0205x2平方誤差

25、為1=2.7796第二問的擬合曲線為yx=111.4940e-0.0903x平方誤差為2=0.4719擬合曲線如以下列圖擬合情況從圖中可以看出曲線二大體符合黑點的分布情況，擬合效果較好。結論：整體插值隨著節(jié)點個數(shù)的增多，多項式的次數(shù)也在升高。高次多項式的插值會出現(xiàn)龍格現(xiàn)象，震蕩劇烈，逼近效果不理想。但是當節(jié)點很多時，低次插值的精度又不夠，所以為了防止這一局限性采用分段低次插值。其中三次樣條插值有良好的收斂性和光滑性，效果較好。數(shù)據(jù)擬合時，可以選擇趨勢相似的函數(shù)形式，在求解相關參量時，將函數(shù)進展適當變換，從而運用最小二乘法得到擬合結果。第四章實驗題目：設計區(qū)間分半求積算法、龍貝格求積算法和自適應

26、辛普森求積算法的程序，觀察n=1,10,100,500 時，積分-11x2cos(nx)dx 的結果，并給出相應的評價. 問題分析：分半求積算法即為復合梯度求積。因為 n 越大，fx=x2cos(nx)在-1,1 內震蕩地越嚴重，自適應辛普森求積算法很難適用，計算復雜度會很高。所以選取辛普森求積算法代替。實驗原理：分半求積算法：將區(qū)間 a,b 等分為 n 個小區(qū)間，分點為xk=a+kh h=b-an , k=0,1,n利用定積分性質，有abfxdx=k=0n-1xkxk+1f(x)dx每個小區(qū)間上均用梯形求積公式，有abfxdx=h2k=0n-1f(xk)+f(xk+1)+k=0n-1-h21

27、2f(k)令Tn=h2k=0n-1f(xk)+f(xk+1)=h2f(a)+f(b)+hk=0n-1f(xk)即為復合梯形公式。龍貝格求積算法：將步長為h的復合梯形公式進展査爾遜外推，就得到龍貝格求積公式。自適應辛普森求積算法：按照被奇函數(shù)在區(qū)間上的變化來安排求積結點的辛普森算法。實驗內容：寫得各個求積算法的函數(shù)。代碼如下：1分半求積算法：function f=CompositeTrapezoida(n)syms x;y=x*x*cos(n*x);m=1:100;a=-1;b=1;for i=1:100 h=(b-a)/m(i); k=0:m(i); xk=a+k*h; yk=subs(y,x

28、k); f(i)=0;for j=1:m(i)f(i)=f(i)+(yk(j)+yk(j+1);end f(i)=h/2*f(i);endend2辛普森算法function f=CompositeSimpson(n)syms x;y=x*x*cos(n*x);m=1:100;a=-1;b=1;for i=1:100 h=(b-a)/(2*m(i); k=0:(2*m(i); xk=a+k*h; yk=subs(y,xk); f(i)=0;for j=1:m(i)f(i)=f(i)+(yk(2*j-1)+4*yk(2*j)+yk(2*j+1);end f(i)=h/3*f(i);endm;f;e

29、nd4龍貝格求積算法function f=Romberg(n)epsilon=0.0001;syms x;y=x*x*cos(n*x);a=-1;b=1;h=b-a;fa=subs(y,a);fb=subs(y,b);T(1)=h/2*(fa+fb);m=1;while 1 h=h/2; k=1:(2m-1); xk=a+k*h; yk=subs(y,xk); su=sum(yk); S(1)=h/2*(fa+fb)+h*su;for i=1:mS(i+1)=S(i)+(S(i)-T(i)/(4i-1);endif(abs(S(m+1)-T(m)epsilon) break;end T=S;

30、m=m+1;endf=S(m+1);end5自適應辛普森求積算法function f=AdaptiveSimpson(n)epsilon=0.001;syms x;y=x*x*cos(n*x);a=-1;b=1;h=b-a;p=a b;p0=p;ep=epsilon;m=0;q=0;f=0;while 1 n1=length(ep); n=length(p0);if(n=1) break; end h=p0(2)-p0(1); k=0:4; yk=subs(y,p0(1)+k*h/4); s0=h/6*(yk(1)+4*yk(3)+yk(5);s1=h/12*(yk(1)+4*yk(2)+yk

31、(3);s2=h/12*(yk(3)+4*yk(4)+yk(5);if(abs(s0-s1-s2)=2 ep=ep(2:n1);end q=q+1;else m=m+1; p0=yk(1),yk(3),p0(2:n);if n1=1 ep=ep(1)/2 ep(1)/2;else ep=ep(1)/2 ep(1)/2 ep(2:n1);endif q=0 p=p0;else p=p(1:q) p0;endendendend各個求積公式下的積分結果代碼如下:clear all;clc;%復合梯形法y11=CompositeTrapezoida(1);y12=CompositeTrapezoida

32、(10);y13=CompositeTrapezoida(100);y14=CompositeTrapezoida(500);1:100;y11;y12;y13;y14%復合辛普森法y21=CompositeSimpson(1);y22=CompositeSimpson(10);y23=CompositeSimpson(100);y24=CompositeSimpson(500);2:2:200;y21;y22;y23;y24%龍貝格法y31=Romberg(1);y32=Romberg(10);y33=Romberg(100);y34=Romberg(500);y31 y32 y33 y34

33、%自適應辛普森法y41=AdaptiveSimpson(1)%以下復雜度太高 不出結果%y42=AdaptiveSimpson(10); %y43=AdaptiveSimpson(100);%y44=AdaptiveSimpson(500);%y41 y42 y43 y44實驗結果：不同積分法得到的結果n110100500復合梯形法0.4782831994-0.1399401910-0.00601370080.0023823806復合辛普森法0.4782672530-0.1401909997-0.00985113540.0072247738龍貝格法0.4782672574-0.14019099

34、850.6112212333-0.2492135619結果評價：龍貝格法的精度高，最接近準確解。因為被積分函數(shù)為fx=x2cos(nx)隨著n的增大，震蕩越劇烈。所以積分面積越來越小。第五章實驗題目：分別取步長 h=0.5,0.1,0.05,0.01, 用龍格-庫塔法求解y=ty13 , y1=1 , 計算到t=2 , 并與準確解yt=(t2+2)332 比較，觀察龍格-庫塔法的收斂性. 分別取步長 h=0.1,0.05,0.01,0.005, 用顯示歐拉法和隱式歐拉法求解y=-50y , y0=100 , 由結果分析算法的穩(wěn)定性. 選擇某常微分方程初值問題的數(shù)值方法計算01sinttdt 的

35、近似值，并保證有4位有效數(shù)字. 問題分析：龍格-庫塔有很多的格式，第一題以常用的四級四階龍格庫塔方法為例。假設，ft=sintt，由于01sinttdt=f(1)-f0第三題可以轉化為求解常微分方程ft=sintt的問題。實驗原理：四級四階龍格-庫塔格式y(tǒng)i+1=yi+h6(k1+2k2+2k3+k4)k1=f(ti,yi)k2=f(ti+h2,yi+h2k1)k3=f(ti+h2,yi+h2k2)k4=f(ti+h,yi+hk3)實驗內容：求解第一題代碼如下：clear all;clc%以四級四階龍格庫塔格式為例t0=1;h=0.5 0.1 0.05 0.01;n=1./h+1;y1(1)=

36、1;y2(1)=1;y3(1)=1;y4(1)=1;t1=1:0.5:2;t2=1:0.1:2;t3=1:0.05:2;t4=1:0.01:2;for j=2:n(1) k1=t1(j-1)*y1(j-1)(1/3); k2=(t1(j-1)+h(1)/2)*(y1(j-1)+(h(1)/2)*k1)(1/3); k3=(t1(j-1)+h(1)/2)*(y1(j-1)+(h(1)/2)*k2)(1/3); k4=t1(j)*(y1(j-1)+h(1)*k3)(1/3); y1(j)=y1(j-1)+h(1)/6*(k1+2*k2+2*k3+k4);endfor j=2:n(2) k1=t2(

37、j-1)*y2(j-1)(1/3); k2=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k1)(1/3); k3=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k2)(1/3); k4=t2(j)*(y2(j-1)+h(2)*k3)(1/3); y2(j)=y2(j-1)+h(2)/6*(k1+2*k2+2*k3+k4);endfor j=2:n(3) k1=t3(j-1)*y3(j-1)(1/3); k2=(t3(j-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k1)(1/3); k3=(t3(j-1)+h(3)/2)*(y3(

38、j-1)+(h(3)/2)*k2)(1/3); k4=t3(j)*(y3(j-1)+h(3)*k3)(1/3); y3(j)=y3(j-1)+h(3)/6*(k1+2*k2+2*k3+k4);endfor j=2:n(4) k1=t4(j-1)*y4(j-1)(1/3); k2=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k1)(1/3); k3=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k2)(1/3); k4=t4(j)*(y4(j-1)+h(4)*k3)(1/3); y4(j)=y4(j-1)+h(4)/6*(k1+2*k2+2*k3

39、+k4);endx=1:0.01:2;syms tyy=(t2+2)/3)(3/2);yy1=subs(yy,x);figure,subplot(2,3,1),plot(x,yy1),title(準確值);subplot(2,3,2),plot(x,yy1,t1,y1),title(h=0.5);subplot(2,3,3),plot(x,yy1,t2,y2),title(h=0.1;subplot(2,3,4),plot(x,yy1,t3,y3),title(h=0.05);subplot(2,3,5),plot(x,yy1,t4,y4),title(h=0.01);subplot(2,3,

40、6),plot(x,yy1,t1,y1,t2,y2,t3,y3,t4,y4),title(所有結果);求解第二題代碼如下：clear all;t0=0;a=0.25; %通過更改a，更改取值區(qū)間syms tyy=exp(-50*t);x=0:0.001:a;yy1=subs(yy,x);h=0.1 0.05 0.01 0.005;n=a./h+1;t1=0:0.1:a;t2=0:0.05:a;t3=0:0.01:a;t4=0:0.005:a;%顯式歐拉法y1(1)=100;y2(1)=100;y3(1)=100;y4(1)=100;for j=2:n(1) y1(j)=y1(j-1)+h(1)

41、*(-50)*y1(j-1);endfor j=2:n(2) y2(j)=y2(j-1)+h(2)*(-50)*y2(j-1);endfor j=2:n(3) y3(j)=y3(j-1)+h(3)*(-50)*y3(j-1);endfor j=2:n(4) y4(j)=y4(j-1)+h(4)*(-50)*y4(j-1);endfigure,plot(x,yy1,t1,y1,t2,y2,t3,y3,t4,y4),.gtext(h=0.1),gtext(h=0.05),gtext(h=0.01),gtext(h=0.005);%隱式歐拉法y5(1)=100;y6(1)=100;y7(1)=100

42、;y8(1)=100;for j=2:n(1) y5(j)=y5(j-1)/(50*h(1)+1);endfor j=2:n(2) y6(j)=y6(j-1)/(50*h(2)+1);endfor j=2:n(3) y7(j)=y7(j-1)/(50*h(3)+1);endfor j=2:n(4) y8(j)=y8(j-1)/(50*h(4)+1);endfigure,plot(x,yy1,t1,y5,t2,y6,t3,y7,t4,y8),.gtext(h=0.1),gtext(h=0.05),gtext(h=0.01),gtext(h=0.005);%一樣步長下的穩(wěn)定性比較figuresub

43、plot(2,2,1),plot(x,yy1,k,t1,y1,b,t1,y5,r),title(h=0.1);subplot(2,2,2),plot(x,yy1,k,t2,y2,b,t2,y6,r),title(h=0.05);subplot(2,2,3),plot(x,yy1,k,t3,y3,b,t3,y7,r),title(h=0.01);subplot(2,2,4),plot(x,yy1,k,t4,y4,b,t4,y8,r),title(h=0.005);求解第三題代碼如下：%防止t=0作分母，所以用隱式歐拉法clear all;t0=0;h=0.0005 0.0001 0.00005;

44、n=1./h+1;y1(1)=5;y2(1)=5;y3(1)=5;y4(1)=5;t1=0:0.0005:1;t2=0:0.0001:1;t3=0:0.00005:1;for j=2:n(1) y1(j)=y1(j-1)+h(1)*sin(t1(j)/t1(j);endfor j=2:n(2) y2(j)=y2(j-1)+h(2)*sin(t2(j)/t2(j);endfor j=2:n(3) y3(j)=y3(j-1)+h(3)*sin(t3(j)/t3(j);endformat longy(1)=y1(n(1)-y1(1);y(2)=y2(n(2)-y2(1);y(3)=y3(n(3)-y

45、3(1);y實驗結果：第一題結果不同步長下的四級龍格-庫塔的解顯示歐拉法求解結果不同步長下的顯式歐拉法的解隱式歐拉法求解結果不同步長下的隱式歐拉法的解顯式與隱式的穩(wěn)定性比較結果不同步長下顯式與隱式的穩(wěn)定性比較第三題結果不同步長下的積分結果步長 h0.00050.00010.00005積分結果0.9460434318390180.9460751436654500.946079107079003題目要求保證有4位有效數(shù)字，所以所求結果為：0.9461結果分析：四級四階龍格庫塔法具有大范圍的收斂性。四級四階龍格庫塔法每計算一步需要計算四次的值，計算有一定的復雜性。當步長h為0.1和0.05時顯式歐拉

46、法不具有穩(wěn)定性。隱式歐拉法在各個步長下的穩(wěn)定情況相似，均有較好的穩(wěn)定性。第六章實驗題目：使用列主元高斯消去法解希爾伯特矩陣 H=(hij)nn 方程組 Hx=b , 考察給定的 n=5,10,20,50,100 及相應的 b 時結果的變化，分析其中的原因并給出結論，其中hij=1/(i+j-1)問題分析：希爾伯特矩陣是病態(tài)的，b的微小變化會造成解發(fā)生很大變化?？坍嬀仃嚥B(tài)程度的量為條件數(shù)。條件數(shù)越大，矩陣性能越差?？梢杂嬎銞l件數(shù)來分析解變化很大的現(xiàn)象。實驗原理：高斯列主元素消元法：將矩陣H按列取絕對值最大項進展行行調換，將其以下局部消為0，不斷循環(huán)直到H化為上三角的過程。實驗內容：高斯列主元素

47、消元法函數(shù)代碼如下：function A b=GaussElimination(A,b)n=size(A,1);for i=1:n-1 q=max(abs(A(i:n,i);for j=i:nif(abs(A(j,i)=q) k=j;endend u=A(i,i:n);v=b(i); A(i,i:n)=A(k,i:n);b(i)=b(k); A(k,i:n)=u;b(k)=v;for l=i+1:n A(l,i)=A(l,i)/A(i,i); b(l)=b(l)-A(l,i)*b(i);for s=i+1:n A(l,s)=A(l,s)-A(l,i)*A(i,s);endendendend計算

48、不同階數(shù)下的希爾伯特矩陣的線性方程的解，并通過更改b，對其做微小調整，觀察解的變化代碼如下：clear all;clc;n=5 10 20 50 100;%5階a=1;H=zeros(n(a);for i=1:n(a) for j=1:n(a) H(i,j)=1/(i+j-1);endendb=ones(1,n(a);b1=b+0.0001;HH,bb=GaussElimination(H,b);HH1,bb1=GaussElimination(H,b1);for i=n(a):-1:1 s=0;s1=0;for j=i+1:n(a) s=s+HH(i,j)*x1(1,j);s1=s1+HH1

49、(i,j)*xx1(1,j);end x1(1,i)=(bb(1,i)-s)/HH(i,i);xx1(1,i)=(bb1(1,i)-s1)/HH1(i,i);end%10階a=2;H=zeros(n(a);for i=1:n(a) for j=1:n(a) H(i,j)=1/(i+j-1);endendb=ones(1,n(a);b1=b+0.0001;HH,bb=GaussElimination(H,b);HH1,bb1=GaussElimination(H,b1);for i=n(a):-1:1 s=0;s1=0;for j=i+1:n(a) s=s+HH(i,j)*x2(1,j);s1=

50、s1+HH1(i,j)*xx2(1,j);end x2(1,i)=(bb(1,i)-s)/HH(i,i);xx2(1,i)=(bb1(1,i)-s1)/HH1(i,i);end%20階a=3;H=zeros(n(a);for i=1:n(a) for j=1:n(a) H(i,j)=1/(i+j-1);endendb=ones(1,n(a);b1=b+0.0001;HH,bb=GaussElimination(H,b);HH1,bb1=GaussElimination(H,b1);for i=n(a):-1:1 s=0;s1=0;for j=i+1:n(a) s=s+HH(i,j)*x3(1,j);s1=s1+HH1(i,j)*xx3(1,j);end x3(1,i)=(bb(1,i