物理・天文で出てくる偏微分方程式は大抵は 2 階である。で、普通は空間3次 元と時間1次元の4次元空間で定義されるわけだが、この講義ではいくつかの特 殊な場合を除いて空間1次元時間1次元の2次元で考える。 これは、そうしておかないとプログラムを書くのも、また計算機のほうも大変 になるからである。

放物型方程式ってのは何か？というのは皆さん知っているはずなので（だよね？） 厳密な定義は省くが、要するに以下の形に書ける（移流）拡散方程式のことで ある。

ここで x, t はそれぞれ空間、時間を表す変数であり、 u は方程式が 記述する量である。拡散方程式として見ればなにかの濃度ということになるし、 熱伝導の方程式としてみれば温度なりエンタルピーなりということになる。 K は一次の係数で、これは普通の拡散方程式や熱伝導方程式では 0 である。 これが 0 でないのはどういう場合かというのは後でちょっと考える。 D が 普通の拡散係数ということになる。今は D が空間・時間に依存しない場合 を考える。

依存しない場合はもちろん変数分離で解けるので、数値計算することに意味が あるのはD が空間・時間に依存する場合だが、まあ、とりあえずは答がわか る場合を計算してみることにしよう。

放物型方程式の場合には、初期条件と境界条件を与えないと解が定まらない。 以下では、

偏微分方程式の数値計算の方法は基本的には有限要素法と差分法に大別される。 あ、あとスペクトラル法というのもあるが、これはちょっとおいておく。 有限要素法の考え方はだいぶややこしいので、まず差分法を考える。

差分法とは、基本的には微分方程式に出てくる空間微分や時間微分の項を、差 分で置き換えるものである。差分というのはそもそも何かというのを説明 しないといけない。

なお、差分のとり方は一意的ではないことに注意して欲しい。例えば、 でもいいし、 といったものも考えられる。さらに、もっ と沢山の点を使うことも原理的には可能である。

さて、拡散方程式なので2階微分がいる。これはどう作ればいいかというわけ だが、これはむしろ一般に高階微分を差分で表す、つまり数値微分の方法を説 明しておくほうが簡単であろう。

の 点での関数値 をつかって、点 での 階導関数の近似値を求める考え方は以下 のようなものである。

1. 各点で を満たす 次補間多項式 を作る。
2. を 回微分する。
3. 微分した結果の式に を代入する。

初期条件は のところ にだけ値がある 関数的なものということにする。

なお、面倒なので ということにする。時間の単位が と思えば同 じことである。

     1://     program parabolic1
     2:
     3:#define _USE_MATH_DEFINES
     4:#include <iostream>
     5:#include <unistd.h>
     6:using namespace std;
     7:#include <cmath>
     8:#include <gr.h>
     9:
    10:void initparam(double  u[],
    11:               int & nx,
    12:               double &dx,
    13:               double & dt,
    14:               int & niter,
    15:               int &gmode)
    16:{
    17:    cerr << "Enter nx, tmax, dt:";
    18:    double tmax;
    19:    cin >> nx >> tmax >> dt;
    20:    dx = 1.0/nx;
    21:    niter = (tmax+dt/2)/dt;
    22:    cerr << "Enter graphics mode 0: no G\n"
    23:         << "                    1: animation\n"
    24:         << "                    2: no-erase)\n";
    25:    cin >> gmode;
    26:    cout << "nx = " << nx<< " dt = " << dt << " niter = "<< niter<< "\n";
    27:    cout << "gmode = " <<  gmode << "\n";
    28:    int i;
    29:    for (i=0;i<=nx; i++) u[i] = 0;
    30:    u[(nx+1)/2] = nx;
    31:}
    32:
    33:void push_system(double u[],
    34:                 double unew[],
    35:                 int nx,
    36:                 double dx,
    37:                 double dt)
    38:{
    39:      double lambda = dt/(dx*dx);
    40:      int i;
    41:      for (i=1; i<nx; i++){
    42:          unew[i] = u[i] + lambda*(u[i-1]-2*u[i]+u[i+1]);
    43:      }
    44:      for(i=1; i<nx; i++){
    45:         u[i] = unew[i];
    46:      }
    47:}
    48:
    49:      
    50:void show_graphics(double u[],
    51:                   int nx,
    52:                   double dx,
    53:                   int gmode,
    54:                   double xmin,
    55:                   double xmax,
    56:                   double ymin,
    57:   double ymax,
    58:                   int first)
    59:{
    60:    if (gmode == 0) return;
    61:    double x[nx+1];
    62:    if (gmode == 1) {
    63:gr_clearws();
    64:    }
    65:    int i;
    66:    for(i=0;i<nx+1;i++){
    67:x[i]=(i*dx);
    68:    }
    69:    gr_setwindow(xmin,xmax,ymin,ymax);
    70:    gr_axes(0.25*(xmax-xmin), 0.25*(ymax-ymin), 0, 0, 1, 1, -0.01);
    71:    gr_polyline(nx+1,x, u);
    72:    if (gmode == 1) usleep(30000);
    73:    gr_updatews();
    74:}
    75:      
    76:int main()
    77:{
    78:#define NMAX 10001
    79:    static double u[NMAX], unew[NMAX];
    80:    double dx, dt;
    81:    int nx, niter, gmode, iter;
    82:    initparam(u, nx, dx, dt, niter, gmode);
    83:    show_graphics(u, nx, dx, gmode, 0.0, 1.0, 0.0, 10.0, 1);
    84:    for(iter = 1; iter <=niter; iter++){
    85:         push_system(u, unew, nx, dx, dt);
    86:         show_graphics(u, nx, dx, gmode,0.0, 1.0, 0.0, 10.0, 0);
    87:    }
    88:    cerr << "Enter any key to finish:";
    89:    char c;
    90:    cin>>c;
    91:    return 0;
    92:
    93:}

