直接三角分解LU法

silence123

%直接三角分解法LU;
format rational
' J# F% A% U! g) Wdisp('请输入矩阵的阶数');
: b& c# X  p8 o, D( `& A1 k! EN=input('N=');
3 g% u; t  I/ ~for i=1:N
8 Z  P. l9 u7 j0 a( n+ [    for j=1:N
+ ~/ o) h5 _$ n, N3 G/ r        A(i,j)=input('A(i,j)=');
    end
( T( u) v3 b- `3 r6 f. n9 {9 fend
L=eye(N);
: B. a5 C7 M% W- w9 Dfor i=1:N
* C3 k& h8 [0 N$ {' }! N    U(1,i)=A(1,i);
    L(i,1)=A(i,1)/U(1,1);
end
for r=2:N
    for i=r:N
6 F0 q: m: ?, ?9 V& s; Z" Y        U(r,i)=A(r,i);
4 D' t% j/ H, W! _# a! F; b        for k=1:r-1
8 U2 o* k3 O; O7 Y3 P            U(r,i)=U(r,i)-L(r,k)*U(k,i);
        end
        L(i,r)=A(i,r);
! Q# y3 |% n! i2 |        for k=1:r-1
4 g: i% z* j' N  o            L(i,r)=L(i,r)-L(i,k)*U(k,r);
- M" a" y3 A* ^) k/ H( J        end
  j$ F8 f9 C0 e9 M        L(i,r)=L(i,r)/U(r,r);
    end
end
disp('A=');
disp(A);
disp('L=');
' o! x! f: T: Z# Idisp(L);
0 v) Q4 M+ k4 Odisp('U=')
4 E0 G* F7 P. Y* G- p% X6 t$ jdisp(U);

0.9清1.8清2.7清