3.2.4.2 例題5.1

Euler 法を用いて、例1 の初期値問題を解くプロ グラムを作って収束を調べよ。

reidai5-1.c

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

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

このプログラムをコンパイルして¹² 実行すると、分割数 $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$ に近付いていくはずです。