数値計算ライブラリLapackを使って線形連立一次方程式
\(
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
\)
をLU分解を用いて数値的に解くプログラムを載せます。

詳しくは説明しませんが、逆行列を用いて解く方法と比べ、
行列が疎である時(行列要素にゼロが多い時)に比較的計算量が減らせます。
また、逆行列を求める際に生じる桁落ちの問題を回避することが出来ます。

解法

連立一次方程式
\(
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
\)
を行列\(A\)のLU分解を利用して解きます。
ここで、未知なのはベクトル\(\mathbf{x}\)、既知なのは行列\(A\)とベクトル\(\mathbf{b}\)です。

\(A\)がLU分解されているとします。
すなわち、行列\(A\)に適当な操作を施して、

\(
\begin{align}
A=LU
\end{align}
\)
と書けているとします。ここで、行列\(L, U\)は下三角行列、上三角行列であり、
\(
\begin{align}
L=
\left(
\begin{array}{*4{>{\displaystyle}c}}
1&0&\cdots &0 \\[1em]
l_{2,1}&1&\cdots &0 \\[0.5em]
\vdots&\vdots&\ddots&\vdots \\[0.3em]
l_{s,1}&l_{s,2}&\cdots &1
\end{array}
\right) ,~~~
U=
\left(
\begin{array}{*4{>{\displaystyle}c}}
u_{1,1}&u_{1,2}&\cdots &u_{1,s} \\[1em]
0&u_{2,2}&\cdots &u_{2,s} \\[0.5em]
\vdots&\vdots&\ddots&\vdots \\[0.3em]
0 & 0 &\cdots &u_{s,s}
\end{array}
\right)
\end{align}
\)

という行列です。
もし、行列\(A\)がLU分解されていれば、

\(
\begin{align}
& A\mathbf{x}=\mathbf{b} \\
& LU\mathbf{x} = \mathbf{b} \\
& L\mathbf{w} = \mathbf{b}~~~…(*)
\end{align}
\)

と書けます。ここで、\(\mathbf{w} \equiv U\mathbf{x}\)と置きました。
まず、式(*)を解いて、\(\mathbf{w}\)を求めます。これは、下三角行列に対する問題なので、前進代入を用いて簡単に解くことが出来ます。

続いて、
\(
U\mathbf{x}=\mathbf{w}
\)
を後退代入を利用して解\(\mathbf{x}\)を求めます。

以上がLU分解を利用して連立一次方程式を解く方法です。
数値計算のルーチンは大きく2つのステップに分かれており、

1. 行列\(A\)のLU分解を求める(ルーチン dgetrf)
2. LU分解された行列\(A\)を用いて連立一次方程式を解く(ルーチン dgetrs)

というステップです。
これは、例えば問題
\(
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
\)
と
\(
\begin{equation}
A\mathbf{x}=\mathbf{b’}
\end{equation}
\)
の右辺だけが変わる複数の問題を解きたい場合なんかに便利です。この問題の場合、行列\(A\)は変わらないため、LU分解した結果を両方の問題に流用することが出来るのです。

プログラム

プログラムは以下の通りです。
いつもと同じように、ワーク配列を減らす為だけに導入したルーチンを挟んでいます。

上のプログラムをlapackと一緒にコンパイルして動かすと、端末上で

という結果が得られるかと思います。
これは、連立一次方程式
\(
\begin{align}
\left(
\begin{array}{ccc}
1& 3& 5 \\
0& 3& 1 \\
6& 2& 5
\end{array}
\right)
\left(
\begin{array}{ccc}
x \\
y \\
z
\end{array}
\right)
=
\left(
\begin{array}{ccc}
33 \\
10 \\
66
\end{array}
\right)
\end{align}
\)
を解いた結果です。
解析解(\(x=7, y=2, z=4\))は、例えばwolframを利用して求められます。
{{1,3,5},{0,3,1},{6,2,5}}.{x,y,z}=={33,10,66} –wolfram alpha

Tweet