數(shù)值計(jì)算B大作業(yè)

1、課 程 設(shè) 計(jì)課程名稱： 數(shù)值計(jì)算B 設(shè)計(jì)題目： 數(shù)值計(jì)算B大作業(yè) 學(xué) 號(hào)： 姓 名： 完成時(shí)間： 題目一：多項(xiàng)式插值某氣象觀測(cè)站在8：00（AM）開始每隔10分鐘對(duì)天氣作如下觀測(cè)，用三次多項(xiàng)式插值函數(shù)（Newton）逼近如下曲線，插值節(jié)點(diǎn)數(shù)據(jù)如上表，并求出9點(diǎn)30分該地區(qū)的溫度（x=10）。x12345678y22.523.324.421.7025.228.524.825.4二、數(shù)學(xué)原理假設(shè)有n+1個(gè)不同的節(jié)點(diǎn)及函數(shù)在節(jié)點(diǎn)上的值（x，y），（x，y），插值多項(xiàng)式有如下形式： （1）其中系數(shù)（i=0,1,2n）為特定系數(shù)，可由插值樣條（i=0,1,2n）確定。根據(jù)均差的定義，把x看成a,b上的

2、一點(diǎn),可得f(x)= f（）+f（）fx, = f+fx, （）fx, ，x= fx, ，x+ fx, ，x(x-x)綜合以上式子，把后一式代入前一式，可得到： f(x)= f+f（）+ f（）（）+ fx, ，x（）(x-x)+ fx, ，x,= N（x）+其中N（x）= f+f（）+ f（）（）+ fx, ，x（）(x-x) （2）= f(x)- N（x）= fx, ，x, （3） =（）(x-x)Newton插值的系數(shù)（i=0,1,2n）可以用差商表示。一般有 （k=0,1,2，n ） （4）把（4）代入（1）得到滿足插值條件N（i=0,1,2，n）的n次Newton插值多項(xiàng)式N（x）=

3、f（）+f（）+f（）（）+f（）（）（）.其中插值余項(xiàng)為： 介于之間。三、程序設(shè)計(jì)function y,A,C,L=newdscg(X,Y,x,M) % y為對(duì)應(yīng)x的值，A為差商表，C為多項(xiàng)式系數(shù)，L為多項(xiàng)式% X為給定節(jié)點(diǎn)，Y為節(jié)點(diǎn)值，x為待求節(jié)點(diǎn)n=length(X); m=length(x); % n為X的長(zhǎng)度f(wàn)or 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); e

4、nd 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,poly(X(k); d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z); %輸出y值endL(k,:)=poly2sym(C); %輸出多項(xiàng)式>> syms M,X=1,3,5,7;Y=22.5,24.4,25.2,24.8;x=10;>> y,A,C,L=newdscg(X,Y,x,M)y = 21.7313A = 22.5000 0 0

5、0 24.4000 0.9500 0 0 25.2000 0.4000 -0.1375 0 24.8000 -0.2000 -0.1500 -0.0021C = -0.0021 -0.1187 1.4521 21.1688 L = - x3/480 - (19*x2)/160 + (697*x)/480 + 3387/1604、 結(jié)果分析和討論 對(duì)于不超過三次的插值多項(xiàng)式，x如果選取1,3,5,7這三個(gè)點(diǎn)能夠得到較好的三次插值多項(xiàng)式L=-0.0021x3-0.1187x2+1.4521x+21.1688。當(dāng)x=10時(shí)，也即9點(diǎn)30分時(shí)的溫度為21.7317度，結(jié)果分析知此值應(yīng)是偏小的。對(duì)于選取

