大连理工大学矩阵与数值分析报告上机作业.docx
- 文档编号:17471053
- 上传时间:2023-07-26
- 格式:DOCX
- 页数:51
- 大小:440.45KB
大连理工大学矩阵与数值分析报告上机作业.docx
《大连理工大学矩阵与数值分析报告上机作业.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵与数值分析报告上机作业.docx(51页珍藏版)》请在冰点文库上搜索。
大连理工大学矩阵与数值分析报告上机作业
矩阵与数值分析上机作业
学校:
大连理工大学
学院:
班级:
姓名:
学号:
授课老师:
注:
编程语言Matlab
程序:
Norm.m函数
functions=Norm(x,m)
%求向量x的范数
%m取1,2,inf分别表示1,2,无穷范数
n=length(x);
s=0;
switchm
case1%1-范数
fori=1:
n
s=s+abs(x(i));
end
case2%2-范数
fori=1:
n
s=s+x(i)^2;
end
s=sqrt(s);
caseinf%无穷-范数
s=max(abs(x));
end
计算向量x,y的范数
Test1.m
clearall;
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
clearall;
clc;
n=100;%区间
h=2*10^(-15)/n;%步长
x=-10^(-15):
h:
10^(-15);
%第一种原函数
f1=zeros(1,n+1);
fork=1:
n+1
ifx(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);
fork=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
functiony=QinJS(a,x)
%y输出函数值
%a多项式系数,由高次到零次
%x给定点
n=length(a);
s=a
(1);
fori=2:
n
s=s*x+a(i);
end
y=s;
计算p(x):
test3.m
clearall;
clc;
x=1.6:
0.2:
2.4;%x=2的邻域
disp('x=2的邻域:
');x
a=[1-18144-6722016-40325376-46082304-512];
p=zeros(1,5);
fori=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;
fork=1:
nk
pk(k)=QinJS(a,xk(k));
end
plot(xk,pk,'-r');
xlabel('x');ylabel('p(x)');
运行结果:
x=2的邻域:
x=
1.60001.80002.00002.20002.4000
相应多项式p值:
p=
1.0e-003*
-0.2621-0.000500.00050.2621
p(x)在
[1.95,20.5]上的图像
程序:
LU分解,LUDecom.m
function[L,U]=LUDecom(A)
%不带列主元的LU分解
N=size(A);
n=N
(1);
L=eye(n);U=zeros(n);
fori=1:
n
U(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);
end
fori=2:
n
forj=i:
n
z=0;
fork=1:
i-1
z=z+L(i,k)*U(k,j);
end
U(i,j)=A(i,j)-z;
end
forj=i+1:
n
z=0;
fork=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分解,PLUDecom.m
function[P,L,U]=PLUDecom(A)
%带列主元的LU分解
[m,m]=size(A);
U=A;
P=eye(m);
L=eye(m);
fori=1:
m
forj=i:
m
t(j)=U(j,i);
fork=1:
i-1
t(j)=t(j)-U(j,k)*U(k,i);
end
end
a=i;b=abs(t(i));
forj=i+1:
m
ifb b=abs(t(j)); a=j; end end ifa~=i forj=1: m c=U(i,j); U(i,j)=U(a,j); U(a,j)=c; end forj=1: m c=P(i,j); P(i,j)=P(a,j); P(a,j)=c; end c=t(a); t(a)=t(i); t(i)=c; end U(i,i)=t(i); forj=i+1: m U(j,i)=t(j)/t(i); end forj=i+1: m fork=1: i-1 U(i,j)=U(i,j)-U(i,k)*U(k,j); end end end L=tril(U,-1)+eye(m); U=triu(U,0); (1) (2)程序: Test4.m clearall; clc; forn=5: 30 x=zeros(n,1); A=-ones(n); A(: n)=ones(n,1); fori=1: n A(i,i)=1; forj=(i+1): (n-1) A(i,j)=0; end x(i)=1/i; end disp('当n=');disp(n); disp('方程精确解: '); x b=A*x;%系数b disp('利用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); fori=1: n InvAL(: i)=U\(L\I(: i)); end disp('利用LU分解的A的逆矩阵: '); InvAL End 运行结果: (1)只列出n=5,6,7的结果 当n=5 方程精确解: x=1.0000 0.5000 0.3333 0.2500 0.2000 利用LU分解方程组的解: xLU= 1.0000 0.5000 0.3333 0.2500 0.2000 利用PLU分解方程组的解: xPLU= 1.0000 0.5000 0.3333 0.2500 0.2000 当n=6 方程精确解: x= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 利用LU分解方程组的解: xLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 利用PLU分解方程组的解: xPLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 当n=7 方程精确解: x= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 利用LU分解方程组的解: xLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 利用PLU分解方程组的解: xPLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 (2)只列出n=5,6,7时A的逆矩阵的结果 当n=5 A的准确逆矩阵: InvA= 0.5000-0.2500-0.1250-0.0625-0.0625 00.5000-0.2500-0.1250-0.1250 000.5000-0.2500-0.2500 0000.5000-0.5000 0.50000.25000.12500.06250.0625 利用LU分解的A的逆矩阵: InvAL= 0.5000-0.2500-0.1250-0.0625-0.0625 00.5000-0.2500-0.1250-0.1250 000.5000-0.2500-0.2500 0000.5000-0.5000 0.50000.25000.12500.06250.0625 当n=6 A的准确逆矩阵: InvA= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0313 00.5000-0.2500-0.1250-0.0625-0.0625 000.5000-0.2500-0.1250-0.1250 0000.5000-0.2500-0.2500 00000.5000-0.5000 0.50000.25000.12500.06250.03130.0313 利用LU分解的A的逆矩阵: InvAL= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0313 00.5000-0.2500-0.1250-0.0625-0.0625 000.5000-0.2500-0.1250-0.1250 0000.5000-0.2500-0.2500 00000.5000-0.5000 0.50000.25000.12500.06250.03130.0313 当n=7 A的准确逆矩阵: InvA= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156 00.5000-0.2500-0.1250-0.0625-0.0313-0.0313 000.5000-0.2500-0.1250-0.0625-0.0625 0000.5000-0.2500-0.1250-0.1250 00000.5000-0.2500-0.2500 000000.5000-0.5000 0.50000.25000.12500.06250.03130.01560.0156 利用LU分解的A的逆矩阵: InvAL= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156 00.5000-0.2500-0.1250-0.0625-0.0313-0.0313 000.5000-0.2500-0.1250-0.0625-0.0625 0000.5000-0.2500-0.1250-0.1250 00000.5000-0.2500-0.2500 000000.5000-0.5000 0.50000.25000.12500.06250.03130.01560.0156 程序: Cholesky分解: Cholesky.m functionL=Cholesky(A) N=size(A); n=N (1); L=zeros(n); L(1,1)=sqrt(A(1,1)); fori=2: n L(i,1)=A(i,1)/L(1,1); end forj=2: n s1=0; fork=1: j-1 s1=s1+L(j,k)^2; end L(j,j)=sqrt(A(j,j)-s1); fori=j+1: n s2=0; fork=1: j-1 s2=s2+L(i,k)*L(j,k); end L(i,j)=(A(i,j)-s2)/L(j,j); end end 计算Ax=b;Test5.m clearall;clc; forn=10: 20 A=zeros(n,n); b=zeros(n,1); fori=1: n forj=1: n A(i,j)=1/(i+j-1); end b(i,1)=i; end disp('n=');disp(n); disp('方程组原始解');x0=A\b disp('利用Cholesky分解的方程组的解'); L=Cholesky(A) x=L'\(L\b) end 运行结果: 只列出了n=10,11的结果 n=10 方程组原始解 x0= 1.0e+008* -0.0000 0.0010 -0.0233 0.2330 -1.2108 3.5947 -6.3233 6.5114 -3.6233 0.8407 利用Cholesky分解的方程组的解 x= 1.0e+008* -0.0000 0.0010 -0.0233 0.2330 -1.2105 3.5939 -6.3219 6.5100 -3.6225 0.8405 n= 11 方程组原始解 x0= 1.0e+009* 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4039 3.2863 -4.7869 4.2260 -2.0685 0.4305 利用Cholesky分解的方程组的解 x= 1.0e+009* 0.0000 -0.0002 0.0046 -0.0563 0.3668 -1.3972 3.2716 -4.7669 4.2094 -2.0608 0.4290 程序: (1)House.m functionu=House(x) n=length(x); e1=eye(n,1); w=x-norm(x,2)*e1; u=w/norm(w,2); (2)Hou_A.m functionHA=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.m clearall; clc; A=[1234; -12sqrt (2)sqrt(3); -22exp (1)pi; -sqrt(10)2-37; 0275/2]; HA=Hou_A(A) 运行结果: H= 0.2500-0.2500-0.5000-0.79060 -0.25000.9167-0.1667-0.26350 -0.5000-0.16670.6667-0.52700 -0.7906-0.2635-0.52700.16670 00001.0000 HA= 4.0000-2.58111.4090-6.5378 0.00000.47300.8839-1.7805 0.0000-1.05411.6576-3.8836 0.0000-2.8289-4.6770-4.1078 02.00007.00002.5000 程序: Jacobi迭代: Jaccobi.m function[x,n]=Jaccobi(A,b,x0) %--·方程组系数阵A %--·方程组右端顶b %--初始值x0 %--求解要求精确度eps %--迭代步数控制M %--·返回求得的解x %--·返回迭代步数n M=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+f n=1;%迭代次数 err=norm(x-x0,inf) while(err>=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 clearall; clc; A=[5-1-3; -124; -3415]; 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.583
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 数值 分析 报告 上机 作业