大连理工大学矩阵与数值分析MATLAB上机实验

矩阵与数值分析之MATLAB试验学生姓名：陈贵鹏学号：21606075专业班级：结构工程建研1602班所在单位：智能结构研究所2017年6月4日2016级工科硕士研究生《矩阵与数值分析》课程数值实验题目教学班号：建研1602班任课教师：董波学生姓名：陈贵鹏学号：21606075院系：建设工程学部结构工程系一、设21141NNjSj，分别编制从小到大和从大到小的顺序程序计算100100001000000,,SSS，并指出有效位数。从小到大求和程序：clc;clear;formatlongS=0;N=100;%此处N分别取100,10000,1000000进行计算forj=1:NS=S+1/(4*j^2-1);endS从小到大求和程序计算结果：N从小到大求和程序得到的SN真实值𝑆𝑁=𝑁2𝑁+1计算值有效位数1000.4975124378109450.49751243781094515100000.4999750012499370.4999750012499371510000000.4999997500001340.49999975000012513从大到小求和程序：clc;clear;formatlongS=0;N=1000000;%此处N分别取100,10000,1000000进行计算forj=N:(-1):1S=S+1/(4*j^2-1);endS从大到小求和程序计算结果N从大到小求和程序得到的SN真实值𝑆𝑁=𝑁2𝑁+1有效位数1000.4975124378109450.49751243781094515100000.4999750012499370.4999750012499371510000000.4999997500001250.49999975000012515分析：由2111111=()412212121NNNjjNSjjjN可以获得真实值。当N很大时，可以发现从小到大计算与从大到小的计算结果出现了差异。这是因为从大到小求和时，变量在分母上，整个算式是由小到大进行求和的，这样求和避免了“大数吃小数”现象，因而计算结果更优，具有更高的有效位数。二、解线性方程组1．分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组12343100113100,0121000130xxxx迭代法计算停止的条件为：6)()1(3110maxkjkjjxx．Jacobi迭代程序:clc;clear;formatlongA=[-3,1,0,0;1,-3,1,0;0,1,-2,1;0,0,1,-3];b=[-1;0;0;0];D=diag(diag(A));%求对角矩阵DL=-tril(A,-1);%求下三角矩阵LU=-triu(A,1);%求上三角矩阵UB=D\(L+U);f=D\b;x0=[0;0;0;0];%取零向量作为初始向量x=B*x0+f;whilenorm(x-x0,'inf')=1e-6x0=x;x=B*x0+f;endxJacobi迭代程序运算结果：Gauss-Seidel迭代程序:clc;clear;formatlongA=[-3,1,0,0;1,-3,1,0;0,1,-2,1;0,0,1,-3];b=[-1;0;0;0];D=diag(diag(A));%求A的对角矩阵DL=-tril(A,-1);%求下三角矩阵LU=-triu(A,1);%求上三角矩阵UG=(D-L)\U;f=(D-L)\b;x0=[0;0;0;0];%取零向量作为初始向量x=G*x0+f;whilenorm(x-x0,'inf')=1e-6x0=x;x=G*x0+f;endxGauss-Seidel迭代程序运算结果:分析由计算结果可见，Jacobi迭代法和Gauss-Seidel求解的结果非常接近。为此我们用x的计算值与真实值之差的2-范数作为误差来衡量两种方法的优劣，在命令行输入norm(x-A\b,2)，Jacobi迭代法的误差结果为1.172109617472753e-06，Gauss-Seidel迭代的误差结果为4.812414751512495e-07，这表明Gauss-Seidel迭代法具有更高的计算精度。这是因为Jacobi迭代法对已计算出来的信息未加以充分利用，一般来说，后面的计算值比前面的计算值要准确些，而Gauss-Seidel迭代考虑了这一点，因此具有更高的精度。2.用Gauss列主元消去法、QR方法求解如下方程组：123424311282006,504023900516xxxxGauss列主元消去法程序:clc;clear;formatlonga=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5];%系数矩阵b=[12;6;23;16];[n,m]=size(a);nb=length(b);det=1;fork=1:n-1amax=0;fori=k:nifabs(a(i,k))amaxamax=abs(a(i,k));r=i;endendifamax1e-10return;endifrkforj=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfori=k+1:nm=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);fork=n:-1:1forj=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endx=x'Gauss列主元消去法程序运算结果:QR方法程序：clc;clear;formatlongA=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5];%系数矩阵B=[12;6;23;16];n=length(B);Q=eye(4);H=A;fori=1:n-1b=H(:,1);c=norm(b,2);b(1)=b(1)-c;b1=transpose(b);Q1=eye(4-i+1)-2/(b1*b)*b*b1;H1=Q1*H;Q1=blkdiag(eye(i-1),Q1);Q=Q*transpose(Q1);ifin-1H1(1,:)=[];H1(:,1)=[];H=H1;endendR=transpose(Q)*A;x=R\(Q\B)QR方法程序运算结果：分析从计算结果来看，如果不计舍入误差，Gauss列主元消去法与QR方法得到的结果几乎相同。但是从编程的复杂程度来看，个人感觉Gauss列主元消去法难度更大，QR方法的语句相对简洁。三、非线性方程的迭代解法1、用Newton迭代法、割线法求方程22cos60xxfxex的根，计算停止的条件为：6110kkxxNewton迭代法程序clc;clear;formatlongx0=1;%设初始值f=inline('exp(x)+2^(-x)+2*cos(x)-6');diff=inline('exp(x)-2^(-x)*log(2)-2*sin(x)');x=x0-feval(f,x0)/feval(diff,x0);whileabs(x-x0)=1e-6x0=x;x=x0-(exp(x0)+2^(-x0)+2*cos(x0)-6)/(exp(x0)-2^(-x0)*log(2)-2*sin(x0));endxNewton迭代法程序运算结果割线法程序clc;clear;formatlongxm=1;%设初始值xn=0.5;%设初始值f=inline('exp(x)+2^(-x)+2*cos(x)-6');whileabs(xn-xm)=1e-6x=xn-(feval(f,xn)*(xn-xm))/(feval(f,xn)-feval(f,xm));xm=xn;xn=x;endx割线法程序运算结果分析在MATLAB命令窗口输入x=solve('exp(x)+2^(-x)+2*cos(x)-6=0')，可以得到方程的根为x=1.8293836019338488171362129468142。对比Newton法和割线法求得的结果，可以发现Newton法求得的结果更精确。这是因为割线法采用了导数的近似式'11()()()kkkkkfxfxfxxx，即用割线代替切线，这引入了误差。2、生成一个列和为1的100阶随机矩阵，编写幂法程序求其最大特征值。幂法程序clc;clear;formatlongA=rand(100,100);%生成100阶随机矩阵fori=1:100%使随机矩阵；列和为1A(:,i)=A(:,i)/sum(A(:,i));endN=100;%最大迭代次数ep=1e-6;%精度要求n=length(A);u=ones(n,1);k=0;m1=0;whilek=Nv=A*u;m=max(abs(v));u=v/m;ifabs(m-m1)epbreak;endm1=m;k=k+1;endm幂法程序运算结果分析如果不计舍入误差，待求列和为1的n阶随机矩阵的最大特征值恒等于1，实际计算结果可能是1，或者1.000000000000001，1.000000000000002等。四、数值积分分别用复化梯形公式和三点Gauss型公式计算积分821dxx要求误差不超过510。复化梯形公式程序clc;clear;formatlongsymstm=int(1/t,2,8);%真实值a=2;b=8;n=300;h=(b-a)/n;sum=0;f=inline('1/x');fori=1:n-1sum=sum+f(a+i.*h);endT=h/2*(f(a)+2*sum+f(b))复化梯形公式程序计算结果分析本程序将区间[2,8]进行了300等分，在命令行输入abs(T-m)，可得计算值与精确值的误差为7.812416996655358e-06，误差未超过10-5，满足题目要求。三点Gauss型公式程序clc;clear;fun=inline('1./x');%被积函数a=2;b=8;tol=1e-5;%精度n=2;%n次代数精度，n=(3-1)*2+1symsxp=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^n*factorial(n));tk=roots(p);%求积节点Ak=zeros(n+1,1);fori=1:n+1xkt=tk;xkt(i)=[];pn=poly(xkt);fp=@(x)polyval(pn,x)/polyval(pn,tk(i));Ak(i)=quadl(fp,-1,1,tol);%求积系数endxk=(b-a)/2*tk+(b+a)/2;fun=fcnchk(fun,'vectorize');fx=fun(xk)*(b-a)/2;ql=sum(Ak.*fx)三点Gauss型公式程序计算结果分析Gauss型求积公式未限制等距节点，是具有最高代数精度的求积公式。本例取三点，因此n=2，代数精度为2n+1=5.五、插值与逼近1．给定1,1上的函数2511xxf，构造等距节点分别取5n，8n，10n的Lagrange插值多项式。画出原函数的图像以及插值函数的图像,观察当n增大时的逼近效果.求Lagrange插值多项式程序clc;clear;N=5;%这个地方分别取N=5,8,10，即可求得在这三种情况下的Lagrange插值多项式x=-1:2/N:1;y=1./(1+5*x.^2);n=length(x);x0=sym('x');f=0;fork=1:nl=1;fori=1:nifi~=kl=l*(