6、不同的插值節(jié)點(diǎn)，能夠得到不同的插值多項(xiàng)式，誤差也不盡相同。5、 完成題目的體會(huì)與收獲 牛頓插值法的重要一點(diǎn)就是對(duì)插值節(jié)點(diǎn)的選取，通過本題的編程很好的加深了對(duì)其概念的理解以及提高了應(yīng)用牛頓插值法的能力，學(xué)會(huì)了運(yùn)用Matlab軟件對(duì)牛頓插值法相關(guān)問題進(jìn)行編程求解，對(duì)Matlab計(jì)算方法與程序編輯更加熟悉。使我對(duì)這類問題的理解有了很大的提升。題目二：曲線擬合在某鋼鐵廠冶煉過程中，根據(jù)統(tǒng)計(jì)數(shù)據(jù)的含碳量與時(shí)間關(guān)系，試用最小二乘法擬合含碳量與時(shí)間t的擬合曲線，并繪制曲線擬合圖形。t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87

7、 4.15 4.37 4.51 4.58 4.02 4.64 2、 數(shù)學(xué)原理從整體上考慮近似函數(shù)同所給數(shù)據(jù)點(diǎn)(i=0,1,m)誤差(i=0,1,m)的大小，常用的方法有以下三種：一是誤差(i=0,1,m)絕對(duì)值的最大值，即誤差 向量的范數(shù)；二是誤差絕對(duì)值的和，即誤差向量r的1范數(shù)；三是誤差平方和的算術(shù)平方根，即誤差向量r的2范數(shù)；前兩種方法簡(jiǎn)單、自然，但不便于微分運(yùn)算 ，后一種方法相當(dāng)于考慮 2范數(shù)的平方，因此在曲線擬合中常采用誤差平方和來 度量誤差(i=0，1，m)的整體大小。數(shù)據(jù)擬合的具體作法是：對(duì)給定數(shù)據(jù) (i=0,1,，m)，在取定的函數(shù)類中,求,使誤差(i=0,1,m)的平方和最小，

8、即 =從幾何意義上講，就是尋求與給定點(diǎn)(i=0,1,m)的距離平方和為最小的曲線。函數(shù)稱為擬合 函數(shù)或最小二乘解，求擬合函數(shù)的方法稱為曲線擬合的最小二乘法。在曲線擬合中，函數(shù)類可有不同的選取方法。三、程序設(shè)計(jì)>> x=0,5,10,15,20,25,30,35,40,45,50,55;>> y=0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64*10(-4);>> p=polyfit(x,y,2)p = 1.0e-04 * -0.0024 0.2037 0.2305>> plot(x,

9、y,'r')4、 結(jié)果分析和討論 通過最小二乘法得到的曲線擬合多項(xiàng)式是p=(-0.0024x2+0.2037x+0.2305)*10-4。利用Matlab軟件調(diào)用最小二乘法函數(shù)即可得到多項(xiàng)式擬合函數(shù)，由于數(shù)據(jù)較少得到的擬合曲線不夠光滑，可能會(huì)存在一定的誤差。擬合曲線總體呈現(xiàn)增函數(shù)趨勢(shì)，與數(shù)據(jù)較為吻合。5、 完成題目的體會(huì)與收獲 曲線擬合較常用的方法之一就包括最小二乘法，因此能夠熟練使用Matlab進(jìn)行數(shù)據(jù)的曲線擬合變得尤為重要。通過完成本題的擬合過程，對(duì)于最小二乘法曲線擬合的操作更加的熟練，能夠較好的完成各類數(shù)據(jù)的擬合。題目三：線性方程組求解分別利用LU分解；平方根法與改進(jìn)平方

10、根法；追趕法求解下列幾個(gè)不同類型的線性方程組，并與準(zhǔn)確值比較結(jié)果，分析誤差產(chǎn)生原因。（1）設(shè)線性方程組 二、數(shù)學(xué)原理 將A分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣U，即：A=LU,其中 L=， U= 解Ax=b 的問題就等價(jià)于要求解兩個(gè)三角形方程組： Ly=b,求y;   Ux=y,求x。 設(shè)A為非奇異矩陣,且有分解式A=LU, L為單位下三角陣,U為上三角陣。 L,U的元素可以有n步直接計(jì)算定出。用直接三角分解法解Ax=b（要求A的所有順序主子式都不為零）的計(jì)算公式： , ,i=2,3,，n. 計(jì)算U的第r行，L的第r列元素(i=2,3,n)：   &#

