矩阵数值分析上机作业

《矩阵与数值分析》课程数值实验报告教学班号：任课老师：学部：专业：姓名：学号：1.方程在x=3.0附近有根，试写出其三种不同的等价形式以构成两种不同的迭代格式，再用这两种迭代求根，并绘制误差下降曲线，观察这两种迭代是否收敛及收敛的快慢newton迭代法functionx=juz1(x1,tol,max)%x1为初值，tol为误差终止条件，max为最大迭代次数xo=x1;T=[];%创建T矩阵储存每次迭代的结果i=1;whileimaxxn=xo-(2*xo^3-5*xo^2-19*xo+42)/(6*xo^2-10*xo-19);T(i)=xn;ifabs(xn-xo)toldisp('迭代法收敛');disp([xn]);iR=xn-T;%计算误差a=linspace(0,0.1,i);plot(a,R,'r');%绘制误差下降曲线title('误差下降曲线');return;elsexo=xn;i=i+1;endenddisp('迭代法不收敛');disp('迭代最大次数的结果为：');xnjuz1(3,10^-5,20)迭代法收敛3.5000i=6弦截法functionx=juz11(x1,x2,tol,max)%x1,x2为初值，tol为误差终止条件，max为最大迭代次数xo=x1;xp=x2;i=1;whileimaxxn=xp-(2*xp^3-5*xp^2-19*xp+42)/(2*xp^2+2*xp*xo+2*xo^2-5*xp-5*xo-19);T(i)=xn;ifabs(xn-xp)toldisp('迭代法收敛');disp([xn]);iR=xn-T;a=linspace(0,1,i);plot(a,R,'r');%绘制误差下降曲线title('误差下降曲线');return;elsexo=xp;xp=xn;i=i+1;endenddisp('迭代法不收敛');disp('迭代最大次数的结果为：');xnjuz11(2.9,3,10^-5,20)迭代法收敛3.5000i=82.用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分，要求绝对误差为，并将计算结果与精确解进行比较：（1）（2）（1）用复化梯形公式求积分functionf=Txing1(a,b,n)%a，b分别为积分下限和上限，n为划分的区间个数whilen10^4h1=(b-a)/n;%h1为积分区间长度x1=zeros(1,n+1);fori=1:n+1x1(i)=a+(i-1)*h1;%计算n+1个节点储存在x1中endy1=(2/3)*(x1.^3).*exp(x1.^2);%计算每个节点对应的函数值to=0;fori=1:nto=to+h1/2*(y1(i)+y1(i+1));%用梯形公式求积分endh2=(b-a)/(n+1);%h2为积分区间长度x2=zeros(1,n+2);fori=1:n+2x2(i)=a+(i-1)*h2;%计算n+2个节点储存在x2中endy2=(2/3)*(x2.^3).*exp(x2.^2);%计算每个节点对应的函数值tn=0;fori=1:n+1tn=tn+h2/2*(y2(i)+y2(i+1));%用梯形公式求积分endifabs(tn-to)0.5*10^(-8)%进行误差判断I=tn%计算结果N=n+1%实际分割区间个数r=exp(4)-tn%计算误差return;elsen=n*2;endenddisp('分割区间太长达不到精度要求');n/2+1tnTxing1(1,2,5)I=54.5982N=5121r=-5.0604e-06复化辛普森求积分functionf=Simpson1(a,b,n)%a,b分别为积分下限和上限%n为分割区间个数whilen10^8T=[];forn=n:n+1h1=(b-a)/n;%将积分区间n等分x1=zeros(1,n+1);fori=1:n+1x1(i)=a+(i-1)*h1;%计算n+1个求积节点值并储存在x中endz1=(2/3)*x1.^3.*exp(x1.^2);%计算n+1个求积节点处的函数值y1=zeros(1,n);fori=1:ny1(i)=a+(2*i-1)*h1/2;%计算n个区间中点值endr1=(2/3)*y1.^3.*exp(y1.^2);%计算n个区间中点对应的函数值t=0;fori=1:n;t=t+h1/6*(z1(i)+4*r1(i)+z1(i+1));%用复化simpson公式求积分endT=[T,t];endifabs(T(2)-T(1))0.5*10^(-8)%进行误差判断I=T(2)%计算结果N=n+1%实际分割区间个数r=exp(4)-T(2)%计算误差return;elsen=2*n;endenddisp('分割区间太长达不到精度要求');T(2)n/2+1Simpson1(1,2,3)I=54.5982N=160r=-2.9586e-08龙贝格公式求积分functionf=romberg1(a,b)t(1,1)=(b-a)/2*(2*a^3*exp(a^2)/3+2*b^3*exp(b^2)/3);n=1;whilen100fork=1:naa=0;fori=0:(power(2,(k-1))-1);aa=aa+2*(a+(2*i+1)*(b-a)/2^k)^3*exp((a+(2*i+1)*(b-a)/2^k)^2)/3;endt(1,k+1)=0.5*t(1,k)+(b-a)/(2^k)*aa;endform=1:nh=1/((4^m)-1);fork=m:nt(m+1,k+1)=((4^m)*t(m,k+1)-t(m,k))*h;endendif(abs(t(n,n)-t(n+1,1+n))0.5*10^(-8))I=t(n,n)return;elsen=n+1;endenddisp('超出最大迭代次数');romberg1(1,2)I=54.5982（2）复化梯形公式求积分functionf=Txing2(a,b,n)%a，b分别为积分下限和上限，n为划分的区间个数whilen10000h1=(b-a)/n;%h为积分区间长度x1=zeros(1,n+1);fori=1:n+1x1(i)=a+(i-1)*h1;%计算n+1个节点储存在x1中endy1=2*x1./(x1.^2-3);%计算每个节点对应的函数值to=0;fori=1:nto=to+h1/2*(y1(i)+y1(i+1));%用梯形公式求积分endh2=(b-a)/(n+1);%h2为积分区间长度x2=zeros(1,n+2);fori=1:n+2x2(i)=a+(i-1)*h2;%计算n+2个节点储存在x2中endy2=2*x2./(x2.^2-3);%计算每个节点对应的函数值tn=0;fori=1:n+1tn=tn+h2/2*(y2(i)+y2(i+1));%用梯形公式求积分endifabs(tn-to)0.5*10^(-8)%进行误差判断I=tn%计算结果N=n+1%实际分割区间个数r=log(6)-tn%计算结果与精确值比较return;elsen=n*2;endenddisp('分割区间太长达不到精度要求');n/2+1tnTxing2(2,3,4)I=1.7918N=1025r=-1.0576e-06复化辛普森公式求积分functionf=Simpson2(a,b,n)%a,b分别为积分下限和上限%n为分割区间个数whilen10000T=[];%生成空矩阵以储存计算结果forn=n:n+1h=(b-a)/(2*n);%将积分区间分成2n份x=zeros(1,2*n+1);fori=1:2*n+1x(i)=a+(i-1)*h;%计算n+1个求积节点值并储存在x中endz=2*x./(x.^2-3);%计算2n+1个求积节点处的函数值t=0;fori=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中endifabs(T(2)-T(1))0.5*10^(-8)%进行误差判断I=T(2)%计算结果N=n+1%实际分割区间个数r=log(6)-t%计算误差return;elsen=2*n;endenddisp('分割区间太长达不到精度要求');T(2)n/2+1Simpson2(2,3,4)I=1.7918N=96r=-4.9476e-09龙贝格公示求积分functionf=romberg2(a,b)t(1,1)=(b-a)/2*(2*a/(a^2-3)+2*b/(b^2-3));n=1;whilen100fork=1:naa=0;fori=0:(power(2,(k-1))-1);aa=aa+2*(a+(2*i+1)*(b-a)/2^k)/((a+(2*i+1)*(b-a)/2^k)^2-3);endt(1,k+1)=0.5*t(1,k)+(b-a)/(2^k)*aa;endform=1:nh=1/((4^m)-1);fork=m:nt(m+1,k+1)=((4^m)*t(m,k+1)-t(m,k))*h;endendif(abs(t(n,n)-t(n+1,1+n))0.5*10^(-8))I=t(n,n)return;elsen=n+1;endenddisp('超出最大迭代次数');romberg2(2,3)I=1.79183.使用带选主元的分解法求解线性方程组,其中，，当时．对于的情况分别求解．精确解为．对得到的结果与精确解的差异进行解释．functionf=LU(n)%生成矩阵AA=zeros(n);fori=1:nforj=1:nA(i,j)=i^(j-1);endend%生成右边项bb=zeros(n,1);b(1)=n;fori=2:nb(i)=(i^n-1)/(i-1);endx=zeros(n,1);y=zeros(n,1);a=zeros(n,1);%储存矩阵A的换行元素c=0;%储存右边项的换行元素fori=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;ifA(i,i)==0disp('A是奇异矩阵，方程组没有唯一解');return;endforj=i+1:nd=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);fork=2:ny(k)=b(k)-A(k,1:k-1)*y(1:k-1);end%求解Ux=yx(n)=y(n)/A(n,n);fork=n-1:-1:1x(k)=(y(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);endxLU(3)x=111LU(7)x=1111111LU(11)x=1.00010.99971.00040.99981.00011.00001.00001.00001.00001.00001.0000进行LU分解时，系数矩阵出现病态，产生误差4.用4阶Runge-kutta法求解微分方程tttteetuuuuu222101)(,101)0(,2(1)令1.0h，使用上述程序执行20步，然后令05.0h，使用上述程序执行40步(2)比较两个近似解与精确解(3)当h减半时，(1)中的最终全局误差是否和预期相符？(4)在同一坐标系上画出两个近似解与精确解．（提示输出矩阵R包含近似解的x和y坐标，用命令plot(R(:,1),R(:,2))画出相应图形．）functionf=Rungek(x,t)uo=x;tn=t;%h=0.1,计算20步，并将计算结果储存