A..6.0.1 [l u]=lu(a) とすると

置換行列を省いて [l u]=lu(a) とすると、 l には $PL$ が、u には $U$ が代入される。

[L U P]=lu(a);
[L0 U0]=lu(a)

とすると

L0 =
    1.0000         0         0         0
    0.5000    1.0000    1.0000         0
    0.3333    1.0000         0         0
    0.2500    0.9000   -0.6000    1.0000

U0 =
    1.0000    0.5000    0.3333    0.2500
         0    0.0833    0.0889    0.0833
         0         0   -0.0056   -0.0083
         0         0         0    0.0004

本当に $\texttt{L0}=PL$, $\texttt{U0}=U$ かチェックする。

norm(U0-U)
norm(L0-P*L)

どちらも0になる。

[L U p]=lu(a,'vector')

とすると

L =
    1.0000         0         0         0
    0.3333    1.0000         0         0
    0.5000    1.0000    1.0000         0
    0.2500    0.9000   -0.6000    1.0000

U =
    1.0000    0.5000    0.3333    0.2500
         0    0.0833    0.0889    0.0833
         0         0   -0.0056   -0.0083
         0         0         0    0.0004

p =
     1     3     2     4

p は置換行列 $P$ の “情報を持った” ベクトルである。

• ベクトル $c$ に対して、$Pc$ を求めるには P*c とする代わりに、 c(p) とすれば良い。
• 行列 $C$ に対して、$P C$ を求めるには P*C とする代わりに、 C(p,:) とすれば良い。

$PA=LU$ となっているとき、

$$Ax=b\quad\Iff\quad PAx=Pb\quad\Iff\quad LUx=Pb$$

であるから、$x$ を求めるには、 U\(L\(P*b))、 あるいは U\(L\b(p)) とする。

解が分かりやすい連立１次方程式を用意しよう。

x=(1:4)';
b=a*x;

解は、もちろん

a\b

とすれば求まるが、[L U P]=lu(a); で求めた L, U, P を使うには

U\(L\(P*b))

とすれば良く、[L U p]=lu(a,'vector'); で求めた L, U, p を使うには

U\(L\b(p))

とすれば良い。

>> x=(1:4)'

x =
     1
     2
     3
     4

>> b=a*x;
>> a\b

ans =
    1.0000
    2.0000
    3.0000
    4.0000

>> U\(L\(P*b))

ans =
    1.0000
    2.0000
    3.0000
    4.0000

>> U\(L\b(p))

ans =
    1.0000
    2.0000
    3.0000
    4.0000