11、160; ， i=r,r+1,n; , i=r+1,n,且rn. 求解Ly=b，Ux=y的計(jì)算公式； 三、程序設(shè)計(jì)function f1=LUresolve(a,b);n,n=size(a); % 行數(shù)z=size(b) % b的行數(shù)l=; % l矩陣u=; % u矩陣for j=1:1:n u(1,j)=a(1,j);endfor i=2:1:n l(i,1)=a(i,1)/u(1,1);endfor i=2:1:(n-1) sum1=0; for k=1:1:(i-1) sum1=l(i,k)*u(k,i)+sum1; end u(i,i)=a(i,i)-sum1; sum2=0; sum

12、3=0; for j=(i+1):1:n for k=1:1:(i-1) sum2=sum2+l(i,k)*u(k,j); sum3=sum3+l(j,k)*u(k,i); end u(i,j)=a(i,j)-sum2; l(j,i)=(a(j,i)-sum3)/u(i,i); endendfor i=1:1:(n-1) l(i,i)=1; l(i,n)=0;endl(n,n)=1;sum4=0;for k=1:1:(n-1) sum4=sum4+l(n,k)*u(k,n);endu(n,n)=a(n,n)-sum4;l %輸出l矩陣u %輸出u矩陣y=lb %輸出向量yx=uy %輸出向量x

13、>> a=4 2 -3 -1 2 1 0 0 0 0 8 6 -5 -3 6 5 0 1 0 0 4 2 -2 -1 3 2 -1 0 3 1 0 -2 1 5 -1 3 -1 1 9 4 -4 2 6 -1 6 7 -3 3 2 3 8 6 -8 5 7 17 2 6 -3 5 0 2 -1 3 -4 2 5 3 0 1 16 10 -11 -9 17 34 2 -1 2 2 4 6 2 -7 13 9 2 0 12 4 0 0 -1 8 -3 -24 -8 6 3 -1;>> b=5;12;3;2;3;46;13;38;19;-21;>> LUresol

14、ve(a,b)z = 10 1l = 1.0000 0 0 0 0 0 0 0 0 0 2.0000 1.0000 0 0 0 0 0 0 0 0 1.0000 0 1.0000 0 0 0 0 0 0 0 0 -2.0000 3.0000 1.0000 0 0 0 0 0 0 -1.0000 1.0000 4.0000 -0.4667 1.0000 0 0 0 0 0 2.0000 1.0000 -5.0000 -0.2667 -1.6441 1.0000 0 0 0 0 0 -1.0000 3.0000 -0.0667 -0.0847 0.2969 1.0000 0 0 0 4.0000

15、-1.0000 6.0000 -0.2667 -3.7627 2.7303 3.1863 1.0000 0 0 1.0000 -4.0000 26.0000 1.2667 5.0000 -0.0182 -18.5356 10.7592 1.0000 0 0 -7.0000 30.0000 4.6000 21.7288 -11.1148 -59.9306 29.2179 0.9856 1.0000u = 1.0e+03 * 0.0040 0.0020 -0.0030 -0.0010 0.0020 0.0010 0 0 0 0 0 0.0020 0.0010 0.0050 0.0100 0.007

16、0 0.0020 0.0030 0.0020 0.0020 0 0 0.0010 0 0.0020 0 -0.0030 -0.0020 0.0010 -0.0010 0 0 0 0.0150 0.0130 0.0310 0.0400 0.0540 0.0630 0.0650 0 0 0 0 -0.0039 0.0155 0.0341 0.0703 0.0927 0.1261 0 0 0 0 0 0.0417 0.0518 0.1728 0.3361 0.5617 0 0 0 0 0 0 0.0062 -0.0297 -0.1214 -0.2672 0 0 0 0 0 0 0 -0.0840 -

