C.2.4.0.2 例題5.1

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

* reidai5-1.f -- 微分方程式の初期値問題を Euler 法で解く
      program rei51
*       開始時刻と終了時刻
      real a,b
      parameter (a=0.0,b=1.0)
*       変数と関数の宣言
      integer N,j
      real t,x,h,x0,f
      external f
*       初期値
      x0 = 1.0
*      区間の分割数 N を入力してもらう
      write(*,*) ' N='
      read(*,*) N
*      小区間の幅
      h=(b-a)/N
*      開始時刻と初期値のセット
      t=a
      x=x0
      write(*,*) t,x
*      Euler 法による計算
      do j=0,N-1
         x=x+h*f(t,x)
         t=t+h
         write(*,*) t,x
      end do
      end
**********************************************************************
*       微分方程式 x'=f(t,x) の右辺の関数 f の定義
      real function f(t,x)
      real t,x
      f=x
      end

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