Subsections

• 7.1 3重対角行列を解く

7 陰解法

前節の議論から、初めに作ったプログラムでまともな答が求まるためには、 ${\Delta t}/{\Delta x^2} < {1}/{2}$ でないといけないということがわかっ た。これは結構厳しい条件である。というのは、空間刻みの数を増やすと、その2 乗に比例して時間刻みの数を増やさないといけないので、全体としては空間刻 みの数の3乗に比例して計算時間がかかるということになるからである。まあ、 最近の計算機は速いので待てなくもないかもしれないが、もうちょっとなんと かならないかという気もする。

実は、なんとかなる方法というのはある。これが陰解法(implicit method)と いうものである。拡散方程式の場合、もっとも簡単な陰解法は、空間微分に $u_j$ ではなく $u_{j+1}$ のほうを使う、つまり、

$$\frac{\partial^2 u}{\partial x^2} \sim \frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}$$ (26)

というものを使うというだけである。時間微分は同じものですます。 差分公式としてまとめて書くと、

$$\frac{u_{i,j+1}-u_{i,j}}{\Delta t} = \frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}$$ (27)

となる。 最初の方法では、 $u_{i,j+1}$ は左辺に単独で現れ、右辺は時刻 $j$ の項だ けだったので、$u_{i,j+1}$を直接計算できた。このように、直接計算できる 方法のことを陽解法 explicit method という。これに対し、上に書いた方法 では時刻 $j+1$の項が右辺に複数現れているので、これは全体として ${\boldmath {u}}_{j+1}$ についての線形連立方程式になっている。このように、新しい時 刻での解が代数方程式を解かないと求められない方法のことを陰解法という。

線形連立方程式を解かないといけないのはもちろん面倒なわけだが、それだけ のことはある。先ほどと全く同じように $\delta$ とか $\lambda$ を定義す ると、両者の関係が今度は

$$\lambda = \frac{1+\alpha \delta}{1+2\delta}$$ (28)

となる。したがって、 $\vert \alpha \vert < 2$ なる任意の $\alpha$ と $\delta > 0$ な る任意の $\delta$ について、 $\vert\lambda\vert<1$ が満たされているのである。 つまり、どんなに $\Delta t$ を大きくとっても、決して不安定にならない。 この性質を無条件安定という。この場合には、安定性ではなくて計算精度の観 点から $\Delta t$ を決められることになり、だいぶましである。

7.1 3重対角行列を解く

陽解法と同じように行列で書いてみると、

$$A{\boldmath {u}}_{j+1} = {\boldmath {u}}_j$$ (29)

で、

$$a_{ii} = 1+2 \delta$$ (30)

$$a_{i,i-1} = a_{i-1,i} = -\delta$$ (31)

ということになる。$A$ は $\delta$ の符号が違うだけで前と同じ三重対角行 列になっている。行列（以下、線形連立方程式のことを単に行列ということ がある）を解くのは例えばガウスの消去法では計算量が一般には次元の3乗に なるので嬉しくないが、三重対角の場合には計算量が次元数に比例する。前進 消去で消さないといけない項が常に1つだけだからである。

具体的には、$a_{ii}$ を ad, $a_{ii-1}$ を al, $a_{ii+1}$ を 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);
}

である。