计算方法作业第六章Word文档格式.docx
《计算方法作业第六章Word文档格式.docx》由会员分享，可在线阅读，更多相关《计算方法作业第六章Word文档格式.docx（16页珍藏版）》请在冰豆网上搜索。
N=size（A）;
n=N
（1）;
（n-1）
for
j=（i+1）:
n
if（A（i,i）==0）
disp（'
对角元素为0！
'
）;
%防止对角元素为0
return;
l
=A（j,i）;
m
=A（i,i）;
A（j,1:
n）=A（j,1:
n）-l*A（i,1:
n）/m;
%消元方程
b（j）=b（j）-l*b（i）/m;
x=SolveUpTriangle（A,b）;
%通用的求上三角系数矩阵线性方程组的函数
XA=A;
%消元后的系数矩阵
（SolveUpTriangle.m）解上三角方程组源程序：
%解上三角方程组
functionx=SolveUpTriangle（A,b）
x（n）=b（n）/A（n,n）;
fork=n-1:
1
s=0;
i=k+1:
s=s+A（k,i）*x（i）;
x（k）=（b（k）-s）/A（k,k）;
结果：
x1=[0,0,0,…..0,0,1]T
x2=[0,0,0,…..0,0,0.3820]T
x3=[0,0,0,…..0,0,0.9900]T
x4=[1,1,1,…..1,1,1]T
Cholesky：
x2=Cholesky（A,b）;
（Cholesky.m）:
Cholesky法源程序
%用Cholesky方法解线性方程组
functionx=Cholesky（A,b）;
N=size（A）;
n=N
（1）;
L（1,1）=sqrt（A（1,1））;
%L为下三角矩阵
forj=2:
L（j,1）=A（1,j）./L（1,1）
fork=2:
i=1:
k-1
s=s+L（k,i）;
L（k,k）=sqrt（A（k,k）-s）;
j=（k+1）:
sum=0;
sum=sum+L（j,i）*L（k,i）;
L（j,k）=（A（j,k）-sum）./L（k,k）;
%定义L的转置
forj=1:
k=1:
LT（j,k）=L（k,j）;
y=SolveDownTriangle（L,b）;
x=SolveUpTriangle（LT,y）;
解上三角方程组源程序：
End
解下三角方程组源程序：
functionx=SolveDownTriangle（A,b）
x
（1）=b
（1）/A（1,1）;
x1=[0,0,0,…..0,0,0.9894]T;
x2=[0,0,0,…..0,0,0.9899]T;
LTDT：
x3=GJCholesky（A,b）;
（GJCholesky.m）:
%用改进Cholesky法求解线性方程组
%记LD为U，U为下三角矩阵，LT为单位上三角矩阵
functionx=GJCholesky（A,b）;
fork=1:
U（k,1）=A（k,1）;
fori=2:
L（k,1）=A（1,k）./U（1,1）;
end
j=k:
s=s+U（j,i）.*L（j,i）;
U（j,k）=A（j,k）-s;
j=k+1:
sum=sum+U（k,i）.*L（j,i）;
L（j,k）=（A（k,j）-sum）./U（k,k）;
D（k,k）=U（k,k）;
L（k,k）=1;
LT（k,k）=1;
LT（k,j）=L（j,k）;
z（i）=y（i）./D（i,i）;
x=SolveUpTriangle（LT,z）;
n=8
x1=[0,0,0,…..0,0,0.9891]T;
x2=[0,0,0,…..0,0,0.9889]T;
n=10
x1=[0,0,0,…..0,0,0.9889]T;
x2=[0,0,0,…..0,0,0.9885]T;
n=12
x1=[0,0,0,…..0,0,0.9886]T;
n=14
x1=[0,0,0,…..0,0,0.9884]T;
x2=[0,0,0,…..0,0,0.9884]T;
n=16
x1=[0,0,0,…..0,0,0.9881]T;
x2=[0,0,0,…..0,0,0.9882]T;
n=18
x1=[0,0,0,…..0,0,0.9878]T;
x2=[0,0,0,…..0,0,0.9879]T;
n=20
x1=[0,0,0,…..0,0,0.9875]T;
x2=[0,0,0,…..0,0,0.9875]T;
追赶法：
x4=zhuigan（A,b）;
追赶法源程序
%用追赶法求解线性方程组
functionM=zhuigan（A,g）
n=length（A）;
L=eye（n）;
U=zeros（n）;
n-1
U（i,i+1）=A（i,i+1）;
U（1,1）=A（1,1）;
n
L（i,i-1）=A（i,i-1）/U（i-1,i-1）;
U（i,i）=A（i,i）-L（i,i-1）*A（i-1,i）;
Y
（1）=g
（1）;
Y（i）=g（i）-L（i,i-1）*Y（i-1）;
M（n）=Y（n）/U（n,n）;
fori=n-1:
-1:
1
M（i）=（Y（i）-A（i,i+1）*M（i+1））/U（i,i）;
（2）
Gauss消去：
n=6;
j=1:
A（i,j）=1/（i+j-1）;
s=s+A（i,j）;
b（i）=s;
[x1,XA]=GaussJordanXQ（A,b）;
Gauss消元法源程序
列主元Gauss消去：
[x2,XA]=GaussLineMainXQ（A,b）;
列主元Gauss消去法源程序
%列主元Gauss消去法解线性方程组
functionx=GaussLinemainXQ（A,b）;
index=0;
main=max（abs（A（1:
n,i）））;
%选取列主元
k=i:
if
（abs（A（k,i）==main）
index=k;
break;
%保存主元所在行的指标
temp=A（i,1:
n）;
A（i,1:
n）=A（index,1:
A（index,1:
n）=temp;
btemp=b（index）;
b（index）=b（i）;
b（i）=btemp;
%交换主元所在行
（A（i,i）==0）
对角线元素为0’）;
l=A（j.i）;
m=A（i,i）;
n）-A（i,1:
n）.*l/m;
b（j）=b（j）-b（i）.*l/m;
2．考虑矩阵
，其中
。
利用Gauss-Jordan方法求解逆矩阵，并验证逆矩阵与A的乘积是否为单位矩阵
Gauss-Jordan方法程序：
functionx=GaussJordan（A）
[n,n]=size（A）;
whiledet（A）==0
break
s=1;
P=zeros（1,n）;
whiles<
=n
max=abs（A（s,s））;
big=0;
form=s:
ifmax<
abs（A（m,s））
max=abs（A（m,s））;
q=m;
big=1;
elsecontinue
end
ifbig==1
P（s）=q;
D=A（s,:
A（s,:
）=A（q,:
A（q,:
）=D;
fori=1:
ifi~=s
forj=1:
ifj~=s
A（i,j）=A（i,j）-A（s,j）*A（i,s）/A（s,s）;
A（i,s）=-A（i,s）/A（s,s）;
fork=1:
ifk<
s
A（s,k）=A（s,k）/A（s,s）;
ifk==s
A（s,k）=1/A（s,s）;
ifk>
A（s,k）=A（s,k）*A（s,s）;
s=s+1;
fori=n:
（-1）:
ifP（i）~=0
p=P（i）;
d=A（:
i）;
A（:
i）=A（:
p）;
p）=d;
A
调用程序为：
formatrat
a=zeros（[n,n]）;
n
a（i,j）=i-j;
End
b=eye（[n,n]）+a
B=GaussJordan（b）
B*b
结果为：
B=
51/106-39/106-23/106-7/1069/10625/106
-41/10675/106-21/106-11/106-1/1069/106
-27/106-23/10687/106-15/106-11/106-7/106
-13/106-15/106-17/10687/106-21/106-23/106
1/106-7/106-15/106-23/10675/106-39/106
15/1061/106-13/106-27/106-41/10651/106
B*b
ans=
1.00000.00000-0.0000-0.0000-0.0000
-0.00001.0000-0.0000-0.00000.0000-0.0000
-0.0000-0.00001.00000.0000-0.00000.0000
-0.0000-0.0000-0.00001.0000-0.0000-0.0000
-0.0000-0.0000-0.0000-0.00001.00000.0000
0.00000.00000.00000.00000.00001.0000
是单位矩阵
升级会员 