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

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

2、收斂的特性。迭代 格式為心)n = 0*1-2-割線法：是牛頓法的改進，具有超線性收斂的特性，收斂階為1.618.迭代格式為xn + 1心)f (暫)-心_ J（叫- S -1）史蒂芬森迭代法：采用不動點迭代進行預(yù)估校正。至少是平方收斂的。迭代格式為y=t,i+1,m;end牛頓迭代法：fun ctio ny=Newt on lterative(x,e)i=2;en=2*e;m(1)=x;while abs(e n)=es=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牛頓割線法：% =甲（yJ兒-兀這里可采用牛頓迭代法的

3、迭代函數(shù)實驗內(nèi)容：1. 寫出該問題的,|函數(shù)代碼如下：fun ctio npy= f(x)syms k;y=(kA2+1)*(k-1)A5;yy=diff(y,k);py(1)=subs(y,k,x);py(2)=subs(yy,k,x);end2. 分別寫出各個迭代法的迭代函數(shù) 代碼如下：二分法：function y=dichotomie(a,b,e)i=2;m(1)=a;while abs(a-b)et=(a+b)/2;s1=f(a);s2=f(b);s3=f(t);if s1(1)*s3(1)v=0b=t;elsea=t;endm(i)=t;i=i+1;endfunction y=Sec

