计算方法报告程序

1、计算方法作业目录1 作业内容32 解题思路和算法组织4解题思路41. Gauss消去法42. LDL和Cholesky分解43. 追赶法54. 迭代法求解线性方程组55. 牛顿插值和三样条插值的比较56. Simpson公式用于求解积分方程57. Romberg方法计算定积分58. 二分法和割线法6程序清单(双栏)61. Gauss消去法62. LDL分解73. cholesky分解74. Jacobi迭代75. Jacobi-sediel迭代函数86. Newton插值函数87. Newton插值后的求解函数98. 三样条插值函数99. 三样条的查找函数910. 三样条计算函数911. 积分

2、方程化为线性方程组的函数1012. 方程求解的脚本1013. Romberg函数1014. 二分法函数1115. 割线法函数113.结果分析121. 列主元Gauss消去122. 矩阵分解123. 追赶法124. 迭代法求解非线性方程组125. Newton插值和三样条插值136. Simpson公式求解积分方程167. Romberg方法求积分178. 二分法和割线法174 心得体会175 主要参考资料181 作业内容2 解题思路和算法组织解题思路1. Gauss消去法 在Matlab中利用列主元的Gauss消去法求解线性方程组，寻找主元的步骤可以借助find函数找到主元的索引坐标，行交换步

3、骤在Matlab中可以借助一个中间矩阵变量方便的实现整行交换，同时列向量中的元素也应交换。消去过程结束后，即可回带解得解向量。值得注意的是，Matlab的索引是从（1,1）开始的，因此程序中应注意循环变量的控制。另外一点是，为了验证结果的正确性，需要提前保留系数矩阵。2. LDL和Cholesky分解对于LDL和Cholesky分解，教材上已经给出了详细的伪代码，只需要将其转换成相应的Matlab语言即可，需要注意Matlab的索引从（1,1）开始，在所要分解的矩阵生成方面，采用两步循环即可。由于Matlab的函数M文件利于模块编程，因此分别建立了LDL和Cholesky分解的子函数 LD和C

4、holesky。事实上，对于题目中的待分解矩阵，LDL和Cholesky分解的结果是相同的。3. 追赶法三带宽矩阵线性方程组的求解不需要像Gauss消去法那么繁琐，TSS方法只需要三个已知系数向量和一个增广向量即可求解，在Matlab中采用横向追赶法求解。编程时，根据教材的详尽的伪代码编译成Matlab语言即可。4. 迭代法求解线性方程组迭代法的Matlab程序也可以根据教材的伪代码得出，此时需要注意D，E，F三个矩阵的得出，Matlab的diag指令可以提取矩阵的对角元素，得到D矩阵，再利用一个循环，消去上三角元素，即可得到E矩阵，计算A-D-E即可得到F矩阵。在迭代过程中，设置一个迭代次数

5、的计数变量，即可比较Jacobi迭代和Jacobi-sediel迭代的速度。由于两种迭代方法均编写成函数形式，初始向量为一个参量，因此可以通过改变初始向量的值，观察初始向量对迭代的影响。5. 牛顿插值和三样条插值的比较 首先编辑一个Matlab程序，获得Newton插值的差商表，将多项式所用的系数储存于一个向量中，再编写一个使用Newton插值多项式计算非插值点值的计算程序（使用秦九韶算法）；三样条插值也是如此，需要注意的是两点，一是M矩阵边界处的处理，实际上，本题中所逼近的函数已经知道，因此两个边界点的二阶导数值可以直接求出。二是需要编写一个查找程序，求得非插值点的位置，然后求解其值。为了运

6、行一个脚本文件分别计算出n=5,10,20时的情况，可以设定一个循环，循环变量i=1:3，则n就可以表示为5*2(i-1).6. Simpson公式用于求解积分方程对于本题，Simpson积分公式设计用到九个节点，据此运用复化Simpson公式可以得到一个9*9的线性方程组。此时系数矩阵和已知向量的得出编写一个名为f71的程序得出 。再运用Jacobi-sediel迭代方法求解该方程组，得到9个节点处的值，运用三次样条插值就可以得到近似的解函数。最后调用三次样条的计算函数，比较其与真解的误差。 7. Romberg方法计算定积分 本题拟编写一个函数M文件，计算出积分结果，并统计Romberg方

7、法达到相应精度时的运算步数，因此误差设定为参量。教材上的伪代码用Matlab实现时需要设定一个死循环，循环里面设定跳出循环的条件。8. 二分法和割线法 本题中积分函数的获得可以运用Matlab的符号运算功能和Int函数获得，运算subs函数计算其在端点处的值，算法实现上，二分法在教材上的伪代码可以通过Matlab方便的实现，并将其编写成一个函数M文件，其误差参数可调，（这里设定的误差比所要求的误差稍微大一点）并且得到其解所在的上下界，在命令窗中显示。割线法的伪代码实现在Matlab中实现也是比较容易的，只不过此时需要用户设定一个循环次数上界，用于规定循环的最大次数。命令窗中得到了解的上下界，则

8、可以选定这两个界限作为割线法的迭代起始值，并且设定误差值，即可得到所求方程的解。程序清单(双栏)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);%查找主元下标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);%求取对角矩阵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);%二次范数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迭代函数 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插值函数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插值后的求解函数 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. 三样条插值函数 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. 三样条的查找函数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. 三样条计算函数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. 积分方程化为线性方程组的函数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函数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. 二分法函数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. 割线法函数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.结果分析 1. 列主元Gauss消去运行method.m得出结果x= -0.491058221221525 -0.050886077442433 0.367257386598483。在命令窗中验证：A*x= 1.0000000000

24、00000 2.000000000000000 3.000000000000000可知所求得的解确实为题目要求的解2. 矩阵分解运行脚本second.m。LDL结果如下：L形如000100110，但是为一个20阶的矩阵，D为20阶的单位阵。Cholesky分解结果为：G形如100110111，为一个20阶的矩阵。3. 追赶法运行脚本文件Third.m。运行结果为x=1,2,3,4,5。4. 迭代法求解非线性方程组运行脚本文件fourth.m。Jacobi迭代运行结果为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。可见本题中Jacobi-sediel迭代次数比Jacobi迭代少一次。5. Newton插值和三样条插值运行script1.m和script2.m文件，可以得到如下结果：Newton插值结果：n=5时：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所求非插值数据点的值为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时：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时：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三样条插值结果：n=5时：二阶导数矩阵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时 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时 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