直接三角分解LU法

silence123

%直接三角分解法LU;
format rational
/ F6 v; J" H$ Z: _# jdisp('请输入矩阵的阶数');
9 m, L  w7 |4 Q# z( i0 s  v7 mN=input('N=');
for i=1:N
9 X6 A$ b( I- u' x& T6 M* d+ d: j    for j=1:N
, ~" w* u* t4 `. I0 w        A(i,j)=input('A(i,j)=');
    end
end
L=eye(N);
7 ?' I  V# E! V/ v0 Lfor i=1:N
6 f/ J$ k) d  W7 z! ~- V" c    U(1,i)=A(1,i);
+ E6 e) z4 ~+ f    L(i,1)=A(i,1)/U(1,1);
end
for r=2:N
% |3 a$ x! v0 k. w/ Y' n1 m: |    for i=r:N
4 U3 A3 {: S* `7 n        U(r,i)=A(r,i);
7 V, J! p! P# r* v5 c9 U& l! C1 L$ i' K        for k=1:r-1
            U(r,i)=U(r,i)-L(r,k)*U(k,i);
        end
+ x/ v. R. i$ q* ~; z/ v, C; `        L(i,r)=A(i,r);
, t; {% H7 n  o, W$ t: t        for k=1:r-1
            L(i,r)=L(i,r)-L(i,k)*U(k,r);
5 V4 z: P- L2 i+ z, }2 N+ ?) }& V        end
        L(i,r)=L(i,r)/U(r,r);
0 v6 V9 X) l/ |% u* D! ]3 R2 O) a    end
; m' S/ l; y/ Q- R# cend
disp('A=');
disp(A);
disp('L=');
disp(L);
disp('U=')
4 C) l) D+ \- e: bdisp(U);

0.9清1.8清2.7清