大连理工大学矩阵与数值分析上机作业.docx
《大连理工大学矩阵与数值分析上机作业.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵与数值分析上机作业.docx(45页珍藏版)》请在冰豆网上搜索。
大连理工大学矩阵与数值分析上机作业
矩阵与数值分析上机作业
学校:
大连理工大学
学院:
班级:
姓名:
学号:
授课老师:
注:
编程语言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
ifbb=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.801