大连理工大学矩阵与数值分析上机作业.docx
《大连理工大学矩阵与数值分析上机作业.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业.docx（51页珍藏版）》请在冰豆网上搜索。
大连理工大学矩阵与数值分析上机作业
矩阵与数值分析上机作业
学校：
大连理工大学
学院:
班级:
姓名:
学号：
授课老师:
注：
编程语言Matlab
1.琴虑计算给定働量的葩址输入向量広』（巾斑…宀产输出||工||“||工|怙㈣心请编制一牛通用程序，并用你編制的程序计算如下询量的范数：
对网1加，wm甚至更大的“计算其范数，你会发现什幺结粟？
你能否修改你的程序使得计算绪果相时赫■确呢？
程序:
Norm.m函数
functions=Norm（x,m）
%求向量x的范数
%mx1,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（if2;
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时
范数:
1
-范数:
10
-范数:
1
x的1-范数29290;x的2-范数：
1.2449;x的无穷-
y的1-范数:
55;y的2-范数:
19.6214;y的无穷n=100时
x的1-范数:
5.1874;x的2-范数：
1.2787;x的无穷
y的1-范数:
5050；
的2-范数:
581.6786；y的无穷-范数:
100
n=1000时
x的1-范数74855;x的2-范数：
1.2822;x的无穷-范数:
1
y的1-范数:
500500;y的2-范数：
1.8271e+004;y的无穷-范数:
1000
2.耆虑砂==呼^其中定51/（0）=此时几期是连绽函戟.用此公式计算当工“―1旷巾U）-缪时的函数值*風出图像.另一方面*哮虑下面算法：
d1+j
1/（/=1tbfjj
1/=1
仙
y=liid/（d—1（
endif
用此算法计％€[-10-0io_is]时的圉数血画出图像.比校一下岌生了什么？
程序
Test2.m
clearall;
clc;
n=100;%区间
h=2*10A（-15）/n;%步长
x=-10A（-15）:
h:
10A（-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;
endend
subplot（2,1,1）;
plot（x,f1,'-r'）;
axis（[-10A（-15）,10A（-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（[-10A（-15）,10A（-15）,-1,2]）;
legend（'第二种算法'）;
运行结果：
農IQ
显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当X[1015,1015]时，
d1x,计算机进行舍入造成d恒等于1,结果函数值恒为1。
3.首先编制一个別用秦丸昶律法计算一个多项式在给定点的函數值的遇用程序，你
的程序邑括揄入多项式的系就以及给定点，嶽出函数值.利用你騙制的粒序计算
p（x）=（^-2）s=泸一18rs+IMj-7-672^+2nitkrs-lO32r4+越芮护-460SX2+230b-512在工=2邻域附近的值-画出此r）在jt€|1-航20-司上的图像・
程序：
秦九韶算法：
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）在X[1.95,20.5]上的图像
4,鎬制计雾给定拒阵冲的U7分解和FL17分解的通用程序、然后用你编制的程序完咸下面两个计鼻任务：
（1）考虑
II
（I
-1
自己取定上eRS并计鼻古=屉.然后用你编制的不选主元和列主无的伽u抽消去法求解谏方程组，记你计算韭的解为氏对和从5刘盹估计计界解的精度*
⑵对仕从訐到30计算其逆雑阵.
程序：
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）;
U（i,j）=A（i,j）-z;
end
forj=i+1:
n
z=0;
fork=1:
i-1z=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
a=i;b=abs（t（i））;
forj=i+1:
m
ifbend
end
ifa~=i
forj=1:
mc=U（i,j）;U（i,j）=U（a,j）;U（a,j）=c;
end
forj=1:
mc=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）;
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.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.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
00.5000-0.2500-0.2500
00.5000-0.5000
0.50000.25000.12500.06250.03130.01560.0156
5.蜻制计算对称正定薛的OioteM洽解的遇用程序，幷用你编制的程序计算如三乩其中丸=（aij）€即".旳=b可以由你自己取忌对必从10到凶验证程序的可命出
程序：
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）A2;
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
-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.0046
-0.0563
0.3668
-1.3972
3.2716
-4.7669
4.2094
-2.0608
0.4290
6.
（1）编制强序其作用是时辅:
入的向董「输出单位向量“彳吏得（mN）』=
⑵编制换薛丹=
序并怎显式的计算出
（3）考虑距阵
汀—加十丘朗旳乘EMW肿山的程序厅扎注意、你的程
1234、
-1乂暑亞
A=-22ff,
-7102-37
k0275/2；
用你编制的程序计算H使得的第一列为“I的形式，并将乩4的结果显示.
程序：
（1）House.m
functionu=House（x）
n=length（x）;
e1=eye（n,1）;
w=x-norm（x,2）*e1;
u=w/norm（w,2）;
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.7906
0
-0.2500
0.9167
-0.1667
-0.2635
0
-0.5000
-0.1667
0.6667
-0.5270
0
-0.7906
-0.2635
-0.5270
0.1667
0
0
0
00
1.0000
HA=
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
7.用代求解下面的方程组*输出迭代每一步的误差||斗—斗||：
5xj—X2—扫=—2
«—叼十2^2十4^3=1
—3g+4电4-15^3=10
程序:
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））;
L=-tril（A,-1）;
U=-triu（A,1）;
%求A的对角矩阵
%求A的下三角阵
%求A的上三角阵
J=D\（L+U）;
f=D\b;
x=J*x0+fn=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.mfunction[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）;
升级会员 