計(jì)算方法報(bào)告程序

1、計(jì)算方法作業(yè)目錄1 作業(yè)內(nèi)容32 解題思路和算法組織4解題思路41. Gauss消去法42. LDL和Cholesky分解43. 追趕法54. 迭代法求解線性方程組55. 牛頓插值和三樣條插值的比較56. Simpson公式用于求解積分方程57. Romberg方法計(jì)算定積分58. 二分法和割線法6程序清單(雙欄)61. Gauss消去法62. LDL分解73. cholesky分解74. Jacobi迭代75. Jacobi-sediel迭代函數(shù)86. Newton插值函數(shù)87. Newton插值后的求解函數(shù)98. 三樣條插值函數(shù)99. 三樣條的查找函數(shù)910. 三樣條計(jì)算函數(shù)911. 積分

2、方程化為線性方程組的函數(shù)1012. 方程求解的腳本1013. Romberg函數(shù)1014. 二分法函數(shù)1115. 割線法函數(shù)113.結(jié)果分析121. 列主元Gauss消去122. 矩陣分解123. 追趕法124. 迭代法求解非線性方程組125. Newton插值和三樣條插值136. Simpson公式求解積分方程167. Romberg方法求積分178. 二分法和割線法174 心得體會(huì)175 主要參考資料181 作業(yè)內(nèi)容2 解題思路和算法組織解題思路1. Gauss消去法 在Matlab中利用列主元的Gauss消去法求解線性方程組，尋找主元的步驟可以借助find函數(shù)找到主元的索引坐標(biāo)，行交換步

3、驟在Matlab中可以借助一個(gè)中間矩陣變量方便的實(shí)現(xiàn)整行交換，同時(shí)列向量中的元素也應(yīng)交換。消去過(guò)程結(jié)束后，即可回帶解得解向量。值得注意的是，Matlab的索引是從（1,1）開(kāi)始的，因此程序中應(yīng)注意循環(huán)變量的控制。另外一點(diǎn)是，為了驗(yàn)證結(jié)果的正確性，需要提前保留系數(shù)矩陣。2. LDL和Cholesky分解對(duì)于LDL和Cholesky分解，教材上已經(jīng)給出了詳細(xì)的偽代碼，只需要將其轉(zhuǎn)換成相應(yīng)的Matlab語(yǔ)言即可，需要注意Matlab的索引從（1,1）開(kāi)始，在所要分解的矩陣生成方面，采用兩步循環(huán)即可。由于Matlab的函數(shù)M文件利于模塊編程，因此分別建立了LDL和Cholesky分解的子函數(shù) LD和C

4、holesky。事實(shí)上，對(duì)于題目中的待分解矩陣，LDL和Cholesky分解的結(jié)果是相同的。3. 追趕法三帶寬矩陣線性方程組的求解不需要像Gauss消去法那么繁瑣，TSS方法只需要三個(gè)已知系數(shù)向量和一個(gè)增廣向量即可求解，在Matlab中采用橫向追趕法求解。編程時(shí)，根據(jù)教材的詳盡的偽代碼編譯成Matlab語(yǔ)言即可。4. 迭代法求解線性方程組迭代法的Matlab程序也可以根據(jù)教材的偽代碼得出，此時(shí)需要注意D，E，F(xiàn)三個(gè)矩陣的得出，Matlab的diag指令可以提取矩陣的對(duì)角元素，得到D矩陣，再利用一個(gè)循環(huán)，消去上三角元素，即可得到E矩陣，計(jì)算A-D-E即可得到F矩陣。在迭代過(guò)程中，設(shè)置一個(gè)迭代次數(shù)

5、的計(jì)數(shù)變量，即可比較Jacobi迭代和Jacobi-sediel迭代的速度。由于兩種迭代方法均編寫(xiě)成函數(shù)形式，初始向量為一個(gè)參量，因此可以通過(guò)改變初始向量的值，觀察初始向量對(duì)迭代的影響。5. 牛頓插值和三樣條插值的比較 首先編輯一個(gè)Matlab程序，獲得Newton插值的差商表，將多項(xiàng)式所用的系數(shù)儲(chǔ)存于一個(gè)向量中，再編寫(xiě)一個(gè)使用Newton插值多項(xiàng)式計(jì)算非插值點(diǎn)值的計(jì)算程序（使用秦九韶算法）；三樣條插值也是如此，需要注意的是兩點(diǎn)，一是M矩陣邊界處的處理，實(shí)際上，本題中所逼近的函數(shù)已經(jīng)知道，因此兩個(gè)邊界點(diǎn)的二階導(dǎo)數(shù)值可以直接求出。二是需要編寫(xiě)一個(gè)查找程序，求得非插值點(diǎn)的位置，然后求解其值。為了運(yùn)

