A..5 不思議な常識: 微分方程式を解くときに逆行列は使わない

前小節までに、LU 分解があれば、 逆行列と同じ計算量で連立1次方程式が解けることを説明したが、 実は微分方程式を解く場合には、現れる大抵の行列が疎行列であることから、 LU 分解を使うのが断然有利である。

前回

>> n=10
>> a=2*eye(n-1,n-1)-diag(ones(n-2,1),1)-diag(ones(n-2,1),-1)

のような行列を例に出しましたが、 毎回こう入力するのが面倒なので、関数を作りましょう、

>> edit Lap1d

とすると、編集ウィンドウが現れます。 そこに次のように入力し、ファイルを保存します。

Lap1d.m

function a=Lap1d(n)
  a=2*eye(n-1,n-1)-diag(ones(n-2,1),1)-diag(ones(n-2,1),-1);

>> a=Lap1d(10)
>> a=Lap1d(20)

前回説明した \ 演算子で連立1次方程式を解いてみます。

\ 演算子を使って解く

>> n=1000
>> tic; Lap1d(n) \ rand(n-1,1); toc
経過時間は 0.069753 秒です。
  
>> n=2000
>> tic; Lap1d(n) \ rand(n-1,1); toc
経過時間は 0.381871 秒です。

>> n=4000
>> tic; Lap1d(n) \ rand(n-1,1); toc
経過時間は 2.136665 秒です。

未知数の個数 n はこれくらいが無難です。 n がある程度大きくなると、メモリーが不足して、

??? メモリが足りません。オプションについては、HELP MEMORY と入力してください。

エラー ==> Lap1d at 2
  a=2*eye(n-1,n-1)-diag(ones(n-2,1),1)-diag(ones(n-2,1),-1);

のようにギブアップされたり、スワップを始めて、 計算自体は行うものの、非常に遅くなります。 Mathematica と違って、MATLAB は簡単に停止できません (Control-C を打ってしばらく待つか、MATLAB を強制的に終了することになります)。

逆行列を用いて連立1次方程式を解いてみましょう。

逆行列 inv() を使って解く

>> n=1000
>> tic; inv(Lap1d(n))*rand(n-1,1); toc
経過時間は 0.300695 秒です。

>> n=2000
>> tic; inv(Lap1d(n))*rand(n-1,1); toc
経過時間は 2.283513 秒です。

>> n=4000
>> tic; inv(Lap1d(n))*rand(n-1,1); toc
経過時間は 16.784775 秒です。

最後の $n=4000$ の場合、8倍程度余計に時間がかかっていることに注目。

一方、LU分解を用いて解いて見ると、

LU分解 lu() を使って解く

>> n=1000
>> tic; [L u p]=lu(lap1d(n)); u\(L\(p*rand(n-1,1))); toc
経過時間は 0.165507 秒です。
>> clear
>> n=2000
>> tic; [L u p]=lu(lap1d(n)); u\(L\(p*rand(n-1,1))); toc
経過時間は 0.483060 秒です。
>> clear
>> n=4000
>> tic; [L u p]=lu(lap1d(n)); u\(L\(p*rand(n-1,1))); toc
経過時間は 3.649599 秒です。

逆行列を使って解くより速く、 \ 演算子を使って解くのと同じオーダーの時間で 解けています。