《矩阵与数值分析》课程数值实验报告

矩阵与数值分析实验报告学院(系)：建设与工程学部专业：岩土工程学生姓名：陶言祺学号：21206012指导教师：孟兆良完成日期：2012-12-15大连理工大学Dalian

1.给定n阶方程组，其中，则方程组有解。对和，分别用Gauss消去法和列主元消去法解方程组，并比较计算结果。Gauss消去法：Matlab编程（建立GS.m文件）：functionx=GS(n)A=[];b=[];fori=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8;b(i)=15;endA(n,n)=6;b(1)=7;b(n)=14;b=b';fork=1:n-1fori=k+1:nm(i,k)=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endendb(n)=b(n)/A(n,n);fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i);endclearx;x=b;disp('AX=b的解x是')end计算结果：在matlab命令框里输出GS（10）得：>>GS(10)AX=b的解x是ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000在matlab命令框里输出GS（84）得：>>GS(84)AX=b的解x是ans=1.0e+008*0.0000………0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00030.0007-0.00130.0026-0.00520.0105-0.02090.0419-0.08360.1665-0.33030.6501-1.25822.3487-4.02635.3684列主元消去法：Matlab编程（建立GLZX.m文件）：functionx=GLZX(n)A=[];b=[];eps=10^-2;fori=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8;b(i)=15;endA(n,n)=6;b(1)=7;b(n)=14;b=b';fork=1:n-1[mainElement,index]=max(abs(A(k:n,k)));index=index+k-1;%indexifabs(mainElement)<epsdisp('列元素太小！！');break;elseifindex>ktemp=A(k,:);temp1=b(k);A(k,:)=A(index,:);b(k)=b(index);A(index,:)=temp;b(index)=temp1;endfori=k+1:nm(i,k)=A(i,k)/A(k,k);%A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endendb(n)=b(n)/A(n,n);fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i);endclearx;x=b;disp('AX=b的解x是')end计算结果：在matlab命令框里输出GLZX（10）得：>>GLZX(10)AX=b的解x是ans=1111111111在matlab命令框里输出GLZX（84）得：>>GLZX(84)AX=b的解x是ans=1.00001.00001.00001.00001.00001.0000………1.00001.00001.00001.00001.00001.00001.00001.0000分析：在n较小时，两种方法均能得到正确解，当n较大后，Gauss消去法计算结果严重偏离准确值成为错解，列主元消去法依然能得到正确解。这是因为Gauss消去法中有小主元做除数，在计算过程中的舍入误差会对解产生极大影响，从而导致错误。列主元消去法则避免了这种情况。2.设有方程组，其中，选取不同的初始向量和不同的右端向量，给定迭代误差要求，用Jacobi和Gauss-Seidel迭代法计算，观测得出的迭代向量序列是否收敛。若收敛，记录迭代次数，分析计算结果并得出你的结论。选定初始向量初始向量和不同的右端向量，如取。将的主对角线元素成倍增长若干次，非主对角元素不变，每次用Jacobi法计算，要求迭代误差满足，比较收敛速度，分析现象并得出你的结论。Matlab编程：Jacobi迭代法：functionx=Jacobi(b,x0)A=zeros(20);fori=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=D\(L+U);f=D\b;x=x0;i=0;ep=10^-5;whilei<1000y=x;x=B*x+f;i=i+1;ifnorm(x-y,inf)<epbreak;endendifi<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endendGauss-Seidel迭代法：functionx=G_S(b,x0)A=zeros(20);fori=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=(D-L)\U;f=(D-L)\b;x=x0;i=0;ep=10^-5;whilei<1000y=x;x=B*x+f;i=i+1;ifnorm(x-y,inf)<epbreak;endendifi<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endend计算结果：设误差要求为1×10^-5任选了两组不同的初始向量和不同的右端向量：b=[1,1,1,1,1,1,1,1,2,3,4,1,4,2,3,1,3,1,4,6]'x0=[2,3,1,4,4,5,5,3,1,3,1,1,1,1,1,1,1,1,1,1]'b=[1,4,5,12,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]'x0=[0,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]'在matlab内运行两个程序得到：第一组：>>Jacobi(b,x0)迭代收敛且迭代次数为i=19>>G_S(b,x0)迭代收敛且迭代次数为i=9第二组：>>Jacobi(b,x0)迭代收敛且迭代次数为i=18>>G_S(b,x0)迭代收敛且迭代次数为i=9分析：两组均收敛，明显在均收敛的情况下Gauss-Seidel迭代法要比Jacobi迭代法收敛速度快。这是因为Gauss-Seidel迭代法对已经计算出来的信息进行了充分利用。(b)Matlab编程：functionx=Jacobi2(n)A=zeros(20);x=diag(A);e=diag(eye(20));fori=1:18A(i,i)=3*n;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3*n;A(20,20)=3*n;A(19,20)=-0.5;A(20,19)=-0.5;b=A*e;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=D\(L+U);f=D\b;i=0;ep=10^-6;whilei<1000y=x;x=B*x+f;i=i+1;ifnorm(x-y,inf)<epbreak;endendifi<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endend计算结果：n代表主元增长倍数，选取不同的n得：>>Jacobi2(1)迭代收敛且迭代次数为i=20>>Jacobi2(2)迭代收敛且迭代次数为i=11>>Jacobi2(3)迭代收敛且迭代次数为i=9>>Jacobi2(4)迭代收敛且迭代次数为i=8>>Jacobi2(10)迭代收敛且迭代次数为i=6分析：随着矩阵A主对角元增长倍数的变大，收敛速度变快，但是这种改善效果会逐渐减小。3.用迭代法求方程的全部根，要求误差限为。首先经过简单计算并结合图形得知该方程的三个单根区间为[-3,-1],[-1,0],[0,1],然后利用二分法求解。Matlab编程：functionx=f(a,b)fa=a^3+3*a^2-1;fb=b^3+3*b^2-1;iffa*fb>0disp('区间不是单根区间')return;endep=0.5*10^-8;whileabs(a-b)>epc=(a+b)/2;f1=a^3+3*a^2-1;f2=c^3+3*c^2-1;j=f1*f2;ifj==0x=c;break;elseifj<0b=c;elsea=c;endx=a;endenddisp('该区间内根为')end计算结果：>>f(-3,-1)该区间内根为ans=-2.879385244101286>>f(-1,0)该区间内根为ans=-0.652703646570444>>f(0,1)该区间内根为ans=0.5320888832211494.给定数据如下表：0123456789100.00.791.532.192.713.033.272.893.063.193.29编制程序求三次样条插值函数在插值中点的样条函数插值，并作点集和样条插值函数的图形，满足的边界条件为（1）（2）(1)Matlab编程：clearallx=[012345678910];f1=0.8;f2=0.2;y=[00.791.532.192.713.033.272.893.063.193.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n-2);yn=[];xn=[];fork=2:nhk=x(n)-x(n-1);h=[hhk];endfork=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[ggk];t=[t,tk];u=[uuk];endfork=2:n-2A(k-1,k)=u(k-1);A(k,k-1)=t(k);endg(1)=g(1)-t(1)*f1;g(n-2)=g(n-2)-u(n-2)*f2;m=A\g';m=m';m=[f1mf2];fork=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2;a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2;b2=y(k+1)/h(k)^2;c1=[1-x(k)];c2=[1-x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfork=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[ynyi];xn=[xnxi];endxn=[xnx(n)];yn=[ynpolyval(s(n-1,:),x(n))];plot(xn,yn,'x');title('3次样条函数插值结果')计算结果：在Matlab命令框内输入上述程序得到插值中点的样条函数值Y和样条插值函数的图形：Y=0.398564254606255Y=1.168428726968726Y=1.871470837518841Y=2.478187922956004Y=2.873277470657021Y=3.213702194414573Y=3.084413751685133Y=2.919892798844714Y=3.149765052935493Y=3.222296989411632(2)Matlab编程：clearallx=[012345678910];f1=0;f2=0;y=[00.791.532.192.713.033.272.893.063.193.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n);yn=[];xn=[];fork=2:nhk=x(n)-x(n-1);h=[hhk];endfork=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[ggk];t=[t,tk];u=[uuk];endfork=1:n-2A(k+1,k+2)=u(k);A(k+1,k)=t(k);endA(1,2)=1;A(n,n-1)=1;g0=3*(y(2)-y(1))/h(1)-h(1)*f1/2;gn=3*(y(n)-y(n-1))/h(n-2)+h(n-2)*f2/2;g=[g0,g,gn];m=A\g';m=m';fork=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2;a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2;b2=y(k+1)/h(k)^2;c1=[1-x(k)];c2=[1-x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfork=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[ynyi];xn=[xnxi];endxn=[xnx(n)];yn=[ynpolyval(s(n-1,:),x(n))];plot(xn,yn,'x');title('3次样条函数插值结果')计算结果：在命令框内输入上述程序得到：Y=0.398428148708663Y=1.168465553874013Y=1.871459635795274Y=2.478195902944736Y=2.873256752425371Y=3.213777087353492Y=3.084134898160016Y=2.920933320005332Y=3.145881821815571Y=3.236789392730741再次输入plot(x,y,'x');可以得到点集的图形5.对下列数据作三次多项式拟合，取权系数，给出拟合多项式的系数、平方误差，并作离散数据和拟合多项式的图形。-1.0-0.50.00.51.01.52.0-4.447-0.4520.5510.048-0.4470.5494.552Matlab编程：clearallx=[-1,-0.5,0,0.5,1,1.5,2];y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];a=polyfit(x,y,3)m=polyval(a,x)-y;e=norm(m,2);ep=e^2xi=[x(1):0.02:x(7)];yi=polyval(a,xi);plot(xi,yi,'x');计算结果：在Matlab命令框内输入上述程序得到：a=1.9991-2.9977-0.00000.5491ep=2.1762e-005数组a为拟合多项式的系数，ep为平方误差。拟合多项式的图形：再输入plot(x,y,'x');得到：离散数据图形：6.用复化梯形公式、复化Simpson公式和复化Gauss-Legendre公式计算下列积分的近似值，使绝对误差限为，并将计算结果与精确解作比较以及比较各种算法的计算量。（1），（2）.复化梯形公式：Matlab编程：functions=tixing(f,a,b,t)ep=0.5*10^-7;n=2;s=0;whileabs(s-t)>eph=(b-a)/n;s1=0;fork=1:n-1x=a+h*k;s1=s1+feval(f,x);endn=n+1;s=h*(feval(f,a)+feval(f,b)+2*s1)/2;endn=n-1end计算结果：在matlab命令框内输入tixing(@(x)1/x,1,2,log(2))得到积分（1）的计算量（用分段数n刻画）以及计算结果：n=1119ans=0.693147230473649在matlab命令框内输入tixing(@(x)1/(1+x^2),0,1,pi/4)得到积分（2）的计算量（用分段数n刻画）以及计算结果：n=913ans=0.785398113411584复化Simpson公式：Matlab编程：functions=simpson(f,a,b,t)ep=0.5*10^-7;N=2;s=0;whileabs(s-t)>eph=(b-a)/(2*N);s1=0;s2=0;fork=1:Nx=a+h*(2*k-1);s1=s1+feval(f,x);endfork=1:(N-1)x=a+h*2*k;s2=s2+feval(f,x);endN=N+1;s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;endN=N-1end计算结果：在matlab命令框内输入simpson(@(x)1/x,1,2,log(2))得到积分（1）的计算量（用分段数N刻画）以及计算结果：（值得注意的是，simpson公式的分段数没有刻画出每段中点值的计算量）N=15ans=0.693147219033552在matlab命令框内输入simpson(@(x)1/(1+x^2),0,1,pi/4)得到积分（2）的计算量（用分段数N刻画）以及计算结果：N=4ans=0.785398125614677复化Gauss-Legendre公式:Matlab编程：functiony=getLegendre(n,x)%返回Legendre多项式ifn==1y=x;elseifn==2y=1.5*x^2-0.5;elsey=(2-1/n)*x*getLegendre(n-1,x)-(n-1)/n*getLegendre(n-2,x);endendfunctionquad=GL(f,a,b,t)%复化Gauss-Legendre公式计算ep=0.5*10^-7;n=1;quad=0;symsx;whileabs(quad-t)>epg=getLegendre(n,x);temp=sym2poly(g);X=roots(temp);g=diff(g);g=g*g;N=length(X);A=zeros(N,1);whileN>0A(N)=2/(1-X(N)*X(N))/subs(g,'x',X(N));N=N-1;endN=length(X);T=zeros(N,1);M=zeros(1,N);T=((a+b)/2)+((b-a)/2)*X;fori=1:NM(i)=feval(f,T(i));endquad=((b-a)/2)*(M*A);n=n+1;endn=n-1end计算结果：在matlab命令框内输入GL(@(x)1/x,1,2,log(2))得到积分（1）的计算量（用Legendre多项式的次数n刻画）以及计算结果：n=5ans=0.693147157853040在matlab命令框内输入GL(@(x)1/(1+x^2),0,1,pi/4)得到积分（2）的计算量（用Legendre多项式的次数n刻画）以及计算结果：n=5ans=0.785398159971188分析：通过比较发现，精确性：复化Gauss-Legendre公式>复化Simpson公式>复化梯形公式，计算量：复化梯形公式>复化Simpson公式>复化Gauss-Legendre公式。7.用4阶Runge-Kutta法求解微分方程：分别用和计算，并比较两个近似值与精确解Matlab编程：functionR=rk4(f,a,b,ya,h)N=(b-a)/h;T=zeros(1,N+1);Y=zeros(1,N+1);y=zeros(1,N+1);T=a:h:b;Y(1)=ya;y(1)=ya;g=@(x)x^2*exp(-x)+cos(2*x);fork=1:Nk1=feval(f,T(k),Y(k));k2=feval(f,T(k)+h/2,Y(k)+h*k1/2);k3=feval(f,T(k)+h/2,Y(k)+h*k2/2);k4=feval(f,T(k)+h,Y(k)+h*k3);Y(k+1)=Y(k)+h*(k1+2*k2+2*k3+k4)/6;y(k+1)=feval(g,T(k));endR=[T'Y'y'];end计算结果：>>rk4(@(x,y)-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x),0,2,1,0.1)ans=01.0000000000000001.0000000000000000.1000000000000000.9891149457346890.9891149520216010.2000000000000000.9538102071638440.9538102241260040.3000000000000000.8920092348744880.8920092547710330.4000000000000000.8039579098467340.8039579167128680.5000000000000000.6919349976895810.6919349707962980.6000000000000000.5599300266012580.5599299434705230.7000000000000000.4132941027862110.4132939417580310.8000000000000000.2583712723870450.2583710147337330.9000000000000000.1021196981119560.1021193296967981.000000000000000-0.048266907857711-0.0482673953757001.100000000000000-0.185726497573598-0.1857271059806701.200000000000000-0.303673326189410-0.3036740503876741.300000000000000-0.396309195058830-0.3963100231414661.400000000000000-0.458891377646222-0.4588922913431091.500000000000000-0.487948660815802-0.4879496362664781.600000000000000-0.481438680925772-0.4814396897284351.700000000000000-0.438841797588724-0.4388428080670581.800000000000000-0.361189039871900-0.3611900184962071.900000000000000-0.251024083624937-0.2510249965207042.000000000000000-0.112301673441313-0.112302487917161>>rk4(@(x,y)-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x),0,2,1,0.05)ans=01.0000000000000001.0000000000000000.0500000000000000.9973822387440070.9973822388392770.1000000000000000.9891149517069640.9891149520216010.1500000000000000.9747024180466210.9747024185951700.2000000000000000.9538102234221290.9538102241260040.2500000000000000.9262576101293000.9262576108323360.3000000000000000.8920092542879130.8920092547710330.3500000000000000.8511664782800030.8511664782750310.4000000000000000.8039579175102620.8039579167128680.4500000000000000.7507296708922890.7507296689740730.5000000000000000.6919349741765920.6919349707962980.5500000000000000.6281234442518160.6281234390656750.6000000000000000.5599299507993260.5599299434705230.6500000000000000.4880631790988550.4880631693061170.7000000000000000.4132939543127370.4132939417580310.7500000000000000.3364434031691730.3364433875845240.8000000000000000.2583710335803170.2583710147337330.8500000000000000.1799628163369100.1799627940374300.9000000000000000.1021193555949370.1021193296967980.9500000000000000.0257442363982430.0257442068041841.0000