大连理工大学矩阵与数值分析上机作业

1 矩阵与数值分析上机作业 学校： 大连理工大学 学院： 班级： 姓名： 2 学号： 授课老师： 注：编程语言注：编程语言 Matlab 程序：程序： Norm.mNorm.m 函数函数 function s=Norm(x,m) %求向量 x 的范数 %m 取 1,2,inf 分别 表示 1,2，无穷范数 n=length(x); 3 s=0; switch m case 1 %1-范数 for i=1:n s=s+abs(x(i); end case 2 %2-范数 for i=1:n s=s+x(i)2; end s=sqrt(s); case inf %无穷-范数 s=max(abs(x); end 计算向量计算向量 x x，y y 的范数的范数 Test1.mTest1.m clear 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); 4 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=10n=10 时时 x 的 1-范数:2.9290；x 的 2-范数:1.2449； x 的无穷-范数:1 y 的 1-范数:55； y 的 2-范数:19.6214； y 的无穷-范数:10 n=100n=100 时时 x 的 1-范数:5.1874；x 的 2-范数: 1.2787； x 的无穷-范数:1 y 的 1-范数:5050； y 的 2-范数:581.6786； y 的无穷-范数:100 5 n=1000n=1000 时时 x 的 1-范数:7.4855； x 的 2-范数:1.2822； x 的无穷-范数:1 y 的 1-范数: 500500； y 的 2-范数:1.8271e+004；y 的无穷-范数:1000 程序程序 Test2.mTest2.m clear all; clc; n=100;%区间 h=2*10(-15)/n;%步长 x=-10(-15):h:10(-15); %第一种原函数 f1=zeros(1,n+1); for k=1:n+1 if x(k)=0 f1(k)=log(1+x(k)/x(k); else f1(k)=1; 6 end end subplot(2,1,1); plot(x,f1,-r); axis(-10(-15),10(-15),-1,2); legend(原图); %第二种算法 f2=zeros(1,n+1); for k=1:n+1 d=1+x(k); if(d=1) f2(k)=log(d)/(d-1); else f2(k)=1; end end subplot(2,1,2); plot(x,f2,-r); axis(-10(-15),10(-15),-1,2); legend(第二种算法); 运行结果：运行结果： 7 显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当时，10,10 1515 x ，计算机进行舍入造成 恒等于 1，结果函数值恒为 1。xd1d 程序：程序： 秦九韶算法秦九韶算法:QinJS.m:QinJS.m function y=QinJS(a,x) %y 输出函数值 %a 多项式系数，由高次到零次 %x 给定点 n=length(a); s=a(1); for i=2:n 8 s=s*x+a(i); end y=s; 计算计算 p(x):test3.mp(x):test3.m clear all; clc; x=1.6:0.2:2.4;%x=2 的邻域 disp(x=2 的邻域：);x a=1 -18 144 -672 2016 -4032 5376 -4608 2304 -512; p=zeros(1,5); for i=1:5 p(i)=QinJS(a,x(i); end disp(相应多项式 p 值：);p xk=1.95:0.01:20.5; nk=length(xk); pk=zeros(1,nk); k=1; for k=1:nk pk(k)=QinJS(a,xk(k); end plot(xk,pk,-r); xlabel(x);ylabel(p(x); 9 运行结果：运行结果： x=2 的邻域： x = 1.6000 1.8000 2.0000 2.2000 2.4000 相应多项式 p 值： p = 1.0e-003 * -0.2621 -0.0005 0 0.0005 0.2621 p(x)在在1.95,20.5上的图像上的图像x 10 程序：程序： LULU 分解，分解，LUDecom.mLUDecom.m function L,U=LUDecom(A) %不带列主元的 LU 分解 N = size(A); n = N(1); L=eye(n);U=zeros(n); for i=1:n U(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1); end for i=2:n for j=i:n z=0; for k=1:i-1 z=z+L(i,k)*U(k,j); end 11 U(i,j)=A(i,j)-z; end for j=i+1:n z=0; for k=1:i-1 z=z+L(j,k)*U(k,i); end L(j,i)=(A(j,i)-z)/U(i,i); end end PLUPLU 分解，分解，PLUDecom.mPLUDecom.m function P,L,U =PLUDecom(A) %带列主元的 LU 分解 m,m=size(A); U=A; P=eye(m); L=eye(m); for i=1:m for j=i:m t(j)=U(j,i); for k=1:i-1 t(j)=t(j)-U(j,k)*U(k,i); end end 12 a=i;b=abs(t(i); for j=i+1:m if b=eps) x0=x; x=J*x0+f n=n+1; err=norm(x-x0,inf) if(n=M) disp(Warning: 迭代次数太多，可能不收敛?); return; end end Gauss_SeidelGauss_Seidel 迭代：迭代：Gauss_Seidel.mGauss_Seidel.m function x,n=Gauss_Seidel(A,b,x0) %-Gauss-Seidel 迭代法解线性方程组 %-方程组系数阵 A %-方程组右端项 b %-初始值 x0 %-求解要求的精确度 eps %-迭代步数控制 M 28 %-返回求得的解 x %-返回迭代步数 n eps=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+f n=1; %迭代次数 err=norm(x-x0,inf) while(err=eps) x0=x; x=G*x0+f n=n+1; err=norm(x-x0,inf) if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; end end 解方程组，解方程组，test7.mtest7.m clear all; 29 clc; A=5 -1 -3; -1 2 4; -3 4 15; b=-2;1;10; disp(精确解); x=Ab disp(迭代初始值); x0=0;0;0 disp(Jacobi 迭代过程：); xj,nj=Jaccobi(A,b,x0); disp(Jacobi 最终迭代结果：); xj disp(迭代次数); nj disp(Gauss-Seidel 迭代过程：); xg,ng=Gauss_Seidel(A,b,x0); disp(Gauss-Seidel 最终迭代结果：); xg disp(迭代次数); ng 运行结果：运行结果： 精确解精确解 x = 30 -0.0820 -1.8033 1.1311 迭代初始值迭代初始值 x0 = 0 0 0 JacobiJacobi 迭代过程：迭代过程： x = -0.4000 0.5000 0.6667 err = 0.6667 x = 0.1000 -1.0333 0.4533 err = 1.5333 . x = -0.0820 -1.8033 31 1.1311 err = 9.6603e-006 JacobiJacobi 最终迭代结果：最终迭代结果： xj = -0.0820 -1.8033 1.1311 迭代次数 nj = 281 Gauss-SeidelGauss-Seidel 迭代过程：迭代过程： x = -0.4000 0.3000 0.5067 err = 0.5067 x = -0.0360 -0.5313 0.8012 err = 0.8313 x = 32 -0.0256 -1.1151 0.9589 err = 0.5838 . x = -0.0820 -1.8033 1.1311 err = 9.4021e-006 Gauss-SeidelGauss-Seidel 最终迭代结果：最终迭代结果： xg = -0.0820 -1.8033 1.1311 迭代次数 ng = 20 程序：程序： 33 NewtonNewton 迭代法：迭代法：Newtoniter.mNewtoniter.m function 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:maxiter x=x0-fvalue/dfvalue err=abs(x-x0) x0=x; fvalue=subs(f,x0) dfvalue=subs(df,x0); if(err=eps) mf=subs(f,(a+b)/2); if(mf=0) x=mf;n=n+1;return end if(mf*fa0) b=(a+b)/2; else a=(a+b)/2; end iter=iter+1; end x=(a+b)/2; iter=iter+1; 求方程的实根：求方程的实根：test9.mtest9.m 38 clear all; clc; syms x f=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 程序：程序： NewtonNewton 插值：插值：Newtominter.mNewtominter.m function f=Newtominter(x,y,x0) %牛顿插值 x 插值节点 %y 为对应的函数值 %函数返回 Newton 插值多项式在 x_0 点的值 f syms t; 39 if(length(x) = length(y) n = length(x); c(1:n) = 0.0; else disp(x 和 y 的维数不相等！); return; end f = 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); end c(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); 40 end end end 用等距节点做用等距节点做 f(x)f(x)的的 NewtonNewton 插值：插值：test10.mtest10.m n1=5;n2=10;n3=15; x0=0:0.01:1; y0=sin(pi.*x0); x1=linspace(0,1,n1);%等距节点,节点数 5 y1=sin(pi.*x1); f01=Newtominter(x1,y1,x0); x2=linspace(0,1,n2);%等距节点,节点数 10 y2=sin(pi.*x2); f02 = Newtominter(x2,y2,x0); x3=linspace(0,1,n3);%等距节点,节点数 15 y3=sin(pi.*x3); f03= Newtominter(x3,y3,x0); plot(x0,y0,-r)%原图 hold on plot(x0,f01,-g)%5 个节点 plot(x0,f02,-k)%10 个节点 plot(x0,f03,-b)%15 个节点 legend(原图,5 个节点 Newton 插值多项式,10 个节点 Newton 插值多项式,15 个 节点 Newton 插值多项式) 运行结果：运行结果： 41 取不同的节点做牛顿插值。得到结果图像如下： 可以看出原图与插值多项式的图像近似重合，说明插值效果较好。 程序：程序： LagrangeLagrange 插值：插值：Lagrange.mLagrange.m function f,f0 = Lagrange(x,y,x0) %Lagrange 插值 x 为插值结点，y 为对应的函数值，x0 为要计算的点。 %函数返回 L_n(x)表达式 f 和 L_n(x0)的值 f0。 syms t; if(length(x) = length(y) n = length(x); else disp(x 和 y 的维数不相等！); return; 42 end %检错 f = 0.0; for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); %计算 Lagrange 基函数 end; f = f + l; %计算 Lagrange 插值函数 simplify(f); %化简 if(i=n) if(nargin = 3) f0 = subs(f,t,x0); %如果 3 个参数则计算插值点的函数值 else f = collect(f); %如果 2 个参数则将插值多项式展开 f = vpa(f,6); %将插值多项式的系数化成 6 位精度的小数 end end end 用等距节点做用等距节点做 LagrangeLagrange 插值：插值：test11.mtest11.m 43 clear all; clc; n1=5;n2=10;n3=15; x0=-5:0.02:5; y0=1./(1+x0.2); x1=linspace(-5,5,n1);%等距节点,节点数 5 y1=1./(1+x1.2); f1,f01 = Lagrange(x1,y1,x0); x2=linspace(-5,5,n2);%等距节点,节点数 10 y2=1./(1+x2.2); f2,f02 = Lagrange(x2,y2,x0); x3=linspace(-5,5,n3);%等距节点,节点数 15 y3=1./(1+x3.2); f3,f03 = Lagrange(x3,y3,x0); plot(x0,y0,-r)%原图 hold on plot(x0,f01,-b)%5 个节点 plot(x0,f02,-g)%10 个节点 plot(x0,f03,-k)%15 个节点 xlabel(x);ylabel(f(x),L(x); legend(原图,5 个节点 Lagrange 插值多项式,10 个节点 Lagrange 插值多项式 ,15 个节点 Lagrange 插值多项式) 运行结果：运行结果： 选取了 5,10,15 个节点做 Lagrange 插值，得到原图与插值多项式的图像如下： 44 当选取节点数较多时，选取 15 个节点时出现 Runge 现象。 程序：程序： 复合梯形求积：复合梯形求积：Comtrap.mComtrap.m function s=Comtrap(f,a,b,n) %复合梯形求积 s 积分结果 %f 积分函数 %a，b 积分区间 %n 区间个数 h=(b-a)/n; index=(a+h):h:(b-h); s1=sum(subs(f,index); s=(h/2)*(subs(f,a)+2*s1+subs(f,b); 45 复合复合 SimpsonSimpson 求积：求积： function s=Simpson(f,a,b,n) %复合 Simpson 求积 s 积分结果 %f 积分函数 %a，b 积分区间 %n 区间个数 h=(b-a)/(2*n); index1=(a+h):(2*h):(b-h); index2=(a+2*h):(2*h):(b-2*h); s1=sum(subs(f,index1); s2=sum(subs(f,index2); s=(h/3)*(subs(f,a)+4*s1+2*s2+subs(f,b); 计算计算 f(x)f(x)积分：积分：test12.mtest12.m clear all; clc; syms x f=exp(3.*x).*cos(pi.*x); a=0;b=2*pi; disp(精确积分值); I=vpa(int(f,x,a,b),10) n1=50;n2=100;n3=200;n4=500;n5=1000; disp(区间为 50,100,200,500,1000 的复合梯形积分值); T1=vpa(Comtrap(f,a,b,n1),10) et1=abs(I-T1) 46 T2=vpa(Comtrap(f,a,b,n2),10) et2=abs(I-T2) T3=vpa(Comtrap(f,a,b,n3),10) et3=abs(I-T3) T4=vpa(Comtrap(f,a,b,n4),10) et4=abs(I-T4) T5=vpa(Comtrap(f,a,b,n5),10) et5=abs(I-T5) disp(区间为 50,100,200,500,1000 的复合 Simpson 积分值); s1=vpa(Simpson(f,a,b,n1),10) es1=abs(I-s1) s2=vpa(Simpson(f,a,b,n2),10) es2=abs(I-s2) s3=vpa(Simpson(f,a,b,n3),10) es3=abs(I-s3) s4=vpa(Simpson(f,a,b,n4),10) es4=abs(I-s4) s5=vpa(Simpson(f,a,b,n5),10) es5=abs(I-s5) 运行结果：运行结果： 精确积分值精确积分值 I = 35232483.36 复合梯形与复合 Simpson 求积与误差 47 区间个数 n复合梯形求积 T误差 eT复合 Simpson 求 积 误差 eS 5035125341.19107142.1735231407.871075.49 10035204891.227592.1635232416.2467.12 20035225534.986948.3835232479.174.19 50035231369.371113.9935232483.250.11 100035232204.78278.5835232483.360.0 程序：程序： GaussLegendreGaussLegendre 求积：求积：GaussLegendre.mGaussLegendre.m function s = GaussLegendre(f,a,b,n) %-Gauss-Legendre 公式求积分，只给出 2-5 个求积结点 %f 被积函数 %b,a 积分上下限 %n 求积结点个数 n+1 %s 积分结果 ta=(b-a)/2;tb=(a+b)/2; switch n case 1, s=ta*(subs(sym(f),findsym(sym(f),ta*0.5773503+tb)+. subs(sym(f),findsym(sym(f),-ta*0.5773503+tb); case 2, 48 s=ta*(0.5555556*subs(sym(f),findsym(sym(f),ta*0.7745967+tb)+. 0.5555556*subs(sym(f),findsym(sym(f),-ta*0.7745967+tb)+. 0.8888889*subs(sym(f),findsym(sym(f),tb); case 3, s=ta*(0.3478548*subs(sym(f),findsym(sym(f),ta*0.8611363+tb)+. 0.3478548*subs(sym(f),findsym(sym(f),-ta*0.8611363+tb)+. 0.6521452*subs(sym(f),findsym(sym(f),ta*0.3399810+tb)+. 0.6521452*subs(sym(f),findsym(sym(f),-ta*0.3399810+tb); case 4, s=ta*(0.2369269*subs(sym(f),findsym(sym(f),ta*0.9061798+tb)+. 0.2369269*subs(sym(f),findsym(sym(f),-ta*0.9061798+tb)+. 0.4786287*subs(sym(f),findsym(sym(f),ta*0.5384693+tb)+. 0.4786287*subs(sym(f),findsym(sym(f),-ta*0.5384693+tb)+. 0.5688889*subs(sym(f),findsym(sym(f),tb); end GaussChebyshevGaussChebyshev 求积：求积：GaussChebyshev.mGaussChebyshev.m function s=GaussChebyshev(f,n) %GaussChebyshev 求积 s 积分值 %f 积分函数 % 求积结点个数 n+1 49 k=0:n; x=cos(2.*k+1).*pi/(2.*(n+1); a=pi/(n+1); s=a*sum(subs(f,x); 计算计算 f(x)f(x)积分：积分：test13.mtest13.m clear all; clc; syms x f1=x.2; p=sqrt(1-x.2); f2=sin(x)./x; n1=2;n2=3;n3=5; a=0;b=pi/2; disp(第一个实际积分值); I1=vpa(int(f1/p,x,-1,1),8) disp(第二个实际积分值); I2=vpa(int(f2,x,a,b),8) disp(第一个 2,3,5 点 Guass 型积分); I2=vpa(GaussChebyshev(f1,n1-1),8) I3=vpa(GaussChebyshev(f1,n2-1),8) I5=vpa(GaussChebyshev(f1,n3-1),8) disp(第二个 2,3,5 点 Guass 型积分); GL2= vpa(GaussLegendre(f2,a,b,n1-1),8) GL3= vpa(GaussLegendre(f2,a,b,n2-1),8) 50 GL5= vpa(GaussLegendre(f2,a,b,n3-1),8) 运行结果：运行结果： （1）实际积分值：I1 =1.5707963 2,3,5 点GaussChebyshev型积分 I2 =1.5707963 I3 = 1.5707963 I5 =1.5707963 （2）实际积分值 I2 =1.3707622 2,3,5 点GaussLegendre型积分 GL2 =1.370419 GL3 = 1.3707635 GL5 =1.3707622 程序：程序： EulerEuler 法：法：Euler.mEuler.m function T=Euler(f,x0,xn,y0,h) %Euler 法 %x,y 微分方程的解 %f 微分方程 %x0,y0 初始值 %xn 端点 51 %h 步长 n=(xn-x0)/h; %区间个数 x=zeros(1,n+1); y=zeros(1,n+1); x(1)=x0; y(1)=y0; for k=1:n y(k+1)=y(k)+h*feval(f,x(k),y(k); x(k+1)=x(k)+h; end T=x,y; 改进改进 EulerEuler 法：法：TranEuler.mTranEuler.m function T=TranEuler(f,x0,xn,y0,h) %改进 Euler 法 %x,y 微分方程的解 %f 微分方程 %x0,y0 初始值 %xn 端点 %h 步长 n=(xn-x0)/h; %区间个数 x=zeros(1,n+1); y=zeros(1,n+1); x(1)=x0; 52 y(1)=y0; for k=1:n x(k+1)=x(k)+h; z0=y(k)+h*feval(f,x(k),y(k); y(k+1)=y(k)+h/2*(feval(f,x(k),y(k)+feval(f,x(k+1),z0); end T=x,y; 经典经典 Runge-KuttaRunge-Kutta 法：法：R_K4.m function R=R_K4(f,x0,xn,y0,h) %经典，4 阶 RungeKutta 法 %x,y 微分方程的解 %f 微分方程 %x0,y0 初始值 %xn 端点 %h 步长 n=(xn-x0)/h; %区间个数 y=zeros(1,n+1); y(1)=y0