4、a nt(x1,x2,e) i=3;m(1)=x1,m(2)=x2;while abs(x2-x1)=es1=f(x1);s2=f(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史蒂芬森迭代法：Fun cti onp=Stephe nson lterative(x,e) i=2;m(2)=x;en=2*e;whileabs(e n)=ey=fai(x);z=fai(y);t=x-(y-x)A2/(z-2*y+x);en=t-x;x=t;m(i)=t;i=i+1;end p=x,i+1,m

5、; end3.因為經(jīng)常被使用，故可以寫一個函數(shù)代碼如下：fun ctio ny=fai(x)s=f(x); y=x-s(1)/s(2);end4.可以繪制不同的圖形來比較不同迭代法的收斂性和不同初值下的收斂性 代碼如下：clear all ;%相同初始值，不同迭代法下的收斂x1=dichotomie(0,3,1e-10); x2=Newto nIterative(0,1e-10);x3=Seca nt(0,2, 1e-10);x4=Stephe nso nl terative(0,1e-10);x1(2) ,x2 (2),x3 (2),x4 (2) figure,subplot(2,2,1),

6、plot(x1(3:x1(2),title(subplot(2,2,2),plot(x2(3:x2(2),title(二分法）；牛頓迭代法）；subplot(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(

7、1),title(牛頓迭代法)；subplot(2,2,3),plot(x3(4:x3(2)-1)-x3(1)./(x3(3:x3(2)-2)-x3(1),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=Newto nlt

8、erative(-2,1e-10);x10=Newton Iterative(-4,1e-10);x1 仁Newt on Iterative(4,1e-10);x12=Newton Iterative(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

9、)-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:x12(2),title(牛頓迭代法);x13=Seca nt(-1,1, 1e-10);x14=Seca nt(-4,5, 1e-10);x15=Seca nt(0,7, 1e-10);x16=Seca nt(-8,2, 1e-10);x17=Stephe nson Iterative(-1,1e-10);x18=Stephe nson Iterative(-4,1e-10);x19=Stephe nso nI terative(4

10、,1e-10);x20=Stephe nso nlterative(6,1e-10);figure,subplot(1,2,1),Plot(1:x3(2)-2,x3(3:x3(2),1:x13(2)-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,

11、x19(3:x19(2),1:x20(2)-2,x20(3:x20(2),title(蒂芬森迭代法);實驗結(jié)果：1. 各個迭代值分布牛頓迭代法0 k1050100史蒂芬森迭代法1.5 I 1 I0.5 *ffIr f0 EEECE02468圖1.1不同迭代法下的得到的迭代值迭代值的情況如下:二分法牛頓迭代法牛頓割線法史蒂芬森迭代法00001.50000000000.20000000002.35555555560.75000000000.37049180320.33333333330.98161652831.12500000000.50764420760.38071968010.99994600

12、030.93750000000.61461894470.49828334190.99999999951.69738690980.57049963330.98437500000.76155380910.63938062441.81154111860.69427858790.99609375000.85067638570.74116926531.88144821230.78027159970.99902343750.90572974000.8132927871當二分法的初始區(qū)間選為W”引,誤差限為1 X 10 ,牛頓迭代法初值選為 ,誤差限為1X1O-W,牛頓割線法初始點為0”2,誤差限為lxHT

13、”，史蒂芬森 迭代法初始點選為，誤差限為1x10 1,1，迭代情況如圖所示。迭代次數(shù)分別為 38次，100次，140次，9次。故而，史蒂芬森迭代法速度最快，效果最好。2收斂情況9x 10分法-5 -“-10 =cc:010203040牛頓迭代法史蒂芬森迭代法0.4牛頓割線法1 卜_ I0 -1-2 Ec10102030400.20-0.2-0.4圖1.2不同迭代法下迭代值得收斂情況二分法收斂效果較差，牛頓迭代法和牛頓割線法相近，史蒂芬森迭代法收斂 次數(shù)高于1,效果最好3. 不同初值的收斂情況O543210-1-2-3-450100150分法牛頓迭代法2.521.50.5020406080圖1.

14、3二分法，牛頓迭代法下不同初值的收斂情況牛頓割線法史蒂芬森迭代法86420-2-4-6-8圖1.4牛頓割線法，史蒂芬森迭代法下不同初值的收斂情況1. 二分法的五個初始區(qū)間分別為1H 丨-I- J 2. 牛頓迭代法的五個初始值分別為 -2,-4,4/6；3. 牛頓割線法的五個初始區(qū)間分別為1 _ I4. 史蒂芬森迭代法的五個初始值分別為6-14九6；由圖可知，它們最終均達到收斂。收斂性分析及結(jié)論：1. 二分法收斂較慢且不能求解崇根，但算法簡單；此處牛頓法具有了平方收斂； 從迭代次數(shù)上看，牛頓割線法較牛頓法的多，所以收斂性較差，是超線性收 斂；史蒂芬森迭代法收斂效果最好。2. 因為牛頓迭代法是局部

15、的二次收斂，所以要注重初值的選取，本次實驗中選 擇的初值均得到了收斂，效果比較好。牛頓割線法也應(yīng)注意初值的選取。（第三章）實驗題目：1. 區(qū)間1作等距劃分：2Xi. 1 + W h k = 0丄Kn以U為結(jié)點對函數(shù)R進行插值逼近。（1） 分別取-丨，dt；用牛頓插值對 進行逼近，并在同一坐標系下做出函數(shù)的圖形，進行比較。寫出插值函數(shù)對的逼近程度與節(jié)點個數(shù)的關(guān)系，并分析原因；（2） 試用三次樣條插值對 進行逼近，在同一坐標下畫出圖形，觀察樣條插值函數(shù)對 的逼近程度與節(jié)點個數(shù)的關(guān)系；（3）整體插值有何局限性？如何避免？2. 已知一組數(shù)據(jù)如下，求其擬合曲線.表2.1數(shù)據(jù)表0234781106.42

16、108.2 109.5110109.93 110.49 110.59 110.6 110.76111111.22（1）求以上數(shù)據(jù)形如yg=% + cix+巧 的擬合曲線及其平方誤差；b（2） 求以上數(shù)據(jù)形如，的擬合曲線及其平方誤差；（3） 通過觀察（1） （2）的結(jié)果，寫出你對數(shù)據(jù)擬合的認識.問題分析：題目除上述要求之外還有以下幾點：1. 明確整體插值和分段插值的不同。牛頓插值多項式屬于整體插值，三次樣條 插值屬于低次分段插值。2. 將結(jié)果在同一坐標下繪制出。但是為了方便分析節(jié)點個數(shù)對于插值效果的影響，也可以單獨繪制。3. 第二題中為了確定各個參數(shù)的大小，可以進行適當變換，轉(zhuǎn)化為線性，運用最小

17、二乘法，得到擬合。實驗原理：牛頓插值多項式：對于給定的插值節(jié)點 心Vf構(gòu)造次數(shù)不超過：的插值多項式=歡0 +聽仗 FJ十幻& - X-巧）+氣& -咒0）（上-心）（工-1）使其滿足插值條件%（&）=兒（旺）2 015W這樣得到的插值多項式. 稱為插值多項式。系數(shù)：為差商，可以通 過構(gòu)造差商表得到。三次樣條插值：三次樣條插值函數(shù)|S（X）在每個小區(qū)間耳上為三次多 項式；在全進（上存在二階連續(xù)導(dǎo)數(shù)；其次符合插值條件。啊XI臼bl中存在內(nèi) 置的三次樣條插值函數(shù)，命令為 - 實驗內(nèi)容：第一題：1. 牛頓插值函數(shù)的構(gòu)造代碼如下：fun ctio nf=Newto n(xO,yO)%牛頓多項式插值函數(shù)s

18、yms x;SZ=size(xO,2);a(1)=y0(1)；y(:,1)=y0;for j=2:SZn x1=1;for i=1:SZ-j+1;n x2=nx1+j-1;y(i,j)=(y(i,j-1)-y(i+1,j-1)/(xO( nx1)-x0( nx2);n x1= nx1+1;endendf=y(1,1);for j=2:SZff=y(1,j);for i=1:j-1ff=ff*(x-x0(i);endf=f+ff;endend2. 牛頓和三次樣條插值情況及比較:代碼如下：clear all ;clc;syms x;fx=1/(5+xA2);%牛頓多項式插值x0=-1:0.01:1

19、;yO=subs(fx,xO);n 1=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); f3=Newton(:x3,y3); y03=subs(f3,x0); x4=-1:h1 (4):1; y4=subs(fx,x4)

20、; 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); x7=-1:h2(2):1; y7=subs(fx,x7); f7=Newton(:x7,y7); y07=subs(f7,x0); x

21、8=-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,x0,y07,x0,y08,x0,y09,x0,y010),title(%三次樣條插值spline自帶命令x0=-5:0.01:5;yO=s

22、ubs(fx,xO);figure,subplot(2,3,1),plot(x0,y0),title(精確圖);y11=spl in e(x1,y1,x0);subplot(2,3,2),plot(x0,y11),title( y12=spli ne(x2,y2,x0);subplot(2,3,3),plot(x0,y12),title( y13=spli ne(x3,y3,x0);subplot(2,3,4),plot(x0,y13), title( y14=spli ne(x4,y4,x0);subplot(2,3,5),plot(x0,y14),title( y15=spli ne(x5

23、,y5,x0);subplot(2,3,6),plot(x0,y15), title( figure,Plot(x0,y0,x0,y11,x0,y12,x0,y13,x0,y14,x0,y15) ,titl e(第二題：代碼如下：clear all ;clc;x=2 3 4 7 8 10 11 14 16 18 19;所有結(jié)果);龍格現(xiàn)象);n=1);n=5);n=10);n=20);n=25);所有結(jié)果);y=106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2;%Age nda 1A=(x.*x) x on e

24、s(11,1);A1=A*A;y1=A*y;c1=A1y1x1=2:0.1:19;y1=polyval(c1,x1);onplot(x,y, k*,x1,y1,r ),gtext( 曲線一),holdy01= polyval(c1,x);delta1_2=sum(y-y01).A2)%Age nda 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=x1;y2=subs(f,x2);plot(x2,y2, b ),gtext(

25、曲線二)；y02=subs(f,x);delta2_2=sum(y-yO2).A2)實驗結(jié)果：第一題：1. 牛頓插值結(jié)果所有結(jié)果圖2.1不同節(jié)點數(shù)的牛頓插值多項式綜合圖龍格現(xiàn)象圖2.2龍格現(xiàn)象由圖可知，2.三次樣條插值結(jié)果精確圖n=10.22 - r0.8L-0.15 1 1-j10.6V40.14/ 0.4 10.05-yx j00.2 、 / 10-10n=5-5-5-5-505n=250.80.60.40.20-5-5圖2.3不同節(jié)點下的三次樣條插值結(jié)果0.70.60.50.4所有結(jié)果0.3 -0-5-4-3-2-1圖2.4不同節(jié)點數(shù)的三次樣條的綜合圖由圖可知，牛頓插值多項式在較小的時候

26、如陽“犖總訶差值效果良好，當 變大時如就出現(xiàn)了龍格現(xiàn)象，三次樣條在各個子區(qū)間內(nèi) 為三次多項式，擬合效果好.第二題1.第一問的多項式擬合得到的擬合曲線為y(x) = 106.2927 + 0.6264X - 0-0205X2平方誤差為=2.7796第二問的擬合曲線為0.0903y(x) = llL4940e *平方誤差為如=0.47192. 擬合曲線如圖所示112111110109108107106112111110109108107106* Mr10121416 182010121416 1820圖2.5擬合情況從圖中可以看出曲線二大體符合黑點的分布情況，擬合效果較好。結(jié)論：整體插值隨著節(jié)點個

27、數(shù)的增多，多項式的次數(shù)也在升高。高次多項式的插值 會出現(xiàn)龍格現(xiàn)象，震蕩劇烈，逼近效果不理想。但是當節(jié)點很多時，低次插值的 精度又不夠，所以為了避免這一局限性采用分段低次插值。 其中三次樣條插值有 良好的收斂性和光滑性，效果較好。數(shù)據(jù)擬合時，可以選擇趨勢相似的函數(shù)形式， 在求解相關(guān)參量時，將函數(shù)進 行適當變換，從而運用最小二乘法得到擬合結(jié)果。（第四章）實驗題目：設(shè)計區(qū)間分半求積算法、龍貝格求積算法和自適應(yīng)辛普森求積算法的程序， 觀察n= 140400,500時，積分J工皿冋血的結(jié)果，并給出相應(yīng)的評價. 問題分析：1. 分半求積算法即為復(fù)合梯度求積。2. 因為越大，“在內(nèi)震蕩地越嚴重，自適應(yīng)辛普森

28、求 積算法很難適用，計算復(fù)雜度會很高。所以選取辛普森求積算法代替。實驗原理：分半求積算法：將區(qū)間 J等分為個小區(qū)間，分點為b （1= a + kK h =, k =利用定積分性質(zhì)，有fWdxn- 1”紺 + Kxk +1) +每個小區(qū)間上均用梯形求積公式，有cb 旗J /（g -譏令- 1札ji - 1即為復(fù)合梯形公式。龍貝格求積算法：將步長為的復(fù)合梯形公式進行査爾遜外推，就得到龍貝格求積公式。自適應(yīng)辛普森求積算法：按照被奇函數(shù)在區(qū)間上的變化來安排求積結(jié)點的辛 普森算法。實驗內(nèi)容:1. 寫得各個求積算法的函數(shù)。 代碼如下：(1) 分半求積算法：fun cti onf=CompositeTrap

29、ezoida( n)syms x;y=x*x*cos( n*x);m=1:100;a=-1;b=1;for i=1:100h=(b-a)/m(i);k=0:m(i);xk=a+k*h;yk=subs(y,xk);f(i)=0;for j=1:m(i)f(i)=f(i)+(yk(j)+yk(j+1);endf(i)=h/2*f(i);endend(2) 辛普森算法fun cti onf=CompositeSimpso n(n)syms x;y=x*x*cos( n*x);m=1:100;a=-1;b=1;for i=1:100h=(b-a)/(2*m(i); k=0:(2*m(i);xk=a+k

30、*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)+y k(2*j+1);endf(i)=h/3*f(i);endm;f;end(4 )龍貝格求積算法fun cti onf=Romberg (n)epsilo n=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 1h=h/2;n=len gth(pO);k=1:(2Am-1);xk=a+k*h;yk=subs(y,

31、xk);su=sum(yk);S( 1)=h/2*(fa+fb)+h*su;for i=1:mS(i+1)=S(i)+(S(i)-T(i)/(4Ai-1)；endif (abs(S(m+1)-T(m)epsilon) break ;endT=S;m=m+1;endf=S(m+1);end(5)自適應(yīng)辛普森求積算法fun cti onf=AdaptiveSimps on(n)epsilo n=0.001;syms x;y=x*x*cos( n*x);a=-1;b=1;h=b-a;p=a b;p0=p；ep=epsil on;m=0;q=0；f=0;while 1n 1=le ngth(ep);i

32、f (n=1) break ; endh=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(3); s2=h/12*(yk(3)+4*yk+yk(5);if (abs(s0-s1-s2)=2ep=ep(2:n 1);endq=q+1;elsem=m+1;p0=yk(1),yk (3),p0(2: n);if n 1=1ep=ep(1)/2 ep(1)/2;elseep=ep(1)/2 ep(1)/2ep(2:n 1);endif q=0p=p0;else

33、p=p(1:q) p0;endendendend2.各個求積公式下的積分結(jié)果代碼如下:clear all ;clc;%復(fù)合梯形法y11=CompositeTrapezoida(1);y12=CompositeTrapezoida(10);y13=CompositeTrapezoida(100);y14=CompositeTrapezoida(500);1:100;y11;y12;y13;y14%復(fù)合辛普森法y21=CompositeSimps on( 1);y22=CompositeSimpso n(10);y23=CompositeSimpso n(100);y24=CompositeSim

34、pso n(500);2:2:200;y21;y22;y23;y24實驗結(jié)果：%龍貝格法y31=Romberg(1);y32=Romberg(10);y33=Romberg(100);y34=Romberg(500);y31 y32 y33 y34%自適應(yīng)辛普森法y41=AdaptiveSimps on (1)%以下復(fù)雜度太高不岀結(jié)果%y42=AdaptiveSimpso n(10);%y43=AdaptiveSimpso n(100);%y44=AdaptiveSimpso n(500);%y41 y42 y43 y44表3.1不同積分法得到的結(jié)果n110100500復(fù)合梯形法0.478-0

35、.1399401910-0.0060.002復(fù)合辛普森法0.478-0.140-0.0090.0072247738龍貝格法0.4782672574-06112212333-0.2492135619結(jié)果評價:龍貝格法的精度高，最接近精確解。因為被積分函數(shù)為f(Q =幾皿(nx)隨著 的增大，震蕩越劇烈。所以積分面積越來越小。(第五章)實驗題目:1R1.分別取步長山=用龍格-庫塔法求解y =tyy(l) = lt計算到y(tǒng)(t)=二用并與精確解比較，觀察龍格-庫塔法的收斂性.2. 分別取步長-1 11,1 -用顯示歐拉法和隱式歐拉法求解y =- 50y，y（o）= wo ”

36、由結(jié)果分析算法的穩(wěn)定性.r isin t3. 選擇某常微分方程初值問題的數(shù)值方法計算丿。曲的近似值，并保證有 4 位有效數(shù)字.問題分析:1. 龍格-庫塔有很多的格式，第一題以常用的四級四階龍格庫塔方法為例。sin t2. 假設(shè)， ，由于gin tt第三題可以轉(zhuǎn)化為求解常微分方程的問題實驗原理：四級四階龍格-庫塔格式丹+】=兒+評1 +生+ 2心+心) 比 | = f ) hh嶺74 +初+尹)hh心= f(q + 2必+尹 知=f(ti + hlyi 4- hkj實驗內(nèi)容：1. 求解第一題代碼如下：clear all ;clc%以四級四階龍格庫塔格式為例t0=1;h=0.5 0.1 0.05

37、0.01;n=1./h+1;y1(1)=1;y2(1)=1;y3(1)=1;y4(1)=1;t1=1:0.5:2;t2=1:0.1:2;t3=1:0.0 5:2;t4=1:0.01:2;for j=2: n(1)k1=t1(j_1)*y1(j_1)A(1/3);k2=(t1(j-1)+h/2)*(y1(j-1)+(h(1)/2)*k1F(1/3);k3=(t1(j-1)+h/2)*(y1(j-1)+(h(1)/2)*k2)A(1/3);k4=t1(j)*(y1(j-1)+h(1)*k3)A(1/3);y1(j)=y1(j-1)+h(1)/6*(k1+2*k2+2*k3+k4);endfor j

38、=2: n(2)k1= t2(j-1)*y2(j-1)A(1/3);k2=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k1)A(1/3);k3=(t2(j-1)+h(2)/2)*(y2(j-1)+(h(2)/2)*k2)A(1/3);k4=t2(j)*(y2(j-1)+h(2)*k3)A(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)A(1/3);k2=(t3(j-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k1)A(1/3);k3=(t3(j

39、-1)+h(3)/2)*(y3(j-1)+(h(3)/2)*k2)A(1/3);k4=t3(j)*(y3(j-1)+h(3)*k3)A(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)A(1/3);k2=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k1)A(1/3);k3=(t4(j-1)+h(4)/2)*(y4(j-1)+(h(4)/2)*k2)A(1/3);k4=t4(j)*(y4(j-1)+h(4)*k3)A(1/3);y4(j)=y4(j-1)+h(4

40、)/6*(k1+2*k2+2*k3+k4);endx=1:0.01:2;syms tyy=(tA2+2)/3)A(3/2);yy1=subs(yy,x);figure,subplot(2,3,1),plot(x,yy1) ,titl e( subplot(2,3,2),plot(x,yy1,t1,y1) ,titl e( subplot(2,3,3),plot(x,yy1,t2,y2),title( subplot(2,3,4),plot(x,yy1,t3,y3) ,titl e( subplot(2,3,5),plot(x,yy1,t4,y4) ,titl e(精確值)；h=0.5);h=0

41、.1;h=0.05);h=0.01);Subplot(2,3,6),plot(x,yy1,t1,y1,t2,y2,t3,y3,t4,y4),title(2.求解第二題代碼如下：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.0 5:a;t3=0:0.01:a;t4=0:0.005:a;% 顯式歐拉法y1(1)=100;y2(1)=100;y3(1)=100;y4(1)=100;所有

42、結(jié)果;forj=2:n(1)y1(j)=y1(j-1)+h(1)*(-50)*y1(j-1);endforj=2: n(2)y2(j)=y2(j-1)+h(2)*(-50)*y2(j-1);endforj=2: n(3)y3(j)=y3(j-1)+h(3)*(-50)*y3(j-1);endforj=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.

43、005);forj=2: n(1)y5(j)=y5(j-1)/(50*h(1)+1);endforj=2: n(2)y6(j)=y6(j-1)/(50*h (2)+1);endforj=2: n(3)y7(j)=y7(j-1)/(50*h (3)+1);endforj=2: n(4)y8(j)=y8(j-1)/(50*h (4)+1);endy5(1)=100;y6(1)=100;y7(1)=100;y8(1)=100;figure,plot(x,yy1,t1,y5,t2,y6,t3,y7,t4,y8).),gtext(h=0.01),gtext(h=0.005);,t1,y1,b,t1,y5

44、,r),title(h=0.1,t2,y2,b,t2,y6,r),title(h=0.05,t3,y3,b,t3,y7,r),title(h=0.01,t4,y4,b,t4,y8,r),title(h=0.005gtext( h=0.1 ),gtext( h=0.05%相同步長下的穩(wěn)定性比較figuresubplot(2,2,1),plot(x,yy1,ksubplot(2,2,2),plot(x,yy1,ksubplot(2,2,3),plot(x,yy1,ksubplot(2,2,4),plot(x,yy1,k3. 求解第三題代碼如下：%避免t=0作分母，所以用隱式歐拉法clear all

45、 ;t0=0;h=0.0005 0.0001 0.00005;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;endendendfor j=2:n(1) y1(j)=y1(j-1)+h(1)*sin(t1(j)/t1(j);for j=2:n (2) y2(j)=y2(j-1)+h (2) *si n(t2(j)/t2(j);for j=2:n (3)y3(j)=y3(j-1)+h (3) *si n(t3(j)/t3(j);format long y(1)=y1( n( 1

46、)-y1(1);y(2)=y2( n( 2)-y2(1);y(3)=y3( n(3)-y3(1);實驗結(jié)果:1.第一題結(jié)果32.521.51h=0.511.52h=0.111.52h=0.01所以結(jié)果h=0.0511.52精確值2.521.5111.52圖4.1不同步長下的四級龍格-庫塔的解2. 顯示歐拉法求解結(jié)果20001500h=0.11000500-500上=0 01 h=0.005h=0.050.050.10.150.20.25圖4.2不同步長下的顯式歐拉法的解0.25隱式歐拉法求解結(jié)果00一*T90-J80-70h=0.160-50h=0.05-h=0.0140*130-h=0.00

47、5X -20-10- -B-00.050.10.150.2圖4.3不同步長下的隱式歐拉法的解顯式與隱式的穩(wěn)定性比較結(jié)果2000h=0.110000-10000.10.20.30.4h=0.011000500.10.20.30.4h=0.055000-500-10000.4h=0.005100500.10.20.30.40.10.20.3圖4.4不同步長下顯式與隱式的穩(wěn)定性比較3. 第三題結(jié)果表4.1不同步長下的積分結(jié)果步長h0.00050.00010.00005積分結(jié)果0.94680.94600.9463題目要求保證有4位有效數(shù)字，所以所求結(jié)果為：0.9461 結(jié)果分析：四級四階龍格庫塔法具有

48、大范圍的收斂性。 四級四階龍格庫塔法每計算一步 需要計算四次的值，計算有一定的復(fù)雜性。當步長:為0.1和0.05時顯式歐拉法不具有穩(wěn)定性。隱式歐拉法在各個步長 下的穩(wěn)定情況相似，均有較好的穩(wěn)定性。（第六章）實驗題目:使用列主元高斯消去法解希爾伯特矩陣叭7方程組Hx =給定的- 及相應(yīng)的h時結(jié)果的變化，分析其中的原因并給出結(jié)論，其中問題分析：1. 希爾伯特矩陣是病態(tài)的，的微小變化會造成解發(fā)生很大變化。2. 刻畫矩陣病態(tài)程度的量為條件數(shù)。條件數(shù)越大，矩陣性能越差?？梢杂嬎銞l 件數(shù)來分析解變化很大的現(xiàn)象。實驗原理:高斯列主元素消元法：將矩陣按列取絕對值最大項進行行行調(diào)換，將其以 F部分消為即，不斷循

49、環(huán)直到H化為上三角的過程。實驗內(nèi)容：1.高斯列主元素消元法函數(shù)代碼如下：fun cti onA b=GaussElimi natio n(A,b)n=size(A,1);for i=1: n-1q=max(abs(A(i: n,i);for j=i:nif (abs(A(j,i)=q)k=j;endendu=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: nA(l,i)=A(l,i)/A(i,i);b(l)=b(l)-A(l,i)*b(i);for s=i+1: nA(l,s)=A(l,s)-A(l,i)*A(i,s);endendend end2. 計算不同階數(shù)下的希爾伯特矩陣的線性方程的