消元法实验报告13.docx
- 文档编号:23153366
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:10
- 大小:84.07KB
消元法实验报告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分解)
设A为非奇异矩阵,且有分解式A=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。
以下是运行结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 消元法 实验 报告 13