17、0.1669 -0.3495 0 0 0 0 0 0 0 0 -0.9987 -1.8562 0 0 0 0 0 0 0 0 0 -0.7229y = 1.0e+03 * 0.0050 0.0020 -0.0020 0.0120 0.0196 0.0594 0.0058 -0.0718 0.8428 1.8498x = 8.3778 41.2714 10.5278 30.7380 -28.0438 7.3550 -14.8478 3.7273 3.9121 -2.55894、 結(jié)果分析和討論 LU分解所得到的結(jié)果x=(8.3778,41.2714,10.5278,30.7380,-28.043

18、8,7.3550,-14.8478,3.7273,3.9121,-2.5589)T與精確解具有很大的誤差。這是由于系數(shù)矩陣本身的數(shù)值性質(zhì)所決定的，由于計(jì)算過程中并未進(jìn)行選主元的過程，所以由系數(shù)矩陣產(chǎn)生的L和U矩陣就具有了很大的數(shù)值問題。由L和U所求出的x和y就會(huì)產(chǎn)生很大的誤差。這是由矩陣本身的數(shù)值問題所引起的，與算法本身無關(guān)，消除誤差的關(guān)鍵在于計(jì)算過程中需要進(jìn)行選主元。5、 完成題目的體會(huì)與收獲對(duì)于其他解線性方程組的方法來看，LU分解是較為高效的，理解并熟練運(yùn)用LU分解對(duì)于學(xué)習(xí)數(shù)值計(jì)算與解線性方程組都有很大的幫助。在進(jìn)行分解的過程中應(yīng)注意矩陣的數(shù)值問題與如何選取主元的問題，這樣才能得到方程組的

19、精確解，否則將產(chǎn)生非常大的誤差。在進(jìn)行分解時(shí)應(yīng)該格外注意，因?yàn)橄禂?shù)矩陣存在很多的數(shù)值問題。（2）設(shè)對(duì)稱正定陣系數(shù)陣線方程組 2、 數(shù)學(xué)原理1、 平方根法解n階線性方程組Ax=b的choleskly方法也叫做平方根法，這里對(duì)系數(shù)矩陣A是有要求的，需要A是對(duì)稱正定矩陣，根據(jù)數(shù)值分析的相關(guān)理論，如果A對(duì)稱正定，那么系數(shù)矩陣就可以被分解為的形式，其中L是下三角矩陣，將其代入Ax=b中，可得：進(jìn)行如下分解：那么就可先計(jì)算y,再計(jì)算x，由于L是下三角矩陣，是上三角矩陣，這樣的計(jì)算比直接使用A計(jì)算簡(jiǎn)便，同時(shí)你應(yīng)該也發(fā)現(xiàn)了工作量就轉(zhuǎn)移到了矩陣的分解上面，那么對(duì)于對(duì)稱正定矩陣A進(jìn)行Cholesky分解，我再描述

20、一下過程吧：如果你對(duì)原理很清楚那么這一段可以直接跳過的。設(shè)，即其中第1步，由矩陣乘法，故求得一般的，設(shè)矩陣L的前k-1列元素已經(jīng)求出第k步，由矩陣乘法得于是2、 改進(jìn)平方根法在平方根的基礎(chǔ)上，為了避免開方運(yùn)算,所以用計(jì)算；其中，；得 按行計(jì)算的元素及對(duì)元素公式 對(duì)于. 計(jì)算出的第行元素后，存放在的第行相置，然后再計(jì)算的第行元素，存放在的第行.的對(duì)角元素存放在的相應(yīng)位置. 對(duì)稱正定矩陣按分解和按分解計(jì)算量差不多，但分解不需要開放計(jì)算。求解, 的計(jì)算公式分別如下公式。 3、 程序設(shè)計(jì)1、平方根法function x=pfpf(A,b)%楚列斯基分解 求解正定矩陣的線性代數(shù)方程 A=LL 先求LY=

21、b 再用LX=Y 即可以求出解Xn,n=size(A);L(1,1)=sqrt(A(1,1);for k=2:n L(k,1)=A(k,1)/L(1,1);endfor k=2:n-1 L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).2); for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)/L(k,k); endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).2);%解下三角方程組Ly=b 相應(yīng)的遞推公式如下，求出y矩陣y=zeros(n,1);%先生成方程組的因變量的位置，給定y的初始值f

