大连理工大学矩阵与数值研究分析上机作业Word文档格式.docx
《大连理工大学矩阵与数值研究分析上机作业Word文档格式.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值研究分析上机作业Word文档格式.docx（45页珍藏版）》请在冰豆网上搜索。
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；
581.6786；
100
n=1000时
7.4855；
x的2-范数:
1.2822；
500500；
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
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）;
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）;
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的结果
当n=5
x=1.0000
0.5000
0.3333
0.2500
0.2000
xLU=
1.0000
xPLU=
当n=6
0.1667
当n=7
0.1429
（2）只列出n=5,6,7时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
InvAL=
0.50000.25000.12500.06250.0625
当n=6
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
0.50000.25000.12500.06250.03130.0313
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
0.50000.25000.12500.06250.03130.01560.0156
Cholesky分解：
Cholesky.m
functionL=Cholesky（A）
L=zeros（n）;
L（1,1）=sqrt（A（1,1））;
L（i,1）=A（i,1）/L（1,1）;
forj=2:
s1=0;
j-1
s1=s1+L（j,k）^2;
L（j,j）=sqrt（A（j,j）-s1）;
fori=j+1:
s2=0;
s2=s2+L（i,k）*L（j,k）;
L（i,j）=（A（i,j）-s2）/L（j,j）;
计算Ax=b；
Test5.m
forn=10:
20
A=zeros（n,n）;
b=zeros（n,1）;
A（i,j）=1/（i+j-1）;
b（i,1）=i;
n='
方程组原始解'
x0=A\b
利用Cholesky分解的方程组的解'
L=Cholesky（A）
x=L'
\（L\b）
只列出了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分解的方程组的解
-1.2105
3.5939
-6.3219
6.5100
-3.6225
0.8405
n=
11
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
-0.0563
0.3668
-1.3972
3.2716
-4.7669
4.2094
-2.0608
0.4290
（1）House.m
functionu=House（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）;
w=a1-norm（a1,2）*e1;
H=eye（n）-2*u*u'
HA=H*A;
（3）test6.m
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）
Warning:
迭代次数太多，可能不收敛?
return;
Gauss_Seidel迭代：
Gauss_Seidel.m
function[x,n]=Gauss_Seidel（A,b,x0）
%--Gauss-Seidel迭代法解线性方程组
%--方程组系数阵A
%--方程组右端项b
%--初始值x0
%--求解要求的精确度eps
%--迭代步数控制M
%--返回求得的解x
%--返回迭代步数n
M=10000;
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f
x=G*x0+f
迭代次数太多，可能不收敛！
解方程组，test7.m
A=[5-1-3;
-124;
-3415];
b=[-2;
1;
10];
精确解'
x=A\b
迭代初始值'
x0=[0;
0;
0]
Jacobi迭代过程：
[xj,nj]=Jaccobi（A,b,x0）;
Jacobi最终迭代结果：
xj
迭代次数'
nj
Gauss-Seidel迭代过程：
[xg,ng]=Gauss_Seidel（A,b,x0）;
Gauss-Seidel最终迭代结果：
xg
ng
精确解
-0.0820
-1.8033
1.1311
迭代初始值
0
-0.4000
0.6667
err=
0.1000
-1.0333
0.4533
1.5333
...
9.6603e-006
xj=
迭代次数
升级会员 