1.LU分解
在线性代数中已经证明,如果方阵是非奇异的,即的行列式不为0,LU分解总是存在的,它的形式为:
例如,对于四阶矩阵来说,可以分解为:
寻找规律:
- 第一行没有依赖因素,因此可以直接得出,第一列仅仅依赖第一行,所以在得到第一行之后,第一列也可以计算求得.
- 同理,第二主行也是这样,可以先计算主行,在计算主列,也就是按照下图中的,先实线,再虚线的顺序计算.
- 仔细观察,实际上只要先计算出U矩阵主对角线元素
,至于先计算行还是列,无关紧要,因为除了对应行,列的未知lu元素,不在依赖其他参变量.
- 每个原矩阵中的元素
仅仅使用一次,计算出对应的
后,其信息已经扩散到整个计算结构中,可以舍弃。后续计算不再依赖,对计算机程序计算来说,好处是为算法的内存优化,提供的抓手。
按照上面的发现的规律,尝试推导计算公式如下:
按照U的行线性组合
得到:
另外,按照L的列的线性组合
所以得到:
对于内层矩阵,由于时,
,所以,为求
时候的值,可以计算得到:
所以:
由于时,
,所以,为求
时候的值,可以计算得到:
所以, 得到:
根据上述两个公式,在octave上编码实现为:
% find the LU factorization of the matrix
% input:
% a: the matrix need to be factorize
% n: the number of the low or column in the matrix
% output:
% no output
function LU(a,n)
m=eye(n); %genernate the unit matrix.
for j = 1 : n-1
if abs(a(j,j))<eps;
error('zero pivot encountered'); % when the zero pivot happens,end the process
end
for i = j+1 : n
mult = a(i,j)/a(j,j);
m(i,j) = mult;
for k = j:n
a(i,k) = a(i,k) - mult*a(j,k);
end
end
end
disp(' L='); disp(m);
disp(' U='); disp(a);
disp(' LU='); disp(m*a); % to check if the result is right
运行结果为:
上面的程序需要矩阵的每阶顺序主子式都不为0,否则程序会报错.例如:
这个问题有时间在修复.