22、or k=1:n j=1:k-1; y(k)=(b(k)-L(k,j)*y(j)/L(k,k);end%解上三角方程組 LX=Y 遞推公式如下，可求出X矩陣x=zeros(n,1);U=L'%求上對(duì)角矩陣for k=n:-1:1 j=k+1:n; x(k)=(y(k)-U(k,j)*x(j)/U(k,k);end >> A=4,2,-4,0,2,4,0,0 2,2,-1,-2,1,3,2,0 -4,-1,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,3 2,1,-8,-1,22,4,-10,-3 4,3,-3,-4,4,11,1,-4 0,2,5,-3,

23、-10,1,14,2 0,0,6,3,-3,-4,2,19;>> b=0;-6;20;23;9;-22;-15;45;>> x=pfpf(A,b)x = 121.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.01852、改進(jìn)平方根法function x=improvecholesky(A,b,n) %用改進(jìn)平方根法求解Ax=bL=zeros(n,n); %L為n*n矩陣D=diag(n,0); %D為n*n的主對(duì)角矩陣S=L*D;for i=1:n %L的主對(duì)角元素均為1 L(i,i)=1;endf

24、or i=1:n for j=1:n %驗(yàn)證A是否為對(duì)稱正定矩陣if (eig(A)<=0)|(A(i,j)=A(j,i) %A的特征值小于0或A非對(duì)稱時(shí)，輸出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %將A分解使得A=LDLTfor i=2:n for j=1:i-1 S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)'); L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1); end D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:

25、i-1)');endy=zeros(n,1); % x，y為n*1階矩陣x=zeros(n,1);for i=1:n y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)/D(i,i); %通過 LDy=b解得y的值endfor i=n:-1:1 x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n); %通過LTx=y解得x的值end>> A=4,2,-4,0,2,4,0,0 2,2,-1,-2,1,3,2,0 -4,-1,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,3 2,1,-8

26、,-1,22,4,-10,-3 4,3,-3,-4,4,11,1,-4 0,2,5,-3,-10,1,14,2 0,0,6,3,-3,-4,2,19;>> b=0;-6;20;23;9;-22;-15;45;>> n=8;>> x=improvecholesky(A,b,n)x = 121.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.01854、 結(jié)果分析和討論 平方根法和改進(jìn)平方根法求解線性方程組的解為x=（121.1481，-140.1127，29.7515，-60.1528，10

27、.9120，-26.7963，5.4259，-2.0185）T。與精確解相比較也存在很大的誤差，雖然系數(shù)矩陣的對(duì)角元素都大于零，原則上可以不必選擇主元，但由于矩陣的數(shù)值問題較大，不選主元的結(jié)果就是產(chǎn)生很大的誤差，所以在求解的過程中還是應(yīng)該選擇主元以此消除誤差，提高精度。 5、 完成題目的體會(huì)與收獲 對(duì)稱正定矩陣的平方根法及改進(jìn)平方根法是目前解決這類問題的最有效的方法之一，合理利用的話，能夠產(chǎn)生很好的求解效果。改進(jìn)平方根法較平方根法，因?yàn)椴挥眠M(jìn)行開方運(yùn)算，所以具有一定的求解優(yōu)勢(shì)。通過求解此題，使我學(xué)會(huì)了如何運(yùn)用Matlab編程來解決平方根法和改進(jìn)平方根法問題。（3）三對(duì)角形線性方程組 2、 數(shù)學(xué)

28、原理 設(shè)系數(shù)矩陣為三對(duì)角矩陣則方程組Ax=f稱為三對(duì)角方程組。 設(shè)矩陣A非奇異，A有Crout分解A=LU，其中L為下三角矩陣，U為單位上三角矩陣，記可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程組Ly=f得出y，再求解上三角方程組Ux=y。 事實(shí)上，求解三對(duì)角方程組的2追趕法將矩陣三角分解的計(jì)算與求解兩個(gè)三角方程組的計(jì)算放在一起，使算法更為緊湊。其計(jì)算公式為：(*)三、程序設(shè)計(jì)function x=chase(a,b,c,f)%求解線性方程組Ax=f,其中A是三對(duì)角陣%a是矩陣A的下對(duì)角線元素a(1)=0%b是矩陣A的對(duì)角線元素%c是矩陣A的上對(duì)角線元素c(n)=0%f是方程組的