コンパイルには

       g++  parabolic1.cpp -o parabolic1   -std=c++11 -I $GRDIR/include\
         -L  $GRDIR/lib -Wl,-rpath,$GRDIR/lib -lGR

なお、環境変数 GKS_WSTYPE を mp4 にすることで動画ファイルを作成できる。 また、 画面に出す時には環境変数 GKS_DOUBLE_BUF を設定する(値は1とかtrueとか適 当に)と、画面のチラツキがなくなる。

具体的には、拡散方程式なので本来は次第に滑らかになっていくべきなのに、 どんどんガタガタになってしまう。しかも、これは であれば や をどんなに小さくしても起こる。

まず、なぜそういう振動が起きるのかを考え、それから振動を起こさない方法 があるのかどうかを議論していこう。

さて、まず、なぜ振動が起きるかということだが、これは原理的には差分化さ れた方程式の固有値を求めればわかる。

この行列は実対称行列なので固有値分解ができて、適当なユニタリ行列を持っ てきて

このように分解できるので、 とすれば最初の差分近似式は

不安定性が起きないための条件は、全ての固有値 の絶対値が を超えないということである。

断るまでもないと思うが、これは元の偏微分を変数分離して空間方向の関数に ついての常微分方程式を導くのとほとんど同じ操作になっている。これを についての差分方程式と見て解を求めることを考える。

線形差分方程式なので、解は の形である。これを代入すると

真面目にやるには境界条件を満たす固有ベクトルをもとめないといけないが、 面倒なので無限遠境界の場合を考える。この時、 が実数のものはどち らかで無限大に発散するのでよろしくないので、複素数の場合を考える。この 時は になっているので無限遠でも発散しない。

なお、有限の固定境界の場合も結局境界条件を満たすような解を作るためには が複素数でないといけないことがすぐにわかる。このため、 だけを考えればよい。

上の議論からわかるように、不安定条件に関係する の値は実際には だけで、それ以外の からはもっと緩い条件しかでな い。 は に対応するので、これは結局 と で値が逆転するようなパターンである。そのようなパターン（モー ド）が 1 ステップ計算するとどうなるかをもう一度もとの差分式に戻って書き直してみ ると、結局

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

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

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

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

具体的には、 を ad, を al, を ad、右辺の値を b と表して、前進消去の部分が、

        double c = al[i]/ad[i-1];
        ad[i] -= c*au[i-1];
        b[i] -=  c*b[i-1];

        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);
 }

さて、前節の簡単な陰解法は無条件安定で大変結構なのであるが、時間刻みを 大きくとれるとなると、今度は計算精度のほうが気になってくる。空間方向の 微分は2次精度で、テイラー展開の2次の項までが正しく入っているのに対し、 時間微分は1次精度にしかなっていないからである。

なお、上の1次の陽解法のことを1次の前進差分、陰解法のことを1次の後退差 分ということがある。

時間方向の精度をあげる一つの方法は、これまでに見てきた2つの方法を混ぜ て使うことである。つまり、

1. メインプログラムを付け加えて、1次の陰解法のプログラムを完成し、 をいろいろ変えて安定に計算できることを確認せよ。
2. クランク・ニコルソン法が無条件安定であることを証明せよ。
3. クランク・ニコルソン法のプログラムを作り、上と同様に安定であるこ とを確認せよ。
4. クランク・ニコルソン法の誤差が , のそれぞれ何乗になっているか、 , をいろいろ変えてプログラムを走ら せて結果を厳密解と比べることで調べよ。厳密解と比べるには、初期条件がデ ルタ関数では具合が良くないのでどういうものを使うべきか考えてみること。
5. クランク・ニコルソン法が時間方向2次精度であることを証明せよ。
6. 前進差分、クランク・ニコルソン法、後退差分は、 時刻 と での空間微分の混ぜ方を変えているだけである。これらすべてを統一的に計算 できるプログラムを作成し、正しく動くことを確認せよ。
7. さらに、周期境界の場合も扱えるようにプログラムを拡張せよ。3重対 角行列の処理の部分をどう変えればいいか？