7. 練習問題答え・2

もともと、

のような形を初期条件としている。これが一定速度で右に動いて行ってもらいたいと考えているわけだが、実際には

となってしまい、計算が安定してくれない。そこで差分の方法をやりなおして、

$$T_i^{n+1}=T_i^{n}-u\frac{\Delta t}{\Delta x}(T_{i}^n-T_{i-1}^n)$$

としてみる。計算を100ステップごとに結果を書きだすように変更して、さらに時間時間を10から500に変更してみると、

#include <stdio.h>
#include <math.h>
#define u 1.0
#define dt 0.1
#define xmax 200
#define imax 200
#define tmax 500

int main(void){
  double dx=(double)xmax/imax;
  double x;
  double T[xmax+1],cfl;
  int t;
  cfl=u*dt/dx;
  int i;
  for (i=0; i<=imax; i++){
    if (i>=10&&i<50)
      T[i]=1.0;
    else
      T[i]=0.0;
  }
  for (i=1; i<imax; i++)
    printf("%2.2e,",T[i]);
  printf("\n");   
  for (t=1; t<=tmax; t++){
    for (i=1; i<imax; i++)
      T[i]=T[i]-cfl*(T[i]-T[i-1]);
    if (t%100==0){
      for (i=1; i<imax; i++)
        printf("%2.2e,",T[i]);
      printf("\n");
    }
  }  
} 

すると、

となって、かなり安定した計算ができていることに気が付く。しかし本来は青色の角ばった形状であったものが、500ステップ後には丸まった形へと変化してしまっている。これは特に角の部分の値が大きく変化するところが、次第になまっていくことに原因があることは容易に想像できる。空間の離散化をもっと細かくすると、少しは改善される部分があろうが、解決しないことも想像できるだろう。

いずれにせよ、風上側の差分を使って解く方が安定していることがわかった。これを風上差分と呼ぶ。しかし速度$u$が逆方向を向くと、差分表記から変更しないと発散してしまう可能性があるのは、都合が悪い。特に複雑な現象で、右にも左にも移動することがあるケースは、うまく解けないことになる。

差分の方法をやりなおして、

$$T_i^{n+1}=T_i^{n}-u\frac{\Delta t}{2\Delta x}(T_{i+1}^n-T_{i-1})$$

としてみるとどうだろうか?結果は

となり、前進差分よりはマシだが、振動しながら進むという結果が得られた。

以上のように、少し差分形式を変えるだけで、得られる解が大きく変わってしまうことがある。これは差分法だけの問題ではないが、収束計算を行わない場合には特に問題となる場合がある。なお移流項を含むこの問題は、多項式で近似するなどして、特に変化の大きい部分について取り扱う手法(CIP法など)を用いて解決する場合が多い。