数值计算-线性方程组求解（1）-LU分解-MATLAB实现

来源：互联网发布：卓讯数据库编辑：程序博客网时间：2024/06/06 01:47

LU分解1

矩阵三角分解基本定理：设A∈Rn×n，若A的顺序主子式detAk≠0(k=1,2,⋯,n)，则存在唯一的Doolittle分解

$A = L U,$ 其中L为单位下三角矩阵，U为非奇异的上三角矩阵.
$A = L \times U = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1 l 21 l 31 ⋮ l n 1 1 l 32 ⋮ l n 2 1 ⋱ \dots ⋱ l n, n - 1 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ \times ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ u 11 u 12 u 22 u 13 u 23 u 33 \dots \dots ⋱ ⋱ u 1 n u 2 n ⋮ u n - 1, n u n n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥$
依照矩阵乘法规则展开上式可知：
U的第一行求法： $u 1 i = a 1 i (i = 1, 2, \dots, n)$
L的第一列求法： $l 1 i = a i 1 u 11 (i = 1, 2, \dots, n)$
一般地，在求出U的第k−1行和L的第k−1列后，第k步计算公式如下：
$u k i = a k i - \sum j = 1 k - 1 l k j u j i (i = k, k + 1, \dots, n)$
$l i k = a i k - \sum k - 1 j = 1 l i j u j k u k k (i = k + 1, k + 2, \dots, n; k \neq n)$

代码实现：

% FILENAME: LU.mfunction [ L, U ] = LU( A )% This function expresses a matrix A as the product of two essentially triangular matrices, one of them a permutation of a lower triangular matrix and the other an upper triangular matrix.    [N, ~] = size(A);    L = zeros(N);    U = zeros(N);    for i = 1 : N        U(1,i) = A(1, i);                      % Get the first row of U        L(i, 1) = A(i, 1) / U(1,1);            % Get the first column of L    end    for k = 2 : N                              % Get the 2~n rows of U and k-th column of L        for i = 1 : N            sum = 0;            for j = 1 : (k-1)                sum = sum+L(k, j) * U(j, i);            end            U(k, i) = A(k, i) - sum;            % U(k, i) = round(U(k, i), 4);        end        for i = 1 : N            sum = 0;            for j = 1 : (k-1)                sum = sum+L(i, j) * U(j, k);            end            L(i, k) = (A(i, k) - sum) / U(k ,k);            % L(i, k) = round(L(i, k), 4);        end    endend

在求出A的LU分解后，按照如下方式求解方程组：

% FILENAME: LU_Run.m% Your A and b go hereA=[2 1 5; 4 1 12; -2 -4 5];b=[11 27 12]';[L, U] = LU(A);[N,~] = size(A);% computing LY = bY = zeros(N, 1);Y(1) = b(1);for k = 2 : N    sum = 0;    for j = 1 : (k-1)        sum = sum + L(k, j) * Y(j);     end    Y(k) = b(k) - sum;end% computing UX = YX = zeros(N ,1);X(N) = Y(N) / U(N, N);disp(X);

测试结果展示如下：

>> AA =     2     1     5     4     1    12    -2    -4     5>> [L,U]=LU(A)L =     1     0     0     2     1     0    -1     3     1U =     2     1     5     0    -1     2     0     0     4>> L*U     % L*U == Aans =     2     1     5     4     1    12    -2    -4     5>> bb =    11    27    12>> LU_Runx =      1    -1     2>> A*X    % A*X == bans =    11    27    12

此文章中除代码外，均参考于《数值计算方法（第三版）》 ↩

阅读全文

0 0