直接三角分解LU法

silence123

%直接三角分解法LU;
! |9 }6 ?  d8 I! S# W/ q: Aformat rational
disp('请输入矩阵的阶数');
: {- p* E0 F% i3 k/ c" e) B2 CN=input('N=');
( V. W5 k1 d3 b# D' _for i=1:N
9 R' j3 I  S9 {' T( q    for j=1:N
# k( g, s+ t5 X* T  `. |) F8 f% I        A(i,j)=input('A(i,j)=');
    end
end
L=eye(N);
( V! |. o& T: j# h; Vfor i=1:N
    U(1,i)=A(1,i);
    L(i,1)=A(i,1)/U(1,1);
end
. O' b1 W: _. p- m, K' r+ C# W0 yfor r=2:N
" a" q6 A' d! s3 A( N0 `! ~/ f# Q    for i=r:N
        U(r,i)=A(r,i);
& c3 ~6 [& c/ X: o6 J        for k=1:r-1
            U(r,i)=U(r,i)-L(r,k)*U(k,i);
1 c: [7 f/ d2 W0 b$ e4 S' b        end
        L(i,r)=A(i,r);
/ b, ~$ l5 O3 }2 v        for k=1:r-1
9 ], C- K$ V  I% ~- s7 i            L(i,r)=L(i,r)-L(i,k)*U(k,r);
' S' l1 u$ O7 u/ n3 r        end
7 X; b# [4 B. }) `        L(i,r)=L(i,r)/U(r,r);
    end
end
- M0 H0 r9 ^! p+ Q1 t& Q8 A# ?1 B! Mdisp('A=');
disp(A);
disp('L=');
" U) s$ Q2 P1 i$ |! ?/ Fdisp(L);
6 n* Y) ~( I9 L7 xdisp('U=')
disp(U);

0.9清1.8清2.7清