6、行一個(gè)腳本文件分別計(jì)算出n=5,10,20時(shí)的情況，可以設(shè)定一個(gè)循環(huán)，循環(huán)變量i=1:3，則n就可以表示為5*2(i-1).6. Simpson公式用于求解積分方程對(duì)于本題，Simpson積分公式設(shè)計(jì)用到九個(gè)節(jié)點(diǎn)，據(jù)此運(yùn)用復(fù)化Simpson公式可以得到一個(gè)9*9的線性方程組。此時(shí)系數(shù)矩陣和已知向量的得出編寫(xiě)一個(gè)名為f71的程序得出 。再運(yùn)用Jacobi-sediel迭代方法求解該方程組，得到9個(gè)節(jié)點(diǎn)處的值，運(yùn)用三次樣條插值就可以得到近似的解函數(shù)。最后調(diào)用三次樣條的計(jì)算函數(shù)，比較其與真解的誤差。 7. Romberg方法計(jì)算定積分 本題擬編寫(xiě)一個(gè)函數(shù)M文件，計(jì)算出積分結(jié)果，并統(tǒng)計(jì)Romberg方

7、法達(dá)到相應(yīng)精度時(shí)的運(yùn)算步數(shù)，因此誤差設(shè)定為參量。教材上的偽代碼用Matlab實(shí)現(xiàn)時(shí)需要設(shè)定一個(gè)死循環(huán)，循環(huán)里面設(shè)定跳出循環(huán)的條件。8. 二分法和割線法 本題中積分函數(shù)的獲得可以運(yùn)用Matlab的符號(hào)運(yùn)算功能和Int函數(shù)獲得，運(yùn)算subs函數(shù)計(jì)算其在端點(diǎn)處的值，算法實(shí)現(xiàn)上，二分法在教材上的偽代碼可以通過(guò)Matlab方便的實(shí)現(xiàn)，并將其編寫(xiě)成一個(gè)函數(shù)M文件，其誤差參數(shù)可調(diào)，（這里設(shè)定的誤差比所要求的誤差稍微大一點(diǎn)）并且得到其解所在的上下界，在命令窗中顯示。割線法的偽代碼實(shí)現(xiàn)在Matlab中實(shí)現(xiàn)也是比較容易的，只不過(guò)此時(shí)需要用戶設(shè)定一個(gè)循環(huán)次數(shù)上界，用于規(guī)定循環(huán)的最大次數(shù)。命令窗中得到了解的上下界，則