29、右端向量n=length(f);x=zeros(1,n);y=zeros(1,n);d=zeros(1,n);u= zeros(1,n);%預(yù)處理d(1)=b(1);for i=1:n-1u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i);end%追的過程y(1)=f(1)/d(1);for i=2:n y(i)=(f(i)-a(i)*y(i-1)/d(i);end%趕的過程x(n)=y(n);for i=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);end>> a=0,-1,-1,-1,-1,-1,-1,-1,-1,-1;>>

30、; b=4,4,4,4,4,4,4,4,4,4;>> c=-1,-1,-1,-1,-1,-1,-1,-1,-1,0;>> f=7,5,-13,2,6,-12,14,-4,5,-5;>> x=chase(a,b,c,f)x =2.00001.0000 -3.00000.00001.0000 -2.00003.0000 -0.00001.0000 -1.00004、 結(jié)果分析和討論 追趕法求解的結(jié)果為x=（2,1，-3,0,1，-2,3,0,1，-1）T。求解結(jié)果與精確解一樣，這表明追趕法對(duì)于求解三對(duì)角方程組具有非常高的精度，誤差非常小。算法次數(shù)也較少，不選主元

31、也可以有效的算出精確結(jié)果，是一種計(jì)算量少而數(shù)值穩(wěn)定的方法。5、 完成題目的體會(huì)與收獲 通過本題的求解，加深了對(duì)追趕法求解三對(duì)角方程組的算法原理的理解。運(yùn)用追趕法的Matlab編程，并學(xué)會(huì)了又一種求解特殊方程組的方法。在計(jì)算量方面，追趕法有著巨大的優(yōu)勢(shì)，因此在可能的情況下應(yīng)優(yōu)先使用追趕法。加深了對(duì)數(shù)值計(jì)算教材知識(shí)的理解，收獲非常大。題目四：微分方程數(shù)值解在傳染病的傳染理論中，一個(gè)比較初等的微分方程可用于預(yù)測(cè)在任何時(shí)刻人群中受傳染的數(shù)量，只要做了適當(dāng)?shù)暮?jiǎn)化。特別的，假定在一個(gè)固定的人群中，所有的人具有同樣的機(jī)會(huì)被感染，且一感染就保持這種狀態(tài)。假設(shè)表示在時(shí)刻易受感染的人的數(shù)量，表示感染別人的人數(shù)。有

32、理由假設(shè)感染別人的人數(shù)變化的速率與和的乘積成正比，因?yàn)樗俾始热Q于感染別人的人數(shù)也取決于那個(gè)時(shí)刻易感染的人數(shù)。如果人口多的足以假定和是連續(xù)的的變量，則問題表示為，其中是常數(shù)，即為人口總數(shù)。這個(gè)方程就變?yōu)閮H包含：?jiǎn)栴}：一個(gè)地區(qū)假定，又假定時(shí)間用天來衡量，求30天結(jié)束時(shí)感染別人的人口數(shù)量的近似值2、 數(shù)學(xué)原理Eular方法： 一階線性微分方程初值問題 （1）方程離散化：差分和差商 （2） 通過初始值，依據(jù)遞推公式（2）逐步算出就為顯性的Eular方法。隱形Eular方法： （3）公式（3）即為隱式Eular公式。3、 程序設(shè)計(jì)function E2=Euler_2(fun,x0,y0,xN,N)%

33、 向后Euler公式，其中，%fun為一階微分方程的函數(shù)%x0,y0為初始條件%xN為取值范圍的一個(gè)端點(diǎn)%h為區(qū)間步長(zhǎng)%N為區(qū)間的個(gè)數(shù)%x為xN構(gòu)成的向量%y為yN構(gòu)成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N%用迭代法求y(n+1) x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n); for k=1:3 z1=y(n)+h*feval(fun,x(n+1),z0); if abs(z1-z0)<1e-3 break; end z0=z1; end y

