矩阵上机作业

1、矩阵与数值分析课程数值实验报告 1. 方程在x=3.0附近有根，试写出其三种不同的等价形式以构成两种不同的迭代格式，再用这两种迭代求根，并绘制误差下降曲线，观察这两种迭代是否收敛及收敛的快慢newton迭代法function x=juz1(x1,tol,max)xo=x1;T=;i=1;while imax xn=xo-(2*xo3-5*xo2-19*xo+42)/(6*xo2-10*xo-19); T(i)=xn; if abs(xn-xo) juz1(3,10-5,20)迭代法收敛 3.5000i = 6弦截法function x=juz11(x1,x2,tol,max)xo=x1;xp=

2、x2;i=1;while imax xn=xp-(2*xp3-5*xp2-19*xp+42)/(2*xp2+2*xp*xo+2*xo2-5*xp-5*xo-19); T(i)=xn; if abs(xn-xp) juz11(2.9,3,10-5,20)迭代法收敛 3.5000i = 82. 用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分，要求绝对误差为，并将计算结果与精确解进行比较：（1） （2）（1） 用复化梯形公式求积分function f=Txing1(a,b,n)% a，b分别为积分下限和上限，n为划分的区间个数while n104 h1=(b-a)/n;%h1为积分区间长度

3、x1=zeros(1,n+1); for i=1:n+1 x1(i)=a+(i-1)*h1;%计算n+1个节点储存在x1中 end y1=(2/3)*(x1.3).*exp(x1.2);%计算每个节点对应的函数值 to=0; for i=1:n to=to+h1/2*(y1(i)+y1(i+1);%用梯形公式求积分 end h2=(b-a)/(n+1);%h2为积分区间长度 x2=zeros(1,n+2); for i=1:n+2 x2(i)=a+(i-1)*h2;%计算n+2个节点储存在x2中 end y2=(2/3)*(x2.3).*exp(x2.2);%计算每个节点对应的函数值 tn=0

4、; for i=1:n+1 tn=tn+h2/2*(y2(i)+y2(i+1);%用梯形公式求积分 end if abs(tn-to) Txing1(1,2,5)I = 54.5982N = 5121r = -5.0604e-06复化辛普森求积分function f=Simpson1(a,b,n)%a,b分别为积分下限和上限%n为分割区间个数while n108T=;for n=n:n+1h1=(b-a)/n;%将积分区间n等分x1=zeros(1,n+1);for i=1:n+1 x1(i)=a+(i-1)*h1;%计算n+1个求积节点值并储存在x中endz1=(2/3)*x1.3.*exp

5、(x1.2);%计算n+1个求积节点处的函数值y1=zeros(1,n);for i=1:n y1(i)=a+(2*i-1)*h1/2;%计算n个区间中点值endr1=(2/3)*y1.3.*exp(y1.2);%计算n个区间中点对应的函数值t=0;for i=1:n; t=t+h1/6*(z1(i)+4*r1(i)+z1(i+1);%用复化simpson公式求积分endT=T,t;endif abs(T(2)-T(1) Simpson1(1,2,3)I = 54.5982N = 160r = -2.9586e-08龙贝格公式求积分function f=romberg1(a,b)t(1,1)=

6、(b-a)/2*(2*a3*exp(a2)/3+2*b3*exp(b2)/3);n=1;while n100for k=1:n aa=0; for i=0:(power(2,(k-1)-1); aa=aa+2*(a+(2*i+1)*(b-a)/2k)3*exp(a+(2*i+1)*(b-a)/2k)2)/3; end t(1,k+1)=0.5*t(1,k)+(b-a)/(2k)*aa;endfor m=1:n h=1/(4m)-1); for k=m:n t(m+1,k+1)=(4m)*t(m,k+1)-t(m,k)*h; endend if (abs(t(n,n)-t(n+1,1+n) ro

7、mberg1(1,2)I = 54.5982（2）复化梯形公式求积分function f=Txing2(a,b,n)% a，b分别为积分下限和上限，n为划分的区间个数while n10000 h1=(b-a)/n;%h为积分区间长度 x1=zeros(1,n+1); for i=1:n+1 x1(i)=a+(i-1)*h1;%计算n+1个节点储存在x1中 end y1=2*x1./(x1.2-3);%计算每个节点对应的函数值 to=0; for i=1:n to=to+h1/2*(y1(i)+y1(i+1);%用梯形公式求积分 end h2=(b-a)/(n+1);%h2为积分区间长度 x2=

8、zeros(1,n+2); for i=1:n+2 x2(i)=a+(i-1)*h2;%计算n+2个节点储存在x2中 end y2=2*x2./(x2.2-3);%计算每个节点对应的函数值 tn=0; for i=1:n+1 tn=tn+h2/2*(y2(i)+y2(i+1);%用梯形公式求积分 end if abs(tn-to) Txing2(2,3,4)I = 1.7918N = 1025r = -1.0576e-06复化辛普森公式求积分function f=Simpson2(a,b,n)%a,b分别为积分下限和上限%n为分割区间个数while n10000T=;%生成空矩阵以储存计算结果

9、for n=n:n+1h=(b-a)/(2*n);%将积分区间分成2n份x=zeros(1,2*n+1);for i=1:2*n+1 x(i)=a+(i-1)*h;%计算n+1个求积节点值并储存在x中endz=2*x./(x.2-3);%计算2n+1个求积节点处的函数值t=0;for i=1:n; t=t+h/3*(z(2*i-1)+4*z(2*i)+z(2*i+1);%用复化simpson公式求积分endT=T,t;%将分成n份和n+1份结果储存在T中endif abs(T(2)-T(1) Simpson2(2,3,4)I = 1.7918N = 96r = -4.9476e-09龙贝格公示

10、求积分function f=romberg2(a,b)t(1,1)=(b-a)/2*(2*a/(a2-3)+2*b/(b2-3);n=1;while n100for k=1:n aa=0; for i=0:(power(2,(k-1)-1); aa=aa+2*(a+(2*i+1)*(b-a)/2k)/(a+(2*i+1)*(b-a)/2k)2-3); end t(1,k+1)=0.5*t(1,k)+(b-a)/(2k)*aa;endfor m=1:n h=1/(4m)-1); for k=m:n t(m+1,k+1)=(4m)*t(m,k+1)-t(m,k)*h; endend if (abs

11、(t(n,n)-t(n+1,1+n) romberg2(2,3)I = 1.79183. 使用带选主元的分解法求解线性方程组,其中，当时对于的情况分别求解精确解为对得到的结果与精确解的差异进行解释function f=LU(n)%生成矩阵AA=zeros(n);for i=1:n for j=1:n A(i,j)=i(j-1); endend%生成右边项bb=zeros(n,1);b(1)=n;for i=2:n b(i)=(in-1)/(i-1);endx=zeros(n,1);y=zeros(n,1);a=zeros(n,1);%储存矩阵A的换行元素c=0;%储存右边项的换行元素for i

12、=1:n-1 e,m=max(abs(A(i:n,i);%寻找最大元素及其位置 a=A(i,:); A(i,:)=A(m+i-1,:); A(m+i-1,:)=a; c=b(i); b(i)=b(m+i-1); b(m+i-1)=c; if A(i,i)=0 disp(A是奇异矩阵，方程组没有唯一解); return; end for j=i+1:n d=A(j,i)/A(i,i); A(j,i)=d; A(j,i+1:n)=A(j,i+1:n)-d*A(i,i+1:n); endend%求解Ly=by(1)=b(1);for k=2:n y(k)=b(k)-A(k,1:k-1)*y(1:k-

13、1);end%求解Ux=yx(n)=y(n)/A(n,n);for k=n-1:-1:1 x(k)=(y(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endx LU(3)x = 1 1 1 LU(7)x = 1 1 1 1 1 1 1 LU(11)x = 1.0001 0.9997 1.0004 0.9998 1.0001 1.0000 1.0000 1.0000 1.0000 1.0000 1.00004. 用4阶Runge-kutta法求解微分方程(1) 令，使用上述程序执行20步，然后令，使用上述程序执行40步(2) 比较两个近似解与精确解(3) 当减半时，(1)中的最终

14、全局误差是否和预期相符？(4) 在同一坐标系上画出两个近似解与精确解（提示输出矩阵包含近似解的和坐标，用命令plot(R(:,1),R(:,2)画出相应图形）function f=Rungek(x,t)uo=x;tn=t;%h=0.1,计算20步，并将计算结果储存在U1中U1=zeros(20,1);h1=0.1;i=1;while i=20 k1=exp(-2*tn)-2*uo; k2=exp(-2*(tn+0.5*h1)-2*(uo+0.5*h1*k1); k3=exp(-2*(tn+0.5*h1)-2*(uo+0.5*h1*k2); k4=exp(-2*(tn+h1)-2*(uo+h1*

15、k3); un=uo+h1/6*(k1+2*k2+2*k3+k4); U1(i)=un; i=i+1; tn=(i-1)*h1; uo=un;end%h=0.05,计算40步，并将计算结果储存在U2中U3=zeros(40,1);h2=0.05;i=1;uo=x;tn=t;while i Rungek(0.1,0)U1 = 0.1637 0.2011 0.2195 0.2247 0.2207 0.2108 0.1973 0.1817 0.1653 0.1489 0.1330 0.1179 0.1040 0.0912 0.0797 0.0693 0.0601 0.0519 0.0447 0.03

16、85U2 = 0.1637 0.2011 0.2195 0.2247 0.2207 0.2108 0.1973 0.1817 0.1653 0.1489 0.1330 0.1179 0.1040 0.0912 0.0797 0.0693 0.0601 0.0519 0.0447 0.0385U = 0.1637 0.2011 0.2195 0.2247 0.2207 0.2108 0.1973 0.1817 0.1653 0.1489 0.1330 0.1179 0.1040 0.0912 0.0797 0.0693 0.0601 0.0519 0.0447 0.0385R1 = 1.0e-0

17、05 * 0.2114 0.3250 0.3732 0.3791 0.3590 0.3242 0.2825 0.2389 0.1966 0.1575 0.1226 0.0924 0.0667 0.0454 0.0281 0.0142 0.0034 -0.0048 -0.0108 -0.0151R2 = 1.0e-006 * 0.1191 0.1830 0.2098 0.2127 0.2010 0.1811 0.1574 0.1326 0.1087 0.0866 0.0670 0.0499 0.0356 0.0236 0.0140 0.0063 0.0003 -0.0042 -0.0075 -0

18、.0097R3 = 1.0e-005 * -0.1995 -0.3067 -0.3522 -0.3578 -0.3389 -0.3061 -0.2667 -0.2256 -0.1857 -0.1488 -0.1159 -0.0874 -0.0632 -0.0430 -0.0267 -0.0136 -0.0034 0.0044 0.01010.01415. 设为阶的三对角方阵，是一个阶的对称正定矩阵其中为阶单位矩阵。设为线性方程组的真解，右边的向量由这个真解给出。(1) 用Cholesky分解法求解该方程.(2) 用Jacobi迭代法和Gauss-Seidel迭代法求解该方程组，误差设为.其中取

19、值为4，5，6.用Cholesky分解法求解方程：function f=Cholesky(n)%生成Tn矩阵T=zeros(n);for i=1:n T(i,i)=2;endfor i=1:n-1 T(i,i+1)=-1;endfor i=2:n T(i,i-1)=-1;end%生成系数矩阵AI=eye(n);B=T+2*I;C=-I;A=zeros(n2);for i=1:n A(n*(i-1)+1:i*n,n*(i-1)+1:i*n)=B;endfor i=1:n-1 A(i*n+1:(i+1)*n,(i-1)*n+1:i*n)=C; A(i-1)*n+1:i*n,i*n+1:(i+1)*

20、n)=C;end%生成右边项by=ones(n2,1);b=A*y;R=chol(A);%对系数矩阵A进行Cholesky分解x=R(Rb)%求解方程组 Cholesky(4)x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 Cholesky(5)x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.000

21、0 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 Cholesky(6)x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00

22、00 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000Jacobi迭代法求解方程：function f=Jacobi5(x0,n)%生成Tn矩阵T=zeros(n);for i=1:n T(i,i)=2;endfor i=1:n-1 T(i,i+1)=-1;endfor i=2:n T(i,i-1)=-1;end%生成系数矩阵AI=eye(n);B=T+2*I;C=-I;A=zeros(n2);for i=1:n A(n*(i-1)+1:i*n,n*(i-1)+1:i*n)=B;endfor i=1:n-1 A(i*n+1:(i+1)

23、*n,(i-1)*n+1:i*n)=C; A(i-1)*n+1:i*n,i*n+1:(i+1)*n)=C;end%生成右边项by=ones(n2,1);b=A*y;r,s=size(A);xold=x0;C=-A;for i=1:r C(i,i)=0;endfor i=1:r C(i,:)=C(i,:)/A(i,i);endd=zeros(n2,1);for i=1:r d(i)=b(i)/A(i,i);endi=1;while i=500 xnew=C*xold+d; if norm(xnew-xold) a=zeros(16,1); Jacobi5(a,4)Jacobi迭代法收敛x = 1

24、.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000i =87 a=zeros(25,1); Jacobi5(a,5)Jacobi迭代法收敛x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.

25、0000 1.0000 1.0000 1.0000 1.0000i = 126 a=zeros(36,1); Jacobi5(a,6)Jacobi迭代法收敛x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1

26、.0000 1.0000 1.0000 1.0000 1.0000i = 172用Gauss-Seidel迭代法求解function 发=Gauss(xo,n)%生成Tn矩阵T=zeros(n);for i=1:n T(i,i)=2;endfor i=1:n-1 T(i,i+1)=-1;endfor i=2:n T(i,i-1)=-1;end%生成系数矩阵AI=eye(n);B=T+2*I;C=-I;A=zeros(n2);for i=1:n A(n*(i-1)+1:i*n,n*(i-1)+1:i*n)=B;endfor i=1:n-1 A(i*n+1:(i+1)*n,(i-1)*n+1:i*

27、n)=C; A(i-1)*n+1:i*n,i*n+1:(i+1)*n)=C;end%生成右边项by=ones(n2,1);b=A*y;%用Gauss-Seidel法求解C=-A;L=zeros(n2);U=zeros(n2);D=zeros(n2);for i=1:n2 D(i,i)=A(i,i);endfor i=2:n2 L(i,1:i-1)=C(i,1:i-1);endfor i=1:n2-1 U(i,i+1:n2)=C(i,i+1:n2);endB=inv(D-L)*U;f=inv(D-L)*b;i=1;while i500 xn=B*xo+f; if norm(xn-xo) a=zeros(16,1); Gauss(a,4)迭代法收敛x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000i = 46 a=zeros(25,1); Gauss(a,5)迭代法收敛x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0