2.5.1 Euler法

$[t_0,T]$ に属する $t$ について解 $x(t)$ を “知る” ため、

$$t_0<t_1<t_2<\cdots<t_N=T$$

となる $t_n$ ( $n=1,2,\cdots,N$) を取り、$x(t_n)$ を近似的に求める (離散変数法という)ことを目標とする。

そのため

$$x_{n+1}=x_n+\Delta t_n f(x_n,t_n),\quad \Delta t_n:=t_{n+1}-t_n$$ (2.11)

により $\left\{x_n\right\}_{0\le n\le N}$ を定め、 $x_n$ を $x(t_n)$ の近似値として採用するのが、 (前進) Euler法 (the (forward) Euler method) である。 これは

$$x'(t_n)=\lim_{\Delta t_n\to0}\frac{x(t_{n}+\Delta t_n)-x(t_n)}{\Delta t_n}$$

に基づく前進差分近似

$$x'(t_n)\kinji\frac{x(t_{n+1})-x(t_n)}{\Delta t_n}$$ (2.12)

を用いたものである。 つまり、微分方程式 $x'(t)=f(x(t),t)$ を $t=t_n$ で考え、 左辺の微分係数を (F.4) の右辺で置き換え、 $j=n,n+1$ について $x(t_j)$ を $x_j$ と書き換えると、 (C.18) が導出できる。

C言語プログラマー向けの説明

「常微分方程式の初期値問題を解くプログラムの書き方 §2 Euler法, Runge-Kutta法入門 -- 1次元の問題」 特に「2.3.1 Euler法のCプログラム例」 にもう少し詳しい説明がある。

そこからC言語で書かれたサンプル・プログラムを抜き出す。

euler.c

/*
 * euler.c (Euler method for Malthusian model)
 */

#include <stdio.h>

double a = 1.0;

int main(void)
{
  int i, N;
  double t, x, dt;
  double f(double, double), x0;
  double Tmax;
  // 初期値
  x0 = 1.0;
  // 最終時刻
  Tmax = 1.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  dt = Tmax / N;
  // 初期値
  t = 0.0;
  x = x0;
  printf("# t x\n");
  printf("%f %f\n", t, x);
  // Euler 法
  for (i = 0; i < N; i++) {
    x += dt * f(x, t);
    t = (i + 1) * dt;
    printf("%f %f\n", t, x);
  }
  return 0;
}

double f(double x, double t)
{
  return a * x;
}

コンパイルして実行すると $N$ の値をどうするか尋ねてくるので、 $10$, $100$ のような “大きい” 数を入力すると、$t_j$ と $x_j$ ( $j=1,2,\cdots,N$) を表示する。

コンパイル&実行

% cc euler.c
% ./a.out
# N: 100
# t x
0.000000 1.000000
中略
1.000000 2.704814
%

初期値問題の解は $x(t)=e^t$ であるから、 最後の数は $x(1)=e=2.7182818284590\cdots$ に近いはずである。 $N$ を大きくすると、$e$ に近くなることが分かるはずである (試してみよう)。

データファイル作成&解曲線描画

% ./a.out > euler.data
100
% gnuplot

ここで gnuplot の起動メッセージが表示される(省略する)。

gnuplot> plot "euler.data" with lp

これで解のグラフ(解曲線)が表示される。以下は、グラフを画像ファイルとして保存する手順を示す。

gnuplot> set term png
gnuplot> set output "euler.png"
gnuplot> replot
gnuplot> quit
%

これで画像ファイル euler.png ができる。

必修課題 2.5.1   上のプログラムを実行して、解のグラフ (解曲線) を描いてみよ。

Python プログラマー向けの説明     上の euler.c を Python で書き換えると、以下のようになる。

euler.py

# euler.py
# Cプログラムの Python バージョン
# euler.c(http://nalab.mind.meiji.ac.jp/~mk/labo/text/ode-workbook/node10.html)
# Euler法という訳で実用性はない。参考プログラム。

def f(x, t):
    return a * x

a=1
x0=1
Tmax=1
str=input('N='); N=int(str);
dt=Tmax/N

t=0.0; x=x0
print('%f %f' % (t,x))
for i in range(N):
    x += dt * f(x, t)
    t = (i+1) * dt
    print('%f %f' % (t,x))

グラフを描くように書き直すのは簡単である。matplotlib を使ってみよう。

euler2.py

# euler2.py
import numpy as np
import matplotlib.pyplot as plt

def f(x, t):
    return a * x

a=1
x0=1
Tmax=1
N=100; dt=Tmax/N;
t=np.linspace(0.0,Tmax,N+1)
x=np.empty(N+1)

x[0]=x0
#print('%f %f' % (t[0],x[0]))
for i in range(N):
    x[i+1] = x[i] + dt * f(x[i], t[i])
    #print('%f %f' % (t[i+1],x[i+1]))

plt.plot(t,x)
plt.show()

このプログラムでは、 $\{t_n\}_{n=0}^N$, $\{x_n\}_{n=0}^N$ を配列 (numpy の ndarray) t, x に記憶させて、plot している。

タイトルや、座標軸が何を表すか等、 有益な情報を表示するようにプログラムを書き直すのは読者に任せる。

(前進)Euler法は精度が低く、数値的安定性についても “平凡な性能” なので、 実際に利用されることは少ない。 多くのテキストでは、次項で説明する Runge-Kutta 法の利用が勧められている。