- 1、本文档共6页,可阅读全部内容。
- 2、原创力文档(book118)网站文档一经付费(服务费),不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
- 3、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。如您付费,意味着您自己接受本站规则且自行承担风险,本站不退款、不进行额外附加服务;查看《如何避免下载的几个坑》。如果您已付费下载过本站文档,您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“版权申诉”(推荐),也可以打举报电话:400-050-0827(电话支持时间:9:00-18:30)。
查看更多
基于Matlab的n阶非奇异方阵的LU分解实现
基于Matlab 的 n阶方阵的LDU分解实现
引言
矩阵的LDU分解是“矩阵理论与方法”课程中非常重要的一部分。LDU分解在实际工程应用中也非常广泛。LDU分解可以将一个矩阵分解为一个下三角矩阵和一个对角矩阵和一个上三角矩阵的乘积。LDU分解主要应用在数值分析中,用来解线性方程、求反矩阵或计算行列式。
将系数矩阵A转变成等价两个矩阵L和U和对角矩阵的乘积 ,其中L和U分别是下三角和上三角矩阵,D为对角矩阵。当A的所有顺序主子式都不为0时,矩阵A可以分解为A=LDU。即:
Matlab 是很好的处理矩阵的工具。它的功能非常强大,包括创建矩阵,对矩阵求逆,转置等操作非常简单,使其成为图像处理,信号分析等领域常用的工具。Matlab 官方已经包括了对非奇异矩阵的LU分解函数[L,U]=lu(A),为了加深对矩阵分解的理解,本文不采用Matlab 官方的LU分解函数对矩阵A进行LDU分解,而是根据理论推导和编程实现LDU分解。
程序设计
输入合法检验
LU分解需要被分解矩阵A满足如下条件:
矩阵A为方阵
A的顺序主子式
故LU分解需先检验A为n阶方阵,然后检验A的n-1个顺序主子式全不为0,才可进行LU分解。而检验主子式可以在n-1次循环LU分解中进行,故先检验矩阵是否为方阵。
代码如下
n-1次循环LDU分解
LDU分解本质上是高斯消元法。实质上是将A通过初等行变换变成一个上三角矩阵,其变换矩阵就是一个单位下三角矩阵。从下至上地对矩阵A做初等行变换,将对角线左下方的元素变成零,然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵,这一系列单位下三角矩阵的乘积的逆就是L矩阵,它也是一个单位下三角矩阵。LDU分解主要分为两步:1根据高斯消元法对消元,消元矩阵为;2计算以产生下一步迭代的。
根据构造
高斯消元,使第i+1列从第i+2行至n行都为0。构造消元矩阵。首先判断是否为0,为0则无法继续分解,退出;否则继续。
代码如下
计算
计算并储存覆盖,并计算下次循环的主子式。
代码如下
计算并返回和为结果
通过n-1次循环累乘得到L矩阵,并且根据最后得到的矩阵分解出D和U矩阵。代码如下
验证结果
可LDU分解的矩阵验证
构造可以LDU分解的矩阵A
调用MyLU函数,分解A矩阵。
计算L*D*U验证是否为原矩阵A(方法下同)
输入矩阵A是否为空验证
输入矩阵A中没有元素,输入不合法。
输入矩阵A为方阵验证
输入矩阵A维度不合法。
顺序主子式出现0的错误验证
因为计算时出现,导致顺序主子式。该情况为不合法输入。
至此,函数功能和合法性检查全部验证完毕。
心得体会
通过“矩阵理论与方法”的理论指导和Matlab编程的实践经验,我基本掌握了矩阵分解中的LDU分解的推导过程和算法步骤。熟练掌握LDU分解,对今后研究LU分解、Doolittle分解、Crout分解、QR分解等矩阵分解方法的实现有非常大的帮助。对LDU分解的推导过程进行步骤分解和归纳,我将N维矩阵的LDU分解归纳总结为n-1次循环,每步循环进行n-i次元素除法(计算高斯消元系数),2次N维矩阵乘法(计算和)和1次元素乘法(计算顺序主子式的值)。按乘法计算,即时间复杂度为。因为计算过程中需要辅助计算矩阵和,即空间复杂度为。该LDU分解设计还不够快速,占用空间相对较多,是以后改进的方向。
附录:程序源码
function [L,D,U] = MyLU(A)
%check validity
if(isempty(A))%check A if is empty
error(A is empty!);
end
[N,D] = size(A);
if(N ~= D)%check A if is square
error(A is not a square!);
end
%LDU decomposing
L = eye(N);
det_k = A(1,1);
for step_n = 1:N-1
if (det_k == 0)% check if Sequence principal minor appears 0
error(Sequence principal minor is 0);
end
Li = eye(N);
Li_inv = eye(N);
for step_row = step_n+1:N
mod=A(step_row,step_n)/A(step_n,step_n);
Li(step_row,step_n) = mod;
Li_inv(step_row,step_n) = -mod;
end
A = Li_inv * A;
det_k = det_k *
文档评论(0)