2.2.1 Euler法のプログラム

Euler法のプログラムをCで書くには、大きく分けて二つのやり方がある。

1. 配列を使う。

   #define MAXN (1000)
   
     double x[MAXN+1];
   ...
     x[0] = x0;
     t = a;
     for (j = 0; j < N; j++) {
       x[j+1] = x[j] + h * f(t,x[j]);
       t += h;
     }
   

2. 配列を使わないで済ませる。

      double x; // 配列でなくて良い。
   ...
      x = x0;
      t = a;
      for (j = 0; j < N; j++) {
        x += h * f(t,x);
        t += h;           // t = a + (j+1)*h の方が良いかも。
      }
   

test-euler.c

/* test-euler.c -- 微分方程式の初期値問題を Euler 法で解く */

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

int main(void)
{
  /* 開始時刻と終了時刻 */
  double a = 0.0, b = 1.0;
  /* 変数と関数の宣言 */
  int N, j;
  double t,x,h,x0,f(double, double);
  /* 初期値 */
  x0 = 1.0;
  /* 区間の分割数 N を入力してもらう */
  printf("N="); scanf("%d", &N);
  /* 小区間の幅 */
  h = (b-a) / N;
  /* 開始時刻と初期値のセット */
  t=a;
  x=x0;
  printf("t=%f, x=%f\n", t, x);
  /* Euler 法による計算 */
  for (j = 0; j < N; j++) {
    x += h * f(t,x);
    t += h;
    printf("t=%f, x=%f\n", t, x);
  }

  return 0;
}

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

コンパイル＆実行

cc test-euler.c
./a.out

分割数 $N$ を尋ねてきますので、 色々な値を入力して試してみて下さい。 各時刻 $t_j$ における $x_j$ の値 ( $j=0,1,\cdots,N$ ) を画面に出力します。

確認用にいくつかの $N$ の値に対する場合の、$x(1)=x_N$ の値を書いてお きます。$N=10$ の場合 $x_N=2.59374261$ , $N=100$ の場合 $x_N=2.70481372$ , $N=1000$ の場合 $x_N=2.71692038$ , $N=10000$ の場合 $x_N=2.71814346$ , $\cdots$ .

分割数 $N$ が大きくなるほど、真の値 $e=2.7182818284590452\cdots$ に近付いていくはずです。