《大连理工大学矩阵与数值分析上机作业.pdf》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业.pdf（33页珍藏版）》请在三一文库上搜索。
1、矩阵与数值分析上机作业 学校：大连理工大学 学院： 班级： 姓名： 学号： 授课老师： 注：编程语言 Matlab 程序： Norm.m函数 function s=Norm(x,m) % 求向量 x 的范数 %m 取 1,2,inf分别表示 1,2 ，无穷范数 n=length(x); 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，y 的范数 Tes
2、t1.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); disp( y 的 1- 范数 : );disp(Norm(y1,1); disp( y 的 2- 范数 : );disp(Norm(y
3、1,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
4、 时 ); 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 的无穷 - 范数 :1 y 的
5、 1- 范数 :55 ； y的 2-范数 :19.6214 ； y 的无穷 - 范数 :10 n=100 时 x 的 1- 范数 :5.1874 ；x 的 2-范数 : 1.2787； x 的无穷 - 范数 :1 y 的 1- 范数 :5050 ； y的 2-范数 :581.6786 ； y 的无穷 - 范数 :100 n=1000 时 x 的 1- 范数 :7.4855 ； x 的 2-范数 :1.2822 ； x的无穷 - 范数 :1 y 的 1- 范数 : 500500 ； y 的 2- 范数 :1.8271e+004 ；y 的无穷 - 范数 :1000 程序 Test2.m clear
6、 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; 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)
7、/(d-1); else f2(k)=1; end end subplot(2,1,2); plot(x,f2,-r); axis(-10(-15),10(-15),-1,2); legend( 第二种算法 ); 运行结果： 显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当10,10 1515 x时，xd 1， 计算机进行舍入造成d恒等于 1，结果函数值恒为1。 程序： 秦九韶算法 :QinJS.m function y=QinJS(a,x) %y输出函数值 %a多项式系数，由高次到零次 %x给定点 n=length(a); s=a(1); for i=2:n s=s*x+a(i)
8、; end y=s; 计算 p(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
9、,pk,-r); xlabel(x);ylabel(p(x); 运行结果： 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) 在x1.95,20.5上的图像 程序： LU分解，LUDem function L,U=LUDe.(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,
10、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 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 PLU分解，PLUDem function P,L,U =PLUDe.(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
11、 k=1:i-1 t(j)=t(j)-U(j,k)*U(k,i); end end 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_Seidel迭代： Gauss_Seidel.m function x,n=Gauss_Seidel(A,b,x0) %-Gauss-Seidel迭代法解线性方程组 %-方程组系数阵 A %-方程组右端项 b %-初始值 x0 %-求
12、解要求的精确度 eps %-迭代步数控制 M %-返回求得的解 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: 迭代次数太多，可能不收敛！ );
13、 return; end end 解方程组， test7.m clear all ; 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最终迭代结果： );
14、xg disp( 迭代次数 ); ng 运行结果： 精确解 x = -0.0820 -1.8033 1.1311 迭代初始值 x0 = 0 0 0 Jacobi 迭代过程： 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 1.1311 err = 9.6603e-006 Jacobi 最终迭代结果： xj = -0.0820 -1.8033 1.1311 迭代次数 nj = 281 Gauss-Seidel迭代过程： x = -0.4000 0
15、.3000 0.5067 err = 0.5067 x = -0.0360 -0.5313 0.8012 err = 0.8313 x = -0.0256 -1.1151 0.9589 err = 0.5838 . x = -0.0820 -1.8033 1.1311 err = 9.4021e-006 Gauss-Seidel最终迭代结果： xg = -0.0820 -1.8033 1.1311 迭代次数 ng = 20 程序： Newton 迭代法：Newtoniter.m function x,iter,fvalue=Newtoniter(f,df,x0,eps,maxiter) % 牛
16、顿法 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 (m
17、f*fa0) b=(a+b)/2; else a=(a+b)/2; end iter=iter+1; end x=(a+b)/2; iter=iter+1; 求方程的实根：test9.m 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) 运行结果
18、： 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.m function f=Newtominter(x,y,x0) % 牛顿插值 x 插值节点 %y为对应的函数值 % 函数返回Newton 插值多项式在x_0 点的值 f syms t ; if (length(x) = len
19、gth(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); els
20、e% 如果 2 个参数则直接给出插值多项式 f = collect(f); % 将插值多项式展开 f = vpa(f, 6); end end end 用等距节点做f(x) 的 Newton 插值： test10.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 = Newtomin
21、ter(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 插值多项式 ) 运行结果： 取不同的节点做牛顿插值。得到结果图像如下： 可以看
22、出原图与插值多项式的图像近似重合，说明插值效果较好。 程序： Lagrange 插值： Lagrange.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; end% 检错 f = 0.0; for (i = 1:n) l = y(i); for (j =
23、 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 en
24、d end 用等距节点做Lagrange 插值： test11.m 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 y
25、3=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 插值，得到原图与插值多项式的图像如下
26、： 当选取节点数较多时，选取15 个节点时出现Runge现象。 程序： 复合梯形求积：.trap.m function s=.trap(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); 复合 Simpson 求积： function s=Simpson(f,a,b,n) % 复合 Simpson 求积 s 积分结果 %f 积分函数 %a ， b 积分区间 %n 区
27、间个数 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)积分： test12.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;
28、n4=500;n5=1000; disp( 区间为 50,100,200,500,1000的复合梯形积分值 ); T1=vpa(.trap(f,a,b,n1),10) et1=abs(I-T1) T2=vpa(.trap(f,a,b,n2),10) et2=abs(I-T2) T3=vpa(.trap(f,a,b,n3),10) et3=abs(I-T3) T4=vpa(.trap(f,a,b,n4),10) et4=abs(I-T4) T5=vpa(.trap(f,a,b,n5),10) et5=abs(I-T5) disp( 区间为 50,100,200,500,1000的复合 Simps
29、on 积分值 ); 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 求积与误差 区间个数 n 复合梯形求积T 误差 eT 复合 Si
30、mpson求积误差 eS 50 35125341.19 107142.17 35231407.87 1075.49 100 35204891.2 27592.16 35232416.24 67.12 200 35225534.98 6948.38 35232479.17 4.19 500 35231369.37 1113.99 35232483.25 0.11 1000 35232204.78 278.58 35232483.36 0.0 程序： GaussLegendre 求积： GaussLegendre.m function s = GaussLegendre(f,a,b,n) %-Ga
31、uss-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, s=ta*(0.5555556*subs(sym(f),findsym(sym(f),ta*0.7745967+tb)+. 0.5555556*subs(sym(f)
32、,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),-t
33、a*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); e
34、nd GaussChebyshev 求积： GaussChebyshev.m function s=GaussChebyshev(f,n) %GaussChebyshev求积 s 积分值 %f 积分函数 % 求积结点个数n+1 k=0:n; x=cos(2.*k+1).*pi/(2.*(n+1); a=pi/(n+1); s=a*sum(subs(f,x); 计算 f(x)积分： test13.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( 第一个实
35、际积分值 ); 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
36、(f2,a,b,n2-1),8) 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 程序： Euler法： Euler.m function T=Euler(f,x0,xn,y0,h) %Euler
37、法 %x,y 微分方程的解 %f 微分方程 %x0,y0 初始值 %xn端点 %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; 改进 Euler 法：TranEuler.m function T=TranEuler(f,x0,xn,y0,h) % 改进 Euler法 %x,y 微分方程的解 %f 微分方程 %x0,y0 初始值 %xn端点 %h 步长 n=(
38、xn-x0)/h; % 区间个数 x=zeros(1,n+1); y=zeros(1,n+1); x(1)=x0; 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-Kutta法： R_K4.m function R=R_K4(f,x0,xn,y0,h) % 经典， 4 阶 RungeKutta法 %x,y 微分方程的解 %f 微分方程 %x0,y0 初始值 %xn端点
39、%h步长 n=(xn-x0)/h; % 区间个数 y=zeros(1,n+1); y(1)=y0; x=x0:h:xn; for k=1:n k1=feval(f,x(k),y(k); k2=feval(f,x(k)+h/2,y(k)+h/2*k1); k3=feval(f,x(k)+h/2,y(k)+h/2*k2); k4=feval(f,x(k)+h,y(k)+h*k3); y(k+1)=y(k)+h/6*(k1+2*k2+2*k3+k4); end R=x,y; 解该方程： test14.m clear all ； clc; t0=0;x0=2;tn=1; h1=0.1;h2=0.01;
40、h3=0.001; E1=Euler( T14fun,t0,tn,x0,h1); TE1=TranEuler(T14fun,t0,tn,x0,h1); R1=R_K4(T14fun,t0,tn,x0,h1); E2=Euler( T14fun,t0,tn,x0,h2); TE2=TranEuler(T14fun,t0,tn,x0,h2); R2=R_K4(T14fun,t0,tn,x0,h2); E3=Euler( T14fun,t0,tn,x0,h3); TE3=TranEuler(T14fun,t0,tn,x0,h3); R3=R_K4(T14fun,t0,tn,x0,h3); subpl
41、ot(3,1,1) plot(E1(:,1),E1(:,2),*g) hold on plot(TE1(:,1),TE1(:,2),.r) plot(R1(:,1),R1(:,2),ob ) xlabel(t);ylabel(x(t); axis(0 1 0 2); title( 步长 0.1 的 Euler法，改进Euler 法， Runge_Kutta 法 ); legend( Euler, 改进 Euler, Runge-Kutta) subplot(3,1,2) plot(E2(:,1),E2(:,2),*g) hold on plot(TE2(:,1),TE2(:,2),.r) pl
42、ot(R2(:,1),R2(:,2),ob ) xlabel(t);ylabel(x(t); axis(0 1 0 2); title( 步长 0.01 的 Euler法，改进Euler法， Runge_Kutta法 ); legend( Euler, 改进 Euler, Runge-Kutta) subplot(3,1,3) plot(E3(:,1),E3(:,2),*g) hold on plot(TE3(:,1),TE3(:,2),.r) plot(R3(:,1),R3(:,2),ob ) xlabel(t);ylabel(x(t); axis(0 1 0 2); title( 步长 0.001 的 Euler法，改进Euler法， Runge_Kutta 法 ); legend( Euler, 改进 Euler, Runge-Kutta) 其中调用的函数 :T14fun.m function z=T14fun(t,x) z=1/(t+1)2*(t*x-x2); 运行结果图：