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