C.2.5.0.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.f --- 常微分方程式の初期値問題を Euler 法で解く
      program rei52
*       開始時刻と終了時刻
      real a,b
      parameter (a=0.0,b=1.0)
*       初期値
      real x0
*       変数と関数の宣言
      real t,x,h,f,e
      integer N0,N1,N2,N,i
*       自然対数の底 e (=2.7182818284..)
      e=exp(1.0)
*       初期値の設定
      x0=1.0
*       どういう N について計算するか？
*       N0 から N1 まで、N2 刻みで増える N に
      write(*,*) ' FIRSTN,LASTN,STEPN='
      read(*,*) N0,N1,N2
*       計算の開始
      do N=N0,N1,N2
         h=(b-a)/N
*       開始時刻と初期値のセット
         t=a
         x=x0
*       Euler 法による計算
         do i=0,N-1
            x=x+h*f(t,x)
            t=t+h
         end do
         write(*,*) N,abs(e-x)
      end do
      end
**********************************************************************
*       微分方程式 x'=f(t,x) の右辺の関数 f の定義 -- 前と同じだから略

コンパイルして実行すると `` FIRSTN,LASTN,STEPN='' と尋ねてき ます。例えば ``100 1000 50'' と答えると $100$ から $1000$ まで $50$ 刻みで増やして行った値 $100$, $150$, $200$, $\cdots$, $900$, $950$, $1000$ を $N$ として計算して、最終的に得られた誤差を出力します。 実行結果を見てみましょう。

waltz11% reidai5-2
 FIRSTN,LASTN,STEPN=
100 1000 50
  100    1.34680271e-02
  150    9.00602341e-03
  200    6.76608086e-03
  .... (途中略) ...
  900    1.50966644e-03
  950    1.42812729e-03
  1000    1.36137009e-03
waltz11%

``mgraph''コマンドには、対数目盛によるグラフを描く機能(起 動時に -lx, -ly を指定します。ここで ``l'' は ``logarithmic''(=「対数の」) の頭文字です）があります。この場合は、例 えば次のようにしたりすると良いでしょう。

        waltz11% reidai5-2 | mgraph -lx -ly | xplot
        10,10000,100

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

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