消元法实验报告13.docx
《消元法实验报告13.docx》由会员分享，可在线阅读，更多相关《消元法实验报告13.docx（10页珍藏版）》请在冰点文库上搜索。
消元法实验报告13
西京学院数学软件实验任务书
课程名称
数学软件实验
班级
数0901
学号
0912020114
姓名
王斌
实验课题
线性方程组直接三角分解法（Doolittle分解，Grout分解），平方根法（Cholesky分解，LDLT分解）
实验目的
熟悉线性方程组直接三角分解法（Doolittle分解，Grout分解），平方根法（Cholesky分解，LDLT分解）
实验要求
运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成
实验内容
线性方程组直接三角分解法（Doolittle分解，Grout分解）
线性方程组平方根法（Cholesky分解，LDLT分解）
成绩
教师
一．线性方程组的直接三角分解法（Doolittle分解）
设Ａ为非奇异矩阵，且有分解式Ａ=LU，其中L为单位下三角矩阵，U为上三角矩阵，则称此分解为Doolittle分解。
此分解的目的是将Ax=b分解为两步，首先由Ly=b解出y，再由Ux=y解出x。
以下为实现此分解的程序：
function[x,l,u]=Doolittle（A,b）
clc
clearall
formatshort
n=input（'请输入矩阵的阶数:
'）;
A=zeros（n,n）;
fori=1:
n
forj=1:
n
A（i,j）=input（'请输入矩阵中的元素：
'）;
end
end
A
ifdet（A）~=0
fori=1:
n
b（i）=input（'请输入b中的元素：
'）;
end
b=b'
u=zeros（n,n）;
l=eye（n,n）;
u（1,:
）=A（1,:
）;
fori=2:
n
form=1
l（i,m）=A（i,m）/u（1,1）;
fork=2:
n
forj=k:
n
u（k,j）=A（k,j）-sum（l（k,1:
k-1）*u（1:
k-1,j））;
fork=2:
n-1
l（k+1:
n,k）=（A（k+1:
n,k）-sum（l（k+1:
n,1:
k-1）*u（1:
k-1,k）））/u（k,k）;
end
end
end
end
u
l
end
y=zeros（n,1）;
y
（1）=b
（1）;
fork=2:
n
y（k）=b（k）-l（k,1:
k-1）*y（1:
k-1）;
end
x=zeros（n,1）;
x（n）=y（n）/u（n,n）;
fork=n-1:
-1:
1
x（k）=（y（k）-u（k,k+1:
n）*x（k+1:
n））/u（k,k）;
end
end
end
此程序运行时，首先输入一个系数矩阵A，和列向量b如果系数矩阵A是非奇异矩阵，则矩阵最后可求得A分解后的矩阵L和矩阵U，以及列向量x。
以下是运行结果：
二．Cholesky分解法。
设A是n（n>=2）阶实对称矩阵，L是非奇异的下三角矩阵，则称
为矩阵A的Cholesky分解。
其计算步骤为将A分解后，首先由Ly=b,求得y，再由
x=y,求得解向量x。
以下为实现此分解的Matlab程序:
function[x]=pingfg（A,b）
clc
clearall
formatshort
n=input（'请输入矩阵的阶数:
'）;
fori=1:
n
forj=1:
n
A（i,j）=input（'请输入矩阵中的元素'）;
end
end
A
fori=1:
n
b（i）=input（'请输入b中的元素：
'）;
end
b=b'
[n,n]=size（A）;
L=zeros（n,n）;
L（1,1）=sqrt（A（1,1））;
fork=2:
n
L（k,1）=A（k,1）/L（1,1）;
end
fork=2:
n-1
L（k,k）=sqrt（A（k,k）-sum（L（k,1:
k-1）.^2））;
fori=k+1:
n
L（i,k）=（A（i,k）-sum（L（i,1:
k-1）.*L（k,1:
k-1）））/L（k,k）
end
end
L（n,n）=sqrt（A（n,n）-sum（L（n,1:
n-1）.^2））;
%
y=zeros（n,1）;
fork=1:
n
j=1:
k-1;
y（k）=（b（k）-L（k,j）*y（j））/L（k,k）;
end
x=zeros（n,1）;
U=L';
fork=n:
-1:
1
j=k+1:
n;
x（k）=（y（k）-U（k,j）*x（j））/U（k,k）;
end
此程序运行时，首先输入一个系数矩阵A（A的阶数>=2）且A为实对称正定矩阵，和列向量b矩阵最后可求得A分解后的非奇异的下三角矩阵L，以及列向量x。
以下是运行结果：
三．改进平方根分解法（LDLT分解）
因为用Cholesky分解后求x时需要开方，于是可以改进平法根分解法，使其分解时不需要开方。
即将A分解为
从而由Ly=b解得y,再由
=y解得列向量x。
从而解的此线性方程组的解。
以下是实现此过程的Matlab程序：
functionimprovecholesky
clc
clearall
formatshort
n=input（'请输入矩阵的阶数:
'）;
fori=1:
n
forj=1:
n
A（i,j）=input（'请输入矩阵中的元素:
'）;
end
end
A
fori=1:
n
b（i）=input（'请输入b中的元素:
'）;
end
b=b'
L=zeros（n,n）;
D=diag（n,0）;
S=L*D;
fori=1:
n
L（i,i）=1;
end
fori=1:
n
forj=1:
n
if（eig（A）<=0）|（A（i,j）~=A（j,i））
disp（'wrong'）;break;end
end
end
D（1,1）=A（1,1）;
fori=2:
n
forj=1:
i-1
S（i,j）=A（i,j）-sum（S（i,1:
j-1）*L（j,1:
j-1）'）;
L（i,1:
i-1）=S（i,1:
i-1）/D（1:
i-1,1:
i-1）;
end
D（i,i）=A（i,i）-sum（S（i,1:
i-1）*L（i,1:
i-1）'）;
end
L
D
y=zeros（n,1）;
x=zeros（n,1）;
fori=1:
n
y（i）=（b（i）-sum（L（i,1:
i-1）*D（1:
i-1,1:
i-1）*y（1:
i-1）））/D（i,i）;
end
fori=n:
-1:
1
x（i）=y（i）-sum（L（i+1:
n,i）'*x（i+1:
n））;
end
x
此程序运行时，首先输入一个系数矩阵A（A的阶数>=2）且A为实对称正定矩阵，和列向量b矩阵最后可求得A分解后的非奇异的下三角矩阵L和对角矩阵D，以及列向量x。
以下是运行结果：
升级会员 