34、(n+1)=z1;end T=x',y'function z=f(x,y)z=(2*10(-6)*(100000-y)*y;>> Euler_2('f',0,1000,29,29)T = 1.0e+04 * 0 0.1000 0.0001 0.1246 0.0002 0.1551 0.0003 0.1929 0.0004 0.2396 0.0005 0.2972 0.0006 0.3680 0.0007 0.4548 0.0008 0.5605 0.0009 0.6886 0.0010 0.8429 0.0011 1.0271 0.0012 1.24

35、50 0.0013 1.4999 0.0014 1.7943 0.0015 2.1295 0.0016 2.5049 0.0017 2.9182 0.0018 3.3647 0.0019 3.8377 0.0020 4.3287 0.0021 4.8281 0.0022 5.3260 0.0023 5.8128 0.0024 6.2800 0.0025 6.7208 0.0026 7.1300 0.0027 7.5045 0.0028 7.8428 0.0029 8.14494、 結(jié)果分析和討論 用隱式歐拉法求得第30天時(shí)感染別人的人口數(shù)量的近似值為y（29）=81449人。隱式歐拉法較之歐拉

36、法的誤差更小，求取的結(jié)果具有一定的可參考價(jià)值。對(duì)于結(jié)果存在的誤差，只能通過改變算法實(shí)現(xiàn)。5、 完成題目的體會(huì)與收獲 求解初值問題，歐拉法是最簡(jiǎn)單的數(shù)值解法。而隱式歐拉法更進(jìn)一步減少了歐拉法的誤差。隱式歐拉法的原理依舊十分的簡(jiǎn)單。在對(duì)數(shù)值結(jié)果的精確度并非要求很高時(shí)可采用隱式歐拉法，這樣既能快速地得到結(jié)果，又能得到合適的結(jié)果。本題的解決對(duì)我數(shù)值計(jì)算的學(xué)習(xí)有著很大的幫助，使我對(duì)初值問題的理解更加深刻，并熟悉了初值問題的求解方法和求解過程，。題目五：數(shù)值積分在光學(xué)物理中研究矩角位的光繞射經(jīng)常會(huì)用到Fresnel積分：對(duì)于構(gòu)造一個(gè)精確到的值的表格供研究者查詢（Romberg積分算法）。2、 數(shù)學(xué)原理 龍

37、貝格算法是由遞推算法得來的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龍貝格公式。在變步長(zhǎng)的過程中探討梯形法的計(jì)算規(guī)律。設(shè)將求積區(qū)間a，b分為n個(gè)等分，則一共有n+1個(gè)等分點(diǎn)，n。這里用表示復(fù)化梯形法求得的積分值，其下標(biāo)n表示等分?jǐn)?shù)。先考察下一個(gè)字段，其中點(diǎn)，在該子段上二分前后兩個(gè)積分值顯然有下列關(guān)系將這一關(guān)系式關(guān)于k從0到n-1累加求和，即可導(dǎo)出下列遞推公式需要強(qiáng)調(diào)指出的是，上式中的代表二分前的步長(zhǎng)，而梯形法的算法簡(jiǎn)單，但精度低，收斂速度緩慢，如何提高收斂速度以節(jié)省計(jì)算量，自然式人們極為關(guān)心的。根據(jù)梯形法的誤差公式，積分值的截?cái)嗾`差大致與成正比，因此步長(zhǎng)減半后誤差將減至四分之一，既有將上

38、式移項(xiàng)整理，知由此可見，只要二分前后兩個(gè)積分值和相當(dāng)接近，就可以保證計(jì)算保證結(jié)果計(jì)算結(jié)果的誤差很小，這種直接用計(jì)算結(jié)果來估計(jì)誤差的方法稱作誤差的事后估計(jì)法。按上式，積分值的誤差大致等于，如果用這個(gè)誤差值作為的一種補(bǔ)償，可以期望，所得的應(yīng)當(dāng)是更好的結(jié)果。按上式，組合得到的近似值直接驗(yàn)證，用梯形二分前后的兩個(gè)積分值和按式組合，結(jié)果得到辛普生法的積分值。再考察辛普生法。其截?cái)嗾`差與成正比。因此，若將步長(zhǎng)折半，則誤差相應(yīng)的減至十六分之一。既有由此得不難驗(yàn)證，上式右端的值其實(shí)就等于，就是說，用辛普生法二分前后的兩個(gè)積分值和，在按上式再做線性組合，結(jié)果得到柯特斯法的積分值，既有重復(fù)同樣的手續(xù)，依據(jù)斯科特法

