4 波動方程式の差分近似 -- I 君のプログラム

$u_{tt}$, $u_{xx}$ を共に 2 階中心差分で近似してできた差分方程式を 「波動方程式に対する差分法」 で説明してあるが、他にも色々な方法がある。

Friedrichs の差分法と呼ばれるものが有名である (藤沼祐一「波動方程式に対する差分法」 を見よ) が、 今日のI君が説明したようにして、 時間について1階の方程式にして解くのも自然である。

ただ、今日のは、 これまで卒研で学んだことの多くを「おいといて」という感じの内容でした。 ツッコミを箇条書きにしておく。

(a)

$0.00625$ のようなマジック・ナンバーをプログラムに書いてはいけない。 これは X ($=160$) の逆数であろうから、 dt=1.0/X; のようにすべきである。 こうしておけば X の値を書き換えても正常に動作する。

(b)

X や t という変数名は誤解されかねないので、 不適当であろう。 特に x や t は別の意味で使うことになる可能性が高い (偏微分方程式の記述に現れるのだから)。

(c)

double [][]u = new double[X+2][T+2]; のように +2 とするのは、慎重なようで決してそうではない。 正しく理解していれば double [][]u = new double[X+1][T+1]; で必要十分と言うことが分かるはず。

(d)

そもそも double [][]u = new double[X+2][T+2]; のように 2次元配列にするのは拙い。 時間方向には 3 ステップ分だけあれば大丈夫。

ちょっと遊んでみる。

   1 /* fib1.c --- 配列を使う */
   2 
   3 #include <stdio.h>
   4 
   5 #define N (20)
   6 
   7 int main(void)
   8 {
   9   int i;
  10   int a[N+1];
  11   a[0] = 1;
  12   a[1] = 1;
  13   for (i = 2; i <= N; i++)
  14     a[i] = a[i - 1] + a[i - 2];
  15   for (i = 0; i <= N; i++)
  16     printf("a[%2d]=%5d\n", i, a[i]);
  17   return 0;
  18 }

   1 /* fib2.c --- 隣接3項の漸化式だから3つメモリーがあれば十分 */
   2 
   3 #include <stdio.h>
   4 
   5 #define N (20)
   6 
   7 int main(void)
   8 {
   9   int i;
  10   int am2, am1, a;
  11   am2 = 1;
  12   am1 = 1;
  13   printf("a[%2d]=%5d\n", 0, am2);
  14   printf("a[%2d]=%5d\n", 1, am1);
  15   for (i = 2; i <= N; i++) {
  16     a = am1 + am2;
  17     printf("a[%2d]=%5d\n", i, a);
  18     am2 = am1;
  19     am1 = a;
  20   }
  21 
  22   return 0;
  23 }

   1 /* fib3.c --- ベクトル列になったときはコピーが面倒なので、添字でローテート
   2               まあ、リングバッファーみたいなもんです。 */
   3 
   4 #include <stdio.h>
   5 
   6 #define N (20)
   7 
   8 int main(void)
   9 {
  10   int i;
  11   int a[3];
  12   a[0] = 1;
  13   a[1] = 1;
  14   printf("a[%2d]=%5d\n", 0, a[0]);
  15   printf("a[%2d]=%5d\n", 1, a[1]);
  16   for (i = 2; i <= N; i++) {
  17     a[i%3] = a[(i-1)%3] + a[(i-2)%3];
  18     printf("a[%2d]=%5d\n", i, a[i%3]);
  19   }
  20 
  21   return 0;
  22 }

   1 /* fib4.c --- 割り算 % をやめて */
   2 
   3 #include <stdio.h>
   4 
   5 #define N (20)
   6 
   7 int main(void)
   8 {
   9   int i,I,Im1,Im2;
  10   int a[3];
  11   a[0] = 1;
  12   a[1] = 1;
  13   printf("a[%2d]=%5d\n", 0, a[0]);
  14   printf("a[%2d]=%5d\n", 1, a[1]);
  15   Im1 = 1; Im2 = 0;
  16   for (i = 2; i <= N; i++) {
  17     I = Im1 + 1; if (I == 3) I = 0;
  18     a[I] = a[Im1] + a[Im2];
  19     printf("a[%2d]=%5d\n", i, a[I]);
  20     Im2 = Im1;
  21     Im1 = I;
  22   }
  23 
  24   return 0;
  25 }