大连理工大学矩阵与数值分析上机作业Word文档格式.docx
《大连理工大学矩阵与数值分析上机作业Word文档格式.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业Word文档格式.docx（50页珍藏版）》请在冰豆网上搜索。
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:
y2=[1:
y3=[1:
disp（'
n=10时'
）;
x的1-范数:
'
disp（Norm（x1,1））;
x的2-范数:
disp（Norm（x1,2））;
x的无穷-范数:
disp（Norm（x1,inf））;
y的1-范数:
disp（Norm（y1,1））;
y的2-范数:
disp（Norm（y1,2））;
y的无穷-范数:
disp（Norm（y1,inf））;
n=100时'
disp（Norm（x2,1））;
disp（Norm（x2,2））;
disp（Norm（x2,inf））;
disp（Norm（y2,1））;
disp（Norm（y2,2））;
disp（Norm（y2,inf））;
n=1000时'
disp（Norm（x3,1））;
disp（Norm（x3,2））;
disp（Norm（x3,inf））;
disp（Norm（y3,1））;
disp（Norm（y3,2））;
disp（Norm（y3,inf））;
运行结果：
n=10时
2.9290；
1.2449；
x的无穷-范数:
1
55；
y的2-范数:
19.6214；
y的无穷-范数:
10
n=100时
5.1874；
1.2787；
5050；
y的2-范数:
581.6786；
100
n=1000时
7.4855；
x的2-范数:
1.2822；
x的无穷-范数:
500500；
y的2-范数:
1.8271e+004；
1000
程序
Test2.m
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）;
d=1+x（k）;
if（d~=1）
f2（k）=log（d）/（d-1）;
f2（k）=1;
subplot（2,1,2）;
plot（x,f2,'
第二种算法'
显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当时，，计算机进行舍入造成恒等于1，结果函数值恒为1。
秦九韶算法:
QinJS.m
functiony=QinJS（a,x）
%y输出函数值
%a多项式系数，由高次到零次
%x给定点
n=length（a）;
s=a
（1）;
fori=2:
s=s*x+a（i）;
y=s;
计算p（x）:
test3.m
x=1.6:
0.2:
2.4;
%x=2的邻域
x=2的邻域：
x
a=[1-18144-6722016-40325376-46082304-512];
p=zeros（1,5）;
fori=1:
5
p（i）=QinJS（a,x（i））;
相应多项式p值：
p
xk=1.95:
0.01:
20.5;
nk=length（xk）;
pk=zeros（1,nk）;
k=1;
nk
pk（k）=QinJS（a,xk（k））;
plot（xk,pk,'
xlabel（'
x'
ylabel（'
p（x）'
x=
1.60001.80002.00002.20002.4000
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）;
U（1,i）=A（1,i）;
L（i,1）=A（i,1）/U（1,1）;
forj=i:
z=0;
fork=1:
i-1
z=z+L（i,k）*U（k,j）;
end
U（i,j）=A（i,j）-z;
forj=i+1:
z=z+L（j,k）*U（k,i）;
L（j,i）=（A（j,i）-z）/U（i,i）;
PLU分解，PLUDecom.m
function[P,L,U]=PLUDecom（A）
%带列主元的LU分解
[m,m]=size（A）;
U=A;
P=eye（m）;
L=eye（m）;
m
t（j）=U（j,i）;
t（j）=t（j）-U（j,k）*U（k,i）;
a=i;
b=abs（t（i））;
ifb<
abs（t（j））
b=abs（t（j））;
a=j;
ifa~=i
forj=1:
c=U（i,j）;
U（i,j）=U（a,j）;
U（a,j）=c;
c=P（i,j）;
P（i,j）=P（a,j）;
P（a,j）=c;
c=t（a）;
t（a）=t（i）;
t（i）=c;
U（i,i）=t（i）;
U（j,i）=t（j）/t（i）;
U（i,j）=U（i,j）-U（i,k）*U（k,j）;
L=tril（U,-1）+eye（m）;
U=triu（U,0）;
（1）
（2）程序：
Test4.m
forn=5:
30
x=zeros（n,1）;
A=-ones（n）;
A（:
n）=ones（n,1）;
fori=1:
A（i,i）=1;
forj=（i+1）:
（n-1）
A（i,j）=0;
x（i）=1/i;
disp（'
当n='
disp（n）;
方程精确解:
x
b=A*x;
%系数b
利用LU分解方程组的解:
[L,U]=LUDecom（A）;
%LU分解
xLU=U\（L\b）
利用PLU分解方程组的解:
[P,L,U]=PLUDecom（A）;
%PLU分解
xPLU=U\（L\（P\b））
%求解A的逆矩阵
A的准确逆矩阵:
InvA=inv（A）
InvAL=zeros（n）;
%利用LU分解求A的逆矩阵
I=eye（n）;
n
InvAL（:
i）=U\（L\I（:
i））;
利用LU分解的A的逆矩阵:
InvAL
End
（1）只列出n=5,6,7的结果
升级会员 