矩阵与数值分析上机作业学校：大连理工大学学院：班级：姓名：学号：授课老师：注：编程语言Matlab程序：Norm.m函数function s=Norm(x,m)%求向量x的范数%m取1,2,inf分别表示1,2，无穷范数n=length(x);s=0;switch mcase 1 %1-范数for i=1:ns=s+abs(x(i));endcase 2 %2-范数for i=1:ns=s+x(i)^2;ends=sqrt(s);case inf %无穷-范数s=max(abs(x));end计算向量x，y的范数Test1.mclear all;clc;n1=10;n2=100;n3=1000;x1=1./[1:n1]';x2=1./[1:n2]';x3=1./[1:n3]'; y1=[1:n1]';y2=[1:n2]';y3=[1:n3]';disp('n=10时');disp('x的1-范数:');disp(Norm(x1,1));disp('x的2-范数:');disp(Norm(x1,2));disp('x的无穷-范数:');disp(Norm(x1,inf)); disp('y的1-范数:');disp(Norm(y1,1));disp('y的2-范数:');disp(Norm(y1,2));disp('y的无穷-范数:');disp(Norm(y1,inf)); disp('n=100时');disp('x的1-范数:');disp(Norm(x2,1));disp('x的2-范数:');disp(Norm(x2,2));disp('x的无穷-范数:');disp(Norm(x2,inf));disp('y的1-范数:');disp(Norm(y2,1));disp('y的2-范数:');disp(Norm(y2,2));disp('y的无穷-范数:');disp(Norm(y2,inf));disp('n=1000时');disp('x的1-范数:');disp(Norm(x3,1));disp('x的2-范数:');disp(Norm(x3,2));disp('x的无穷-范数:');disp(Norm(x3,inf));disp('y的1-范数:');disp(Norm(y3,1));disp('y的2-范数:');disp(Norm(y3,2));disp('y的无穷-范数:');disp(Norm(y3,inf));运行结果：n=10时x的1-范数:2.9290；x的2-范数:1.2449； x的无穷-范数:1y的1-范数:55； y的2-范数:19.6214； y的无穷-范数:10n=100时x的1-范数:5.1874；x的2-范数: 1.2787； x的无穷-范数:1y的1-范数:5050； y的2-范数:581.6786； y的无穷-范数:100n=1000时x的1-范数:7.4855； x的2-范数:1.2822； x的无穷-范数:1y的1-范数: 500500； y的2-范数:1.8271e+004；y的无穷-范数:1000程序Test2.mclear all;clc;n=100;%区间h=2*10^(-15)/n;%步长x=-10^(-15):h:10^(-15);%第一种原函数f1=zeros(1,n+1);for k=1:n+1if x(k)~=0f1(k)=log(1+x(k))/x(k);elsef1(k)=1;endendsubplot(2,1,1);plot(x,f1,'-r');axis([-10^(-15),10^(-15),-1,2]); legend('原图');%第二种算法f2=zeros(1,n+1);for k=1:n+1d=1+x(k);if(d~=1)f2(k)=log(d)/(d-1);elsef2(k)=1;endendsubplot(2,1,2);plot(x,f2,'-r');axis([-10^(-15),10^(-15),-1,2]); legend('第二种算法');运行结果：显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当]10,10[1515--∈x 时，x d +=1，计算机进行舍入造成d 恒等于1，结果函数值恒为1。

程序：秦九韶算法:QinJS.mfunction y=QinJS(a,x) %y 输出函数值%a 多项式系数，由高次到零次 %x 给定点n=length(a); s=a(1); for i=2:n s=s*x+a(i); end y=s;计算p(x):test3.mclear all ; clc;x=1.6:0.2:2.4;%x=2的邻域disp('x=2的邻域：');xa=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; p=zeros(1,5); for i=1:5p(i)=QinJS(a,x(i)); enddisp('相应多项式p 值：');p xk=1.95:0.01:20.5; nk=length(xk); pk=zeros(1,nk); k=1; for k=1:nkpk(k)=QinJS(a,xk(k)); endplot(xk,pk,'-r');xlabel('x');ylabel('p(x)');运行结果：x=2的邻域：x =1.6000 1.80002.0000 2.2000 2.4000相应多项式p值：p =1.0e-003 *-0.2621 -0.0005 0 0.0005 0.2621 p(x)在 x[1.95,20.5]上的图像程序：LU分解，LUDecom.mfunction [L,U]=LUDecom(A)%不带列主元的LU分解N = size(A);n = N(1);L=eye(n);U=zeros(n);for i=1:nU(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1); endfor i=2:nfor j=i:nz=0;for k=1:i-1z=z+L(i,k)*U(k,j);endU(i,j)=A(i,j)-z;endfor j=i+1:nz=0;for k=1:i-1z=z+L(j,k)*U(k,i);endL(j,i)=(A(j,i)-z)/U(i,i);endendPLU分解，PLUDecom.mfunction [P,L,U] =PLUDecom(A)%带列主元的LU分解[m,m]=size(A);U=A;P=eye(m);L=eye(m);for i=1:mfor j=i:mt(j)=U(j,i);for k=1:i-1t(j)=t(j)-U(j,k)*U(k,i);endenda=i;b=abs(t(i));for j=i+1:mif b<abs(t(j))b=abs(t(j));a=j;endendif a~=ifor j=1:mc=U(i,j);U(i,j)=U(a,j);U(a,j)=c;endfor j=1:mc=P(i,j);P(i,j)=P(a,j);P(a,j)=c;endc=t(a);t(a)=t(i);t(i)=c;endU(i,i)=t(i);for j=i+1:mU(j,i)=t(j)/t(i);endfor j=i+1:mfor k=1:i-1U(i,j)=U(i,j)-U(i,k)*U(k,j);endendendL=tril(U,-1)+eye(m);U=triu(U,0);（1）（2）程序：Test4.mclear all;clc;for n=5:30x=zeros(n,1);A=-ones(n);A(:,n)=ones(n,1);for i=1:nA(i,i)=1;for j=(i+1):(n-1)A(i,j)=0;endx(i)=1/i;enddisp('当n=');disp(n);disp('方程精确解:');xb=A*x; %系数bdisp('利用LU分解方程组的解:');[L,U]=LUDecom(A); %LU分解xLU=U\(L\b)disp('利用PLU分解方程组的解:');[P,L,U] =PLUDecom(A); %PLU分解xPLU=U\(L\(P\b))%求解A的逆矩阵disp('A的准确逆矩阵:');InvA=inv(A)InvAL=zeros(n); %利用LU分解求A的逆矩阵 I=eye(n);for i=1:nInvAL(:,i)=U\(L\I(:,i));enddisp('利用LU分解的A的逆矩阵:');InvALEnd运行结果：（1）只列出n=5,6,7的结果当n= 5方程精确解:x =1.00000.50000.33330.25000.2000利用LU分解方程组的解:xLU =1.00000.50000.33330.25000.2000利用PLU分解方程组的解:xPLU =1.00000.50000.33330.25000.2000当n=6方程精确解:x =1.00000.50000.33330.25000.20000.1667利用LU分解方程组的解: xLU =1.00000.50000.33330.25000.20000.1667利用PLU分解方程组的解: xPLU =1.00000.50000.33330.25000.20000.1667当n= 7方程精确解:x =1.00000.50000.33330.25000.20000.16670.1429利用LU分解方程组的解: xLU =1.00000.50000.33330.25000.20000.16670.1429利用PLU分解方程组的解:xPLU =1.00000.50000.33330.25000.20000.16670.1429（2）只列出n=5,6,7时A的逆矩阵的结果当n= 5A的准确逆矩阵:InvA =0.5000 -0.2500 -0.1250 -0.0625 -0.06250 0.5000 -0.2500 -0.1250 -0.12500 0 0.5000 -0.2500 -0.25000 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0625利用LU分解的A的逆矩阵:InvAL =0.5000 -0.2500 -0.1250 -0.0625 -0.06250 0.5000 -0.2500 -0.1250 -0.12500 0 0.5000 -0.2500 -0.25000 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0625当n= 6A的准确逆矩阵:InvA =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313 利用LU分解的A的逆矩阵:InvAL =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.03130 0.5000 -0.2500 -0.1250 -0.0625 -0.06250 0 0.5000 -0.2500 -0.1250 -0.12500 0 0 0.5000 -0.2500 -0.25000 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0313当n= 7A的准确逆矩阵:InvA =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156利用LU分解的A的逆矩阵:InvAL =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156程序：Cholesky分解：Cholesky.mfunction L=Cholesky(A)N = size(A);n = N(1);L=zeros(n);L(1,1)=sqrt(A(1,1));for i=2:nL(i,1)=A(i,1)/L(1,1);endfor j=2:ns1=0;for k=1:j-1s1=s1+L(j,k)^2;endL(j,j)=sqrt(A(j,j)-s1);for i=j+1:ns2=0;for k=1:j-1s2=s2+L(i,k)*L(j,k);endL(i,j)=(A(i,j)-s2)/L(j,j);endend计算Ax=b；Test5.mclear all;clc;for n=10:20A=zeros(n,n);b=zeros(n,1);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endb(i,1)=i;enddisp('n=');disp(n);disp('方程组原始解');x0=A\bdisp('利用Cholesky分解的方程组的解'); L=Cholesky(A)x=L'\(L\b)end运行结果：只列出了n=10,11的结果n=10方程组原始解x0 =1.0e+008 *-0.00000.0010-0.02330.2330-1.21083.5947-6.32336.5114-3.62330.8407利用Cholesky分解的方程组的解x =1.0e+008 *-0.00000.0010-0.02330.2330-1.21053.5939-6.32196.5100-3.62250.8405n=11方程组原始解x0 =1.0e+009 *0.0000-0.00020.0046-0.05670.3687-1.40393.2863-4.78694.2260-2.06850.4305利用Cholesky分解的方程组的解x =1.0e+009 *0.0000-0.00020.0046-0.05630.3668-1.39723.2716-4.76694.2094-2.06080.4290程序：（1）House.mfunction u=House(x)n=length(x);e1=eye(n,1);w=x-norm(x,2)*e1;u=w/norm(w,2);(2)Hou_A.mfunction HA=Hou_A(A)a1=A(:,1);n=length(a1);e1=eye(n,1);w=a1-norm(a1,2)*e1;u=w/norm(w,2);H=eye(n)-2*u*u'HA=H*A;(3)test6.mclear all;clc;A=[1 2 3 4;-1 2 sqrt(2) sqrt(3);-2 2 exp(1) pi;-sqrt(10) 2 -3 7;0 2 7 5/2];HA=Hou_A(A)运行结果：H =0.2500 -0.2500 -0.5000 -0.7906 0 -0.2500 0.9167 -0.1667 -0.2635 0 -0.5000 -0.1667 0.6667 -0.5270 0 -0.7906 -0.2635 -0.5270 0.1667 00 0 0 0 1.0000 HA =4.0000 -2.5811 1.4090 -6.53780.0000 0.4730 0.8839 -1.78050.0000 -1.0541 1.6576 -3.88360.0000 -2.8289 -4.6770 -4.10780 2.0000 7.0000 2.5000程序：Jacobi迭代：Jaccobi.mfunction [x,n]=Jaccobi(A,b,x0)%--·方程组系数阵A%--·方程组右端顶b%-- 初始值x0%--求解要求精确度eps%--迭代步数控制M%--·返回求得的解x%--·返回迭代步数nM=1000;eps=1.0e-5;D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵J=D\(L+U);f=D\b;x=J*x0+fn=1; %迭代次数err=norm(x-x0,inf)while(err>=eps)x0=x;x=J*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M)disp('Warning: 迭代次数太多，可能不收敛?');return;endendGauss_Seidel迭代：Gauss_Seidel.mfunction [x,n]=Gauss_Seidel(A,b,x0)%--Gauss-Seidel迭代法解线性方程组%--方程组系数阵 A%--方程组右端项 b%--初始值 x0%--求解要求的精确度 eps%--迭代步数控制 M%--返回求得的解 x%--返回迭代步数 neps=1.0e-5;M=10000;D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+fn=1; %迭代次数err=norm(x-x0,inf)while(err>=eps)x0=x;x=G*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endend解方程组，test7.mclear all;clc;A=[5 -1 -3;-1 2 4;-3 4 15];b=[-2;1;10];disp('精确解');x=A\bdisp('迭代初始值');x0=[0;0;0]disp('Jacobi迭代过程：');[xj,nj]=Jaccobi(A,b,x0);disp('Jacobi最终迭代结果：');xjdisp('迭代次数');njdisp('Gauss-Seidel迭代过程：'); [xg,ng]=Gauss_Seidel(A,b,x0);disp('Gauss-Seidel最终迭代结果：'); xgdisp('迭代次数');ng运行结果：精确解x =-0.0820-1.80331.1311迭代初始值x0 =Jacobi迭代过程：x =-0.40000.50000.6667err =0.6667x =0.1000-1.03330.4533err =1.5333...x =-0.0820-1.80331.1311err =9.6603e-006Jacobi最终迭代结果：xj =-0.0820-1.80331.1311迭代次数nj =281Gauss-Seidel迭代过程：x =-0.40000.30000.5067err =0.5067x =-0.0360-0.53130.8012err =0.8313x =-0.0256-1.11510.9589err =0.5838...x =-0.0820-1.80331.1311err =9.4021e-006Gauss-Seidel最终迭代结果：xg =-0.0820-1.80331.1311迭代次数ng =20程序：Newton迭代法：Newtoniter.mfunction [x,iter,fvalue]=Newtoniter(f,df,x0,eps,maxiter) %牛顿法 x得到的近似解%iter迭代次数%fvalue函数在x处的值%f，df被求的非线性方程及导函数%x0初始值%eps 允许误差限%maxiter 最大迭代次数fvalue=subs(f,x0);dfvalue=subs(df,x0);for iter=1:maxiterx=x0-fvalue/dfvalueerr=abs(x-x0)x0=x;fvalue=subs(f,x0)dfvalue=subs(df,x0);if(err<eps)||(fvalue==0),break,endend弦截法：secant.mfunction [x,iter,fvalue]=secant(f,x0,x1,eps,maxiter)%弦截法 x得到的近似解%iter迭代次数%fvalue函数在x处的值%f被求的非线性方程%x0，x1初始值%eps 允许误差限%maxiter 最大迭代次数fvalue0=subs(f,x0);fvalue=subs(f,x1);for iter=1:maxiterx=x1-fvalue*(x1-x0)/(fvalue-fvalue0)err=abs(x-x1)x0=x1;x1=x;fvalue0=subs(f,x0);fvalue=subs(f,x1)if(err<eps)||(fvalue==0),break,endend求方程的实根：test8.mclear all；clc;syms xf=x.^3+2*x.^2+10*x-100;df=diff(f,x,1);eps=10e-6;maxiter=100;disp('Newton迭代初始值');xn1_0=0disp('Newton迭代结果');[xn1,iter_n1,fxn1]=Newtoniter(f,df,xn1_0,eps,maxiter) disp('Newton迭代初始值');xn2_0=5disp('Newton迭代结果');[xn2,iter_n2,fxn2]=Newtoniter(f,df,xn2_0,eps,maxiter) disp('弦截法初始值');xk1_0=0xk1_1=1disp('弦截法迭代结果');[xk1,iter_k1,fxk1]=secant(f,xk1_0,xk1_1,eps,maxiter) disp('弦截法初始值');xk2_0=5xk2_1=6disp('弦截法迭代结果');[xk2,iter_k2,fxk2]=secant(f,xk2_0,xk2_1,eps,maxiter)运行结果：Newton法结果:取两个不同初值0，5程序：二分法：resecm.mfunction [x,iter]=resecm(f,a,b,eps) %二分法 x 近似解%iter 迭代次数%f 求解的方程%a，b 求解区间%eps 允许误差限fa=subs(f,a);fb=subs(f,b);iter=0;if(fa==0)x=a;returnendif(fb==0)x=b;returnendwhile(abs(a-b)>=eps)mf=subs(f,(a+b)/2);if(mf==0)x=mf;n=n+1;returnendif(mf*fa<0)b=(a+b)/2;elsea=(a+b)/2;enditer=iter+1;endx=(a+b)/2;iter=iter+1;求方程的实根：test9.mclear all;clc;syms xf=exp(x).*cos(x)+2;a=0;a1=pi;a2=2*pi;a3=3*pi;b=4*pi;eps=10e-6;[x1,iter1]=resecm(f,a,a1,eps)[x2,iter2]=resecm(f,a1,a2,eps)[x3,iter3]=resecm(f,a2,a3,eps)[x4,iter4]=resecm(f,a3,b,eps)运行结果：[0,pi]区间的根x1 =1.8807；迭代次数iter1 = 20[pi,2*pi]区间的根x2 =4.6941；迭代次数iter2 =20[2*pi,3*pi]区间的根x3 =7.8548；迭代次数iter3 =20[3*pi,4*pi]区间的根x4 =10.9955；迭代次数iter4 =20程序：Newton插值：Newtominter.mfunction f=Newtominter(x,y,x0)%牛顿插值 x插值节点%y为对应的函数值%函数返回Newton插值多项式在x_0点的值fsyms t;if(length(x) == length(y))n = length(x);c(1:n) = 0.0;elsedisp('x和y的维数不相等！');return;endf = y(1);y1 = 0;l = 1;for(i=1:n-1)for(j=i+1:n)y1(j) = (y(j)-y(i))/(x(j)-x(i));endc(i) = y1(i+1);l = l*(t-x(i));f = f + c(i)*l;simplify(f);y = y1;if(i==n-1)if(nargin == 3) %如果3个参数则给出插值点的插值结果 f = subs(f,'t',x0);else%如果2个参数则直接给出插值多项式f = collect(f); %将插值多项式展开f = vpa(f, 6);endendend用等距节点做f(x)的Newton插值：test10.mn1=5;n2=10;n3=15;x0=[0:0.01:1];y0=sin(pi.*x0);x1=linspace(0,1,n1);%等距节点,节点数5y1=sin(pi.*x1);f01=Newtominter(x1,y1,x0);x2=linspace(0,1,n2);%等距节点,节点数10y2=sin(pi.*x2);f02 = Newtominter(x2,y2,x0);x3=linspace(0,1,n3);%等距节点,节点数15y3=sin(pi.*x3);f03= Newtominter(x3,y3,x0);plot(x0,y0,'-r')%原图hold onplot(x0,f01,'-g')%5个节点plot(x0,f02,'-k')%10个节点plot(x0,f03,'-b')%15个节点legend('原图','5个节点Newton插值多项式','10个节点Newton插值多项式','15个节点Newton插值多项式')运行结果：取不同的节点做牛顿插值。