8、可以選定這兩個(gè)界限作為割線法的迭代起始值，并且設(shè)定誤差值，即可得到所求方程的解。程序清單(雙欄)1. Gauss消去法18format longt=rand(3,3);b=1e-8,2,3;-1,3.712,4.623;-2,1.072,5.643;A=b;m=1,2,3;for k=1:2u,p=find(abs(b(k:end,k)=max(abs(b(k:end,k);%查找主元下標(biāo)u=u+k-1;if b(u,k)=0 disp('error! It never be 0!')endt(k,:)=b(k,:);b(k,:)=b(u,:);b(u,:)=t(k,:);l=

9、m(k);m(k)=m(u);m(u)=l;%行交換for i=k+1:3 b(i,k)=b(i,k)/b(k,k); for j=k+1:3 b(i,j)=b(i,j)-b(i,k)*b(k,j); end m(i)=m(i)-b(i,k)*m(k);endendx(3)=m(3)/b(3,3);for k=2:-1:1 x(k)=(m(k)-sum(b(k,k+1:end).*x(k+1:end)/b(k,k)end%回帶求解x2. LDL分解function L,D=LD(n)r=;for i=1:n for j=1:n if i-j<0; A(i,j)=i; else A(i,j

10、)=j; end endend%生成待分解矩陣for k=1:n for p=1:k-1 r(p)=A(p,p)*A(k,p); end A(k,k)=A(k,k)-sum(A(k,1:k-1).*r(1:k-1); if A(k,k)=0 disp('stop fatal error!'); else for i=k+1:n A(i,k)=(A(i,k)-sum(A(i,1:k-1).*r(1:k-1)/A(k,k); end endendD=diag(diag(A);%求取對(duì)角矩陣for i=1:n A(i,i+1:end)=0;endL=A;3. cholesky分解fu

11、nction G=cholesky(n)r=;for i=1:n for j=1:n if i-j<0; A(i,j)=i; else A(i,j)=j; end endendA(1,1)=sqrt(A(1,1);for k=1:n-1 for i=k+1:n A(i,k)=(A(i,k)-sum(A(i,1:k-1).*A(k,1:k-1)/A(k,k); end A(k+1,k+1)=sqrt(A(k+1,k+1)-sum(A(k+1,1:k).2);endfor i=1:n A(i,i+1:end)=0;endG=A;4. Jacobi迭代function x,N=Jacobi(A

12、,b,x0,e,n)x1=;N=0;D=diag(diag(A);Ac=A;x0c=x0;for i=1:n A(i,i:end)=0;endE=A;F=Ac-A-D;E=-E;F=-F;%獲得D,E,F矩陣。G=inv(D)*(E+F);d=inv(D)*b;q=norm(G);%二次范數(shù)x0=G*x0+d;fanshu=norm(x0-x0c,inf);Nx=log(1-q)*e/fanshu)/log(q);Nx=Nx+2;%考慮其他因素，加大迭代上限for k=2:Nx x1=x0; x0=G*x0+d; if norm(x0-x1,inf)<e*(1-q) k=Nx+1; en

13、d N=N+1;endx=x1;5. Jacobi-sediel迭代函數(shù) function x,N=Jacobi-sediel(A,b,x0,e,n)% Jacobi-sediel method % edited by shao% for special usex1=;N=0;D=diag(diag(A)Ac=A;x0c=x0;for i=1:n A(i,i:end)=0;endE=A;F=Ac-A-D;E=-E;F=-F;G=inv(D-E)*F;d=inv(D-E)*b;q=norm(G);x0=G*x0+d;fanshu=norm(x0-x0c,inf);Nx=log(1-q)*e/fa

14、nshu)/log(q);Nx=Nx+2;for k=2:Nx x1=x0; x0=G*x0+d; if norm(x0-x1,inf)<e*(1-q) k=Nx+1; end N=N+1;endx=x1;6. Newton插值函數(shù)function y=Newton(xi,yi,n)y=yi;for k=1:n for i=n+1:-1:k+1 y(i)=(y(i)-y(i-1)/(xi(i)-xi(i-k); endend 7. Newton插值后的求解函數(shù) function y=Ncal(x,x1,y1,m,n)for k=1:my(k)=y1(n+1);for i=n:-1:1 y

15、(k)=y(k)*(x(k)-x1(i)+y1(i);endend8. 三樣條插值函數(shù) function M=Splineself(x,y,n) M=y; for k=1:2 for i=n+1:-1:k+1 M(i)=(M(i)-M(i-1)/(x(i)-x(i-k); end end h=ones(1,n); c=ones(1,n); a=ones(1,n); b=2*ones(1,n+1); h(1)=x(2)-x(1); for i=2:n h(i)=x(i+1)-x(i); c(i)=h(i)/(h(i)+h(i-1); a(i-1)=1-c(i); M(i)=6*M(i+1); e

16、nd duo=diff('1/(1+25*x2)',2); M(1)=2*subs(duo,-1); M(n+1)=2*subs(duo,1); c(1)=0; a(n)=0; A=diag(b)+diag(c,1)+diag(a,-1); M=Tss(A,M,n+1);9. 三樣條的查找函數(shù)function K=FindK(x,xs,n)K=1;for i=1:n if xs<=x(i+1) K=i; i=n+1; else K=i+1; endend10. 三樣條計(jì)算函數(shù)function Y=cacul2(xk,xi,yi,M,N,n)for i=1:N K=Find

17、K(xi,xk(i),n); h=xi(K+1)-xi(K); xbal=xi(K+1)-xk(i); xtao=xk(i)-xi(K); Y(i)=(M(K+1)*xbal3/6+M(K)*xtao3/6+(yi(K)-M(K)*h2/6)*xbal. +(yi(K+1)-M(K+1)*h2/6)*xtao)/h;end11. 積分方程化為線性方程組的函數(shù)function y=f71(t,x,n)y=;for i=1:ny=y;1./(1+t)-x(i);endy=y*diag(1/6*1/8,4/6*1/8,2/6*1/8,4/6*1/8,2/6*1/8,4/6*1/8,2/6*1/8,4

18、/6*1/8,1/6*1/8);y1=eye(9);y=y-y1;12. 方程求解的腳本clc;x=0:0.125:1;t=x;b=-f7(x);A=f71(t,x,9);y,N=Jacobi(A,b',ones(9,1),1e-7,9);yNM=Splineself(x,y,8);xk=0:1/2:1;Y=cacul2(xk,x,y,M,3,8);Yc=1./(xk+1).2;e=Y-Ycem=max(abs(e)13. Romberg函數(shù)function Result,K=Romberg(e)h=1;k=2;T(1)=(h/2)*(f6(1)+f6(0);m=true;while

19、m=trueS=0;x=h/2;S=S+f6(x);x=x+h;while x<1 S=S+f6(x); x=x+h; endT(2)=0.5*(T(1)+h*S);s(2)=T(2)+(1/3)*(T(2)-T(1);h=h/2;T(1)=T(2);s(1)=s(2);if k>2C(2)=s(2)+(1/15)*(s(2)-s(1);C(1)=C(2);endif k>3R(2)=C(2)+(1/63)*(C(2)-C(1);if k>4 if abs(R(2)-R(1)<e m=false; endendR(1)=R(2);endk=k+1;endResul

20、t=R(2);K=k-1;14. 二分法函數(shù)function result,x2,y2=bisection(E)syms x y;I=int('cos(y*cos(x)','x',0,pi);I=I/pi;x0=2;x1=3;f0=subs(I,x0);f1=subs(I,x1);if f0*f1>0 disp('error!')ends=true;while s if abs(f0)<1e-9 result=x0; break end if abs(f1)<1e-9 result=x1; break end x=0.5*(x0

21、+x1); if abs(x1-x)<E result=x; break end f=subs(I,x); if abs(f)<1e-9 result=x; break end if f1*f<0 x0=x; f0=f; else x1=x; f1=f; endendx2=x0;y2=x1;ft=subs(I,result); if ft*f1<0 fprintf('the section is %.8f to %.8f',result,x1) else fprintf('the section is %.8f to %.8f',x0,r

22、esult) end15. 割線法函數(shù)function result=secant(x0,x1,E1,E2,n)syms x y;I=int('cos(y*cos(x)','x',0,pi);I=I/pi;f0=subs(I,x0);f1=subs(I,x1);if f1>f0 t=x1;x1=x0;x0=t; t=f1;f1=f0;f0=t;endfor k=1:n if abs(f1)<E2 result=x1; break end s=f1/f0; x=x1-(x0-x1)*s/(1-s); f=subs(I,x); if abs(x1-x)&

23、lt;E1 result=x; break end if abs(f)>abs(f1) x0=x;f0=f; else x0=x1;f0=f1; x1=x;f1=f; endendif k=nfprintf('sorry,%.2g times iteration can not get a convergent result',n);end3.結(jié)果分析 1. 列主元Gauss消去運(yùn)行method.m得出結(jié)果x= -0.491058221221525 -0.050886077442433 0.367257386598483。在命令窗中驗(yàn)證：A*x= 1.0000000000

24、00000 2.000000000000000 3.000000000000000可知所求得的解確實(shí)為題目要求的解2. 矩陣分解運(yùn)行腳本second.m。LDL結(jié)果如下：L形如000100110，但是為一個(gè)20階的矩陣，D為20階的單位陣。Cholesky分解結(jié)果為：G形如100110111，為一個(gè)20階的矩陣。3. 追趕法運(yùn)行腳本文件Third.m。運(yùn)行結(jié)果為x=1,2,3,4,5。4. 迭代法求解非線性方程組運(yùn)行腳本文件fourth.m。Jacobi迭代運(yùn)行結(jié)果為x = 0.856984839849005 -0.032615489411454 -1.795417505967341 0.88

25、9601489821177。N=7。Jacobi-sediel迭代為x= 0.856975730668451 -0.032619572404447 -1.795416561900068 0.889595204980871。N=6?？梢?jiàn)本題中Jacobi-sediel迭代次數(shù)比Jacobi迭代少一次。5. Newton插值和三樣條插值運(yùn)行script1.m和script2.m文件，可以得到如下結(jié)果：Newton插值結(jié)果：n=5時(shí)：xi= -1.0000 -0.6000 -0.2000 0.2000 0.6000 1.0000，f(xi)= 0.0385 0.1000 0.5000 0.5000

26、0.1000 0.0385,差商向量y= 0.038461538461538 0.153846153846154 1.057692307692308 -1.923076923076925 1.201923076923080 -0.000000000000004所求非插值數(shù)據(jù)點(diǎn)的值為Y= 0.0137 -0.0069 -0.0236 -0.0366 -0.0460 -0.0522 -0.0553 -0.0555 -0.0530 -0.0481 -0.0481 -0.0530 -0.0555 -0.0553 -0.0522 -0.0460 -0.0366 -0.0236 -0.0069 0.013

27、7，最大誤差e= 0.1092n=10時(shí)：xi= -1.0000 -0.8000 -0.6000 -0.4000 -0.2000 0 0.2000 0.4000 0.6000 0.8000 1.0000, f(xi)= 0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.0385,y=100* 0.000384615384615 0.001018099547511 0.002601809954751 0.007918552036199 0.026866515837104 -0.0636312217194

28、57 -0.176753393665158 0.848416289592759 -1.679157239819003 2.209417420814478 -2.209417420814479,Y= 1.2303 1.8044 1.9590 1.8458 1.5787 1.2402 0.8881 0.5604 0.2802 0.0588 0.0588 0.2802 0.5604 0.8881 1.2402 1.5787 1.8458 1.9590 1.8044 1.2303， e= 1.9156n=20時(shí)：xi= -1.0000 -0.9000 -0.8000 -0.7000 -0.6000 -

29、0.5000 -0.4000 -0.3000 -0.2000 -0.1000 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000, f(xi)= 0.0385 0.0471 0.0588 0.0755 0.1000 0.1379 0.2000 0.3077 0.5000 0.8000 1.0000 0.8000 0.5000 0.3077 0.2000 0.1379 0.1000 0.0755 0.0588 0.0471 0.0385,y=100000*0.000000384615385 0.00000

30、0859728507 0.000001583710407 0.000002860070008 0.000005335951507 0.000010377505689 0.000020019018067 0.000027967745829 -0.000075439314408 -0.001011990803028 -0.000643994147381 0.021527804355314 -0.072679339490163 0.113937426075134 -0.035384293812155 -0.283074350497223 0.866915198397751 -1.5922932215

31、46891 2.237536226356743 -2.601786309717143 2.601786309717144,Y=-58.2381 -50.8644 -28.6626 -10.3345 0.0471 3.9451 4.0691 2.6744 1.1353 0.0588 0.0588 1.1353 2.6744 4.0691 3.9451 0.0471 -10.3345 -28.6626 -50.8644 -58.2381，e=58.2781三樣條插值結(jié)果：n=5時(shí)：二階導(dǎo)數(shù)矩陣M= 0.2105 4.0742 -3.8148 -3.8148 4.0742 0.2105,Y= -6.

32、1874 -5.9547 -5.7279 -5.5070 -5.2919 -5.0825 -4.8787 -4.6805 -4.4878 -4.3004 0.0264 0.0328 0.0401 0.0484 0.0578 0.0684 0.0802 0.0934 0.1079 0.1239,最大誤差e=6.2274n=10時(shí) M= 0.2105 0.3536 1.4971 2.4816 18.5767 -46.7883 18.5767 2.4816 1.4971 0.3536 0.2105 ,Y= -0.0711 -0.0585 -0.0463 -0.0345 -0.0232 -0.0122

33、 -0.0017 0.0085 0.0183 0.0277 0.0579 0.0556 0.0533 0.0512 0.0492 0.0472 0.0454 0.0437 0.0422 0.0407 ,e=0.1111n=20時(shí) M= 0.2105 0.3051 0.4694 0.7474 1.2692 2.2174 4.3439 7.7809 15.3016 -4.3719 -57.8141 -4.3719 15.3016 7.7809 4.3439 2.2174 1.2692 0.7474 0.4694 0.3051 0.2105 ,Y=-0.4505 -0.4272 -0.4045 -0.3824 -0.3609 -0.3400 -0.31