コンテンツ

オイラー法

今日の演習

1階常微分方程式の数値解法（オイラー法）

今日の目標

1階常微分方程式を解く方法を知る

常微分方程式の数値解法

次のような常微分方程式の初期値問題を、オイラー法を用いて数値的に解き、gnuplot を用いて結果をグラフにしてもらいます。

次の初期値問題を考えます。

(1) $\begin{equation*} \left\{ \begin{array}{ccc} \displaystyle\frac{dy}{dt} &=& f(t, y) \\ y(0) &=& a \end{array} \right. \end{equation*}$

求めたいのは、 $t>0$ の範囲の関数 $y(t)$ であり、関数 $f(t, y)$ および定数 $a$ は与えられているとします。また、 $f(t, y)$ の値は $t, y$ が与えられればいつでも計算できるものとします。（プログラム中では関数として定義すると便利！）

前回までのgnuplotでのグラフの描画と同様に、コンピューターでは数値を離散的に（とびとびに）しか扱うことができません。そこで、式(1)をコンピューターで扱えるようにする必要があります。式(1)をコンピュータで解く方法には色々とありますが、この演習では代表的なオイラー法（とルンゲ・クッタ法）を用います。

テイラー展開を使った導出

$y(t)$ について、時刻 $t$ から $h>0$ だけ進んだ時刻 $t+h$ における $y(t+h)$ について考えます。テイラー展開を使うと

(2) $\begin{equation*}y(t+h)=\sum_{k=0}^{\infty}\frac{h^k}{k!}\frac{d^ky(t)}{dt^k}\end{equation*}$

と $h$ の多項式で表現できます。ここで、 $h\2$ 以上に関する項は、 $h$ が十分小さいことを考えると、ランダウの記号 $\mathcal{O}$ を用いて次のように表せます。

(3) $\begin{equation*}y(t+h)=y(t)+h\frac{dy(t)}{dt}+\mathcal{O}(h^2)\end{equation*}$

$\mathcal{O}(h^2)$ は残りの項がおおよそ $h^2$ に比例していることを示しています。今我々が考えている式(1)を使うと、 $\frac{dY(t)}{dt}=f(t,Y(t))$ ですので、

(4) $\begin{equation*}y(t+h)=y(t)+hf(t,y(t)+\mathcal{O}(h^2)\end{equation*}$

となります。ここで、 $\mathcal{O}(h^2)$ を無視した式を考えたいのですが、これまで考えていた $y(t)$ とは厳密には $\mathcal{O}(h^2)$ だけ異なるので、等号で表記するのは憚られます。そこで、数値計算解を $Y(t)$ として表現すると、

(5) $\begin{equation*}Y(t+h)=Y(t)+hf(t,Y(t))\end{equation*}$

となります。この式がオイラー法の根幹となる次ステップの予測式になります。ある時刻において $y(t)=Y(t)$ として、次のステップにおける誤差（局所誤差）を考えます。式(4) $,$ 式(5)は

(6) $\begin{equation*}y(t+h)-Y(t+h)=\mathcal{O}(h^2)\end{equation*}$

となります。すなわち、真の解 $y(t+h)$ と数値解 $Y(t+h)$ の間にはおおよそ $h^2$ に比例した局所誤差が生じることを意味しています。局所誤差については、このあとのプログラムで実際に確認していきます。

比較：差分商を用いた導出

高校までの微分の定義を使って導出することもできます。ただし、この場合は誤差に関する議論はできませんので、注意が必要です。 $f(x)$ を $x$ の関数としたとき、 $f(x)$ の $x=a$ における微分係数は

(7) $\begin{equation*}\lim_{h\rightarrow 0} \frac{1}{h}\left\{ f(a+h) - f(a) \right\}\end{equation*}$

で定義されます。ここで、定義中の極限操作を取り払い、 $h$ を有限にとどめた

(8) $\begin{equation*}D = \frac{1}{h}\left\{ f(a+h) - f(a) \right\}\end{equation*}$

を考えると、 $h$ を十分小さな値とすれば、式(8)は式(7)の近似となっていると考えられます。この $D$ のように、関数のいくつかの点における値の差を用いてその関数の微分係数を近似することを差分近似といい、式(8)の右辺のような量を差分商といいます。いまの場合は１階微分係数を近似する１階差分商です。

では、このような差分商を用いて式(1)を離散化してみましょう。まず、式(1)の第一式の微分係数を上記の差分商で置き換えてみます。

(9) $\begin{equation*}\frac{1}{h}\left\{ y(t+h) - y(t) \right\} = f(t, y(t))\end{equation*}$

式(9)と式(1)の第一式は異なる方程式なので、式(1)の第一式の解 y(t) は一般に式(9)を満たしません。そこで、混乱を防ぐために式(9)の $y(t)$ を $Y(t)$ と書き換えましょう。すると、

(10) $\begin{equation*}\frac{1}{h}\left\{ Y(t+h) - Y(t) \right\} = f(t, Y(t))\end{equation*}$

式(10)のように差分を含む方程式を差分方程式といいます。この式を整理すると、式(5)と同じ式を得ることができます。

帰納的解法

　次に、式(1)の $y$ を $Y$ に置き換えた初期条件

(11) $\begin{equation*}Y(0) = a\end{equation*}$

の下で式(5)を解くことを考えます。式(5)に $t = 0$ を代入すると、

(12) $\begin{equation*}Y(h) = Y(0) + hf(0, Y(0))=a+hf(0, a)\end{equation*}$

となり、 $Y(0)$ から直ちに $Y(h)$ が求まります。同様にして式(5)を繰り返し用いると、

(13) $\begin{equation*} \left\{ \begin{array}{rll} Y(2h) &=& Y(h)+hf(h, Y(h)), \\ Y(3h) &=& Y(2h)+hf(2h, Y(2h)), \\ Y(4h) &=& Y(3h)+hf(3h, Y(3h)), \\ &\vdots& \\ Y((j+1)h) &=& Y(jh)+hf(jh, Y(jh)), \\ &\vdots& \\ \end{array} \right. \end{equation*}$

というように、 $j=1, 2, 3, \dots$ とすると帰納的に $Y((j+1)h)$ を $Y(jh)$ から計算できることがわかります。

$h$ をいったん決めると、 $t=jh$ 以外の時刻の $Y$ の値は式(5)からは求めることができません。このように、差分方程式を解くと従属変数はとびとびの時刻で値が定まります。そのようなとびとびの時刻を格子点と呼びます。 $Y$ は格子点でのみ意味があるので、そのことを明示するために、 $Y_j=Y(jh)$ および $t_j = jh$ と数列で表現とすると、式(5)と初期条件は、

(14) $\begin{equation*} \left\{ \begin{array}{rll} Y_{j+1} &=&Y_{j} + hf(t_j, Y_j), \\ Y_0 &=& a \end{array} \right. \end{equation*}$

となり、結局式(1)の常微分方程式の初期値問題が、 $Y_j$ に関する漸化式の問題式(14)に置き換えられました。この方法をオイラー法と呼びます。

オイラー法のアルゴリズム

以下では

(15) $\begin{eqnarray*} &\dot y= -y\\ &y(0)=1 \end{eqnarray*}$

について考えていきます。基本的には前回の課題3、4でやったような漸化式を解くアルゴリズムが流用できます。例えば、前回の課題３の解答例として

#include <stdio.h>
  
double f(int, double); // 漸化式の右辺の定義
 
int main() {
    int j, N = 20;
    double Y = 1.0;
    
    for (j = 0; j < N; j++) {
        printf("%d %lf\n", j, Y); // Yの出力
        Y = f(j, Y); // Y_{j+1} = f(Y_j)
    }
    return 0;
}

/* 漸化式の右辺 */
double f(int n, double Y) {
    return 1 + 1.0/Y;
}

を用いてみましょう。

【注意】後々の説明のため、以下の点を変更しています。

- - 独立変数： $i$ → $j$
  - $i$ に対する従属変数： $x_i$ → $Y_j$
  - 漸化式の右辺を関数 f で表現

漸化式のプログラムとオイラー法のプログラムで変更・追加しなければならないことを列挙します。下の箇条書き先頭の【　】内に書かれている行番号は、この後に書かれているオイラー法のサンプルプログラム（euler.c）の行番号に対応しています。

【7, 8行目】解を求める時刻の上限を Tとし、適当な自然数 N (格子点の数。時間刻み数ともいう) を定める。
【9行目】時間刻み幅 $h=T/(N-1)$ を定義
【8, 13行目】時刻 $t_j=j\times h$ を定義（グラフを描くときの横軸になる。int型からdouble型になるので、諸々チェック！）
【15, 22行目】漸化式を式(5)に変更（main関数の中と関数 f の中）
【10, 15, 16行目】（今後連立方程式を扱うことも念頭において）次の時刻での Y の値である Y_newを定義。
【16行目】for文の最後に、次の時刻のためにYをY_newで更新

次のサンプルプログラムは

$\dot{y}&=-y\\y(0)&=1$

すなわち式(1)における $f(t,y)=-y$ を例にプログラムを作成しています。

オイラー法のサンプルプログラム：euler.c

#include <stdio.h>

double f(double, double); // 差分方程式の右辺の定義
#define a (1.0)

int main() {
    int j, N = 100000;
    double t, T = 10.0;
    double h = T/(N-1);
    double Y = a, Y_new;
    
    for (j = 0; j < N; j++) {
        t = j * h; // 現在の時刻 
        printf("%lf %lf\n", t, Y); // 現在のYの出力
        Y_new = Y + h * f(t, Y); // Y_{j+1} = Y_j + h * f(t, Y_j)
        Y = Y_new; // 次の時刻のためにYを更新
    }
    return 0;
}

/* 差分方程式の右辺 */
double f(double t, double Y) {
    return -Y;
}

今後は、サンプルプログラムを陽には書かずに、以下のようなアルゴリズムで書くことが増えます。お決まりである#include <stdio.h>や変数の定義などが書かれておらず、またほとんど自然言語です。対応する制御構文や関数など、必要に応じて各自でアルゴリズムを元に、プログラミングしていく練習をしていきましょう。

例 euler.cのアルゴリズム

時刻の上限 T（下限は 0 ）, 時間刻み数(格子点数) N を設定
時間刻み幅 h = T/(N-1) を設定
（ここから計算開始）
初期値を与える Y = 1
各格子点での値を j = 0,...,N-1 の順に計算 (for 文）
    t = j*h // 現在の時刻 
    （一定間隔毎に）t と Y の値を画面に表示（printf関数）
    Y_new = Y + h*f(t, Y)
    次の時刻の計算のため、Y を新しく求まった値（Y_new）で置き換え
以上を繰り返す

オイラー法　あれこれ

#defineの利点と欠点・注意点

　第９回で#defineの使い方を習いました。#defineを使うと、プログラム全体にわたる定数を１ヶ所で定義することができます。例えば上にあるeuler.cの場合、NやTおよびhなどはdefine文を使って次のように定義しても良いかもしれません。

euler2.c

#include <stdio.h>
#define N (100000)
#define T (10.0)
#define h (T/(N-1))

double f(double, double); // 差分方程式の右辺の定義

int main() {
    int j;
    double t;
    double Y = 1.0, Y_new;
    
    for (j = 0; j < N; j++) {
        t = j * h; // 現在の時刻 
        printf("%lf %lf\n", t, Y); // 現在のYの出力 
        Y_new = Y + h * f(t, Y); // Y_{j+1} = Y_j + h * f(t, Y_j)
        Y = Y_new; // 次のステップのためにYを更新
    }
    return 0;
}

/* 差分方程式の右辺 */
double f(double t, double Y) {
    return -Y;
}

$\dot{y}=ay, y(0)=1$ のように、微分方程式に定数係数（ここでは $a$ ）が含まれる場合も、#defineを使うと全ての関数に跨って定義できるため有用です。一方、一度定義してしまうと、それ以降で値を変更することができなくなります。すなわち、例えば「実行の度にシミュレーション終了時間Tを入力する」場合には、#defineではなく通常の変数と同様にmain関数内で定義し、scanfで値を指定することが必要です。また#defineでは、プログラムに含まれるトークン全てを変更してしまうため、上の例euler2.cの１６行目

Y_new = Y + h * f(t, Y);

では実際には

Y_new = Y + ((10.0)/((100000)-1)) * f(t, Y);

と置き換わっています。つまり、((10.0)/((100000)-1))と言う計算が毎回（無駄に）実行されていることになります。このような無駄な計算を省きたければ、やはりmain関数内で通常の変数として定義したほうが良いでしょう。 ※ 今回のような場合は大した無駄には思えないかもしれませんが、実際の研究ではもっと複雑で長いプログラムを実行する必要があり、このような無駄の積み重ねが命取りになりかねません。無駄計算をしていないか、注意しましょう。

時間と刻み数・刻み幅

上の例のeuler.cやeuler2.cでは、まず刻み数 N を定義し、次に時間の上限 T を定義し、最後にそれらから刻み幅 h を定義しています。一方、このようにまず刻み数 N の値を陽に定義し、刻み幅 h を間接的に決まるような方法は、科学研究にあまり向いていません。その理由として計算の精度があげられます。精度についての詳しい話は次回する予定ですが、ともかく、刻み幅 h を陽に決めるほうが一般的です。つまり、「時刻の上限Tと刻み幅 h を先に決め、刻み数 N を間接的に決める」というやり方や、「時刻の上限Tと刻み幅 h を先に決め、t = j*h > Tとなったら終了する」というやり方のほうが良いでしょう。実際にどのようなプログラムになるか、各自で考えてみてください（課題2(4)）。

一部分だけ出力

h が小さくなる（N が大きくなる）と、全ての j について t と Y_new を出力すると、結構なファイルサイズになり、PCに負担がかかります。そのため、出力する部分を以下のようなif文で囲むことで、○回に1回出力するようにできます (a % b = [a を b で割った時の余り])。

if (j % 100 == 0) {
    printf("%lf %lf\n", t, Y);
}

この例では作成されるデータファイルの容量は1/100になります。これは、第９回で見たように、グラフはある程度の折れ線数であれば十分なめらかに見えるので、グラフを描く為のデータとして巨大になりすぎることを防ぐ為です。ここで100という数値は適当に与えた物であって、グラフが十分なめらかに見えかつデータ数が小さくなるような適当な値を探す必要があります。（通常の連続なグラフであれば、グラフの点が200点くらいあれば十分滑らかに見えるでしょう。）

課題１

常微分方程式の初期値問題

$\dot{y} = -y + 3\exp(-t), y(0) = 1$ を解け。

オイラー法を用いて時刻 0≦ t ≦20 の範囲での解を求めるプログラムを作成せよ。時刻 0≦ t ≦20 の範囲を20000等分割して(N = 20000)、100ステップ毎に printf関数により、t および y の値を表示し、それを gnuplot でグラフ化せよ。また手計算にて解を求め、数値計算結果が正しいか確認せよ。

課題２

次の常微分方程式の初期値問題をオイラー法を用いて解いて、 $t$ と $y$ の関係を gnuplot でグラフ化せよ。

(1) $\dot{y} = 10 -10t, y(0)=1$ (時刻 0≦ t ≦2 の範囲を20000等分割、100ステップ毎に値を表示)

(2) $\dot{y} = -y + \sin(t), y(0) = 1$ (時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)

(3) $\ds \dot{y} =2\frac{y}{t+1}, y(0) = 1$ (時刻 0≦ t ≦10 の範囲を10000等分割、100ステップ毎に値を表示)

(4) $\dot{y} = y-2y^3, y(0) = 0.1$ もしくは　 $y(0) = -0.1$ 　(時刻 0≦ t ≦20 の範囲で時間刻み幅h=0.001、100ステップ毎に値を表示)
※ (4)はhを先に定め、TとhからNを計算する方式にしてください。

今日学んだ事

オイラー法による常微分方程式に数値解法

Masashi FUJII

計算数学演習第１１回

オイラー法

今日の演習

今日の目標

常微分方程式の数値解法

テイラー展開を使った導出

比較：差分商を用いた導出

帰納的解法

オイラー法のアルゴリズム

オイラー法　あれこれ

#defineの利点と欠点・注意点

時間と刻み数・刻み幅

一部分だけ出力

課題１

課題２

今日学んだ事

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University

オイラー法

今日の演習

今日の目標

常微分方程式の数値解法

テイラー展開を使った導出

比較：差分商を用いた導出

帰納的解法

オイラー法のアルゴリズム

オイラー法 あれこれ

#defineの利点と欠点・注意点

時間と刻み数・刻み幅

一部分だけ出力

課題１

課題２

今日学んだ事

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University

オイラー法　あれこれ