大连理工大学矩阵与数值分析上机作业18478.docx
《大连理工大学矩阵与数值分析上机作业18478.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业18478.docx（52页珍藏版）》请在冰豆网上搜索。
大连理工大学矩阵与数值分析上机作业18478
矩阵与数值分析上机作业
学校：
大连理工大学
学院：
班级：
姓名：
学号：
授课老师：
注：
编程语言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分解，LUDe..m
function[L,U]=LUDe.（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分解，PLUDe..m
function[P,L,U]=PLUDe.（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]=LUDe.（A）;%LU分解
xLU=U\（L\b）
disp（'利用PLU分解方程组的解:
'）;
[P,L,U]=PLUDe.（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
升级会员 