3.2.5.1 例題 5.2

例題5.1 と同じ初期値問題で、色々な分割数 $N$ に対して 問題を時、$t=1$ での誤差の大きさ $\vert x(1)-x_N\vert=\vert e-x_N\vert$ について調べよ。

こういう場合は、色々な $N$ に対して一斉に $\vert e-x_N\vert$ を計算するプログ ラムを作るのがいいでしょう。

reidai5-2.c

/*
 * reidai5-2.c --- 常微分方程式の初期値問題を Euler 法で解く
 * http://nalab.mind.meiji.ac.jp/~mk/program/ode_prog/reidai5-2.c
 */

#include <stdio.h>
#include <math.h>

int main(void)
{
  /* 開始時刻と終了時刻 */
  double a = 0.0, b = 1.0;
  /* 初期値 */
  double x0;
  /* 変数と関数の宣言 */
  double t,x,h,f(double,double),e;
  int N0,N1,r,N,i;
  /* 自然対数の底 e (=2.7182818284..) */
  e = exp(1.0);
  /* 初期値の設定 */
  x0 = 1.0;
  /* どういう N について計算するか？ 
   * N0 から N1 まで、r をかけて増える N に */
  printf("# FIRSTN,LASTN,r=");
  scanf("%d%d%d", &N0,&N1,&r);
  /* 計算の開始 */
  for (N = N0; N<= N1; N *= r) {
    h = (b-a)/N;
    /* 開始時刻と初期値のセット */
    t = a;
    x = x0;
    /* Euler 法による計算 */
    for (i = 0; i < N; i++) {
      x += h * f(t,x);
      t += h;
    }
    printf("%d %e\n", N, fabs(e-x));
  }
  return 0;
}

/* 微分方程式 x'=f(t,x) の右辺の関数 f の定義 */
double f(double t, double x)
{
  return x;
}

コンパイルして実行すると “ FIRSTN,LASTN,r=” と尋ねてきます。 例えば “10 1280 2” と答えると、 $10$ から始めて $2$ 倍ずつしていって $1280$ を超えないまでの値 (つまり $10$, $20$, $40$, $80$, $160$, $320$, $640$, $1280$) を $N$ として計算して、最終的に得られた誤差を出力します。 実行結果を見てみましょう。

コンパイル＆実行

oyabun% gcc -o reidai5-2 reidai5-2.c -lm
oyabun% ./reidai5-2
# FIRSTN,LASTN,r=10 1000 2
10 1.245394e-01
20 6.498412e-02
40 3.321799e-02
80 1.679689e-02
160 8.446252e-03
320 4.235185e-03
640 2.120621e-03
1280 1.061069e-03
oyabun%

出力結果を保存して作ったファイル “reidai5-2.data” の内容を gnuplot でグラフにするには、 以下のようにします。 両側対数目盛によるグラフを描く機能を用いています。

oyabun% ./reidai5-2 > reidai5-2.data
10 1280 2
oyabun%

これで計算結果を reidai5-2.data に記録できた。 この内容を gnuplot でプロットする。

oyabun% gnuplot

        G N U P L O T
        Unix version 3.7
        (以下色々なメッセージ…省略)

gnuplot> set logscale xy
gnuplot> plot "reidai5-2.data" with linespoints

  (以下印刷用のデータ作り)
gnuplot> set term postscript eps color
Terminal type set to 'postscript'
Options are 'eps noenhanced color dashed defaultplex "Helvetica-Ryumin"
	14'
gnuplot> set output "reidai5-2.eps"
gnuplot> replot
gnuplot> quit
oyabun%

図 1: Euler法の誤差 (横軸: 分割数 $N$, 縦軸: 誤差の絶対値)

(2021/4/1追記: ここでは EPS (Encapsulated PostScript) 形式にしていますが、 現在では set term png や set term pdf とする方が良いかも。)

このグラフから、誤差 $=O(N^{-1})\quad(N\to+\infty)$ であることが読み とれます。実はこれは Euler 法の持つ一般的な性質です。

実は Euler 法は収束があまり速くないので、実際には特殊な場合を除いて 使われていません。そこで…