前節の議論から、初めに作ったプログラムでまともな答が求まるためには、 でないといけないということがわかっ た。これは結構厳しい条件である。というのは、空間刻みの数を増やすと、その2 乗に比例して時間刻みの数を増やさないといけないので、全体としては空間刻 みの数の3乗に比例して計算時間がかかるということになるからである。まあ、 最近の計算機は速いので待てなくもないかもしれないが、もうちょっとなんと かならないかという気もする。
実は、なんとかなる方法というのはある。これが陰解法(implicit method)と
いうものである。拡散方程式の場合、もっとも簡単な陰解法は、空間微分に
ではなく のほうを使う、つまり、
(26) |
(27) |
線形連立方程式を解かないといけないのはもちろん面倒なわけだが、それだけ のことはある。先ほどと全く同じように とか を定義す ると、両者の関係が今度は
(28) |
陽解法と同じように行列で書いてみると、
(29) |
(30) | |||
(31) |
具体的には、 を ad, を al, を ad、右辺の値を b と表して、前進消去の部分が、
double c = al[i]/ad[i-1]; ad[i] -= c*au[i-1]; b[i] -= c*b[i-1];というだけである。後退代入も、 au のところだけを消せばいいので
b[i] = (b[i]-au[i]*b[i+1])/ad[i];というふうに一つ下の値を引くだけになる。
一応、3重対角行列を解くプログラム全体を与えておくと、
void solve_tridiagonal(double ad[], double al[], double au[], double b[], int n) { int i; for(i=1;i<n;i++){ double c = al[i]/ad[i-1]; ad[i] -= c*au[i-1]; b[i] -= c*b[i-1]; } b[n-1] /= ad[n-1]; for(i=n-2;i>=0;i--){ b[i] = (b[i]-au[i]*b[i+1])/ad[i]; } }となる。
係数の ad 等を計算して、1ステップ進める関数全体は、
void push_system(double u[], double ad[], double au[], double al[], int nx, double dx, double dt) { double delta = dt/(dx*dx); int i; for(i=0;i<nx-1;i++){ ad[i] = 1 + 2*delta; al[i] = au[i] = -delta; } solve_tridiagonal(ad, al, au, u+1, nx-1); }
である。