大连理工大学矩阵与数值分析上机作业.pdf
《大连理工大学矩阵与数值分析上机作业.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 插值多项式 ) 运行结果: 取不同的节点做牛顿插值。得到结果图像如下: 可以看
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 数值 分析 上机 作业
链接地址:https://www.31doc.com/p-4680510.html