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

《大连理工大学矩阵与数值分析上机作业.doc》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业.doc（33页珍藏版）》请在装配图网上搜索。

矩阵与数值分析上机作业 学校： 大连理工大学 学院： 班级： 姓名： 学号： 授课老师： 注：编程语言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的范数 Test1.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(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的无穷-范数:1 y的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 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)/(d-1); else f2(k)=1; end end subplot(2,1,2); plot(x,f2,-r); axis([-10^(-15),10^(-15),-1,2]); legend(第二种算法); 运行结果： 显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当时，，计算机进行舍入造成恒等于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); 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,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)在[1.95,20.5]上的图像 程序： LU分解，LUDe..m 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,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分解，PLUDe..m 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 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 %--求解要求的精确度 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: 迭代次数太多，可能不收敛！); return; end end 解方程组，test7.m clear all; clc; A=[5 -1 -3; -1 2 4; -3 4 15]; b=[-2;1;10]; disp(精确解); x=A\b 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 = -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.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) %牛顿法 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*fa<0) 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) 运行结果： [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) == 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); 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 = 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插值多项式) 运行结果： 取不同的节点做牛顿插值。得到结果图像如下： 可以看出原图与插值多项式的图像近似重合，说明插值效果较好。 程序： 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 = 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 用等距节点做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 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插值，得到原图与插值多项式的图像如下： 当选取节点数较多时，选取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 区间个数 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;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的复合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求积与误差 区间个数n 复合梯形求积T 误差eT 复合Simpson求积 误差eS 50 35125341.19 107142.17 35231407.87 1075.49 100 35204891.2 27592.16 35232416.24 67.12