39、的誤差公式可進(jìn)一步導(dǎo)出龍貝格公式應(yīng)當(dāng)注意龍貝格公式已經(jīng)不屬于牛頓柯特斯公式的范疇。在步長(zhǎng)二分的過程中運(yùn)用公式加工三次，就能將粗糙的積分值逐步加工成精度較高的龍貝格，或者說，將收斂緩慢的梯形值序列加工成熟練迅速的龍貝格值序列，這種加速方法稱龍貝格算法。3、 程序設(shè)計(jì)function lbg(fx,a,t,k,e)% fx為要求的積分函數(shù)% a為要求的積分的下限% t為要求的積分的上限% k為要求的積分的最大次數(shù)% e為要求積分的結(jié)束精度T=zeros(k,k); % T為龍貝格算法遞推表T(1,1)=(t-a)*(1+fx(t)/2;for i=1:k m=0; for j=1:2(i-1) m

40、=m+fx(a+(2*j-1)*(t-a)/(2i); end T(i+1,1)=0.5*T(i,1)+(t-a)*m/2i; for n=1:i T(i+1,n+1)=(4n*T(i+1,n)-T(i,n)/(4n-1); end if abs(T(i+1,i+1)-T(i,i)<=e & i>=4 break; else ; endendfor i=1:k if T(i,1)=0 j=i; break; else ; endendif j=k error('所求次數(shù)不夠或不可積')else ;endT=T(1:j-1,1:j-1)disp('所求

41、的積分值為：') % 輸出求得的積分值disp(T(j-1,1)function fx = f(x)fx=cos(pi*x2/2);function fx = f(x)fx=sin(pi*x2/2);>> a=0;t=0.1;k=100;e=10(-4);fx=(x)cos(pi*x2/2);>> lbg(fx,a,t,k,e)T = 0.1000 0 0 0 0 0.1000 0.1000 0 0 0 0.1000 0.1000 0.1000 0 0 0.1000 0.1000 0.1000 0.1000 0 0.1000 0.1000 0.1000 0.10

42、00 0.1000所求的積分值為： 0.1000>> a=0;t=0.2;k=100;e=10(-4);fx=(x)cos(pi*x2/2);>> lbg(fx,a,t,k,e)T = 0.1998 0 0 0 0 0.1999 0.1999 0 0 0 0.1999 0.1999 0.1999 0 0 0.1999 0.1999 0.1999 0.1999 0 0.1999 0.1999 0.1999 0.1999 0.1999所求的積分值為： 0.1999>> a=0;t=0.3;k=100;e=10(-4);fx=(x)cos(pi*x2/2);>

43、;> lbg(fx,a,t,k,e)T = 0.2985 0 0 0 0 0.2992 0.2994 0 0 0 0.2993 0.2994 0.2994 0 0 0.2994 0.2994 0.2994 0.2994 0 0.2994 0.2994 0.2994 0.2994 0.2994所求的積分值為： 0.2994>> a=0;t=0.4;k=100;e=10(-4);fx=(x)cos(pi*x2/2);>> lbg(fx,a,t,k,e)T = 0.3937 0 0 0 0 0.3965 0.3974 0 0 0 0.3972 0.3975 0.3975 0 0 0.3974 0.3975 0.3975 0.3975 0 0.3975 0.397