MATLABでオイラー法を実装してみよう。

こんにちは、arcadia13です。

14回目のMATLAB記事となる今回は「オイラー法」を扱います。

オイラー法とは数値計算を使って常微分方程式の解ーここでいう解とは特殊解のことーを求めるアルゴリズムの一つです。

数学的に理解しやすく、またプログラムでの実装も簡単です。

本記事は微分方程式の知識を前提にして進めます

1. 前準備

オイラー法の理屈を説明する前に、いくつか準備をします。

1.1 常微分方程式のおさらい

まずは、オイラー法に深く関わる「一階常微分方程式」のおさらいです。

これは、ある変数 $x$ とその関数 $y=y(x)$、および $y$ の一次導関数 $\frac{dy}{dx}$ によって構成される常微分方程式のことです。

その一般表現を以下に示します。

\frac{dy}{dx}=f(x,y)

$\frac{dy}{dx}$ が、$x,y$ の関数 $f$ とイコールで結ばれています。関数という表現をしていますが、要するに $x,y$ を使った何かしらの数式ですよ、という意味です。

言うまでもありませんが、この常微分方程式の解は $y(x)$ です。

これが今回の主役となります。

1.2 数値計算の重要性

二つ目として、微分方程式には数値計算がとても重要だという話をします。

11回の時、微分方程式には手計算で解けるものがすごく少ないという問題があると書きました。

そして、それに対抗するべく人間が生み出したのが「コンピュータの数値計算を使って微分方程式を解く」という方法です。これを、微分方程式を数値的に解くといいます。

数値計算で求めた解は近似解と呼ばれ、また本物の解のことを真の解と呼びます¹。

今日では、真の解よりも近似解を利用する場合が多いです。理由はいくつかありますが、一番大きな理由は解けない微分方程式がとても多いという事実です。

また、近年では過去と比較にならないほど高スペックなコンピュータが誕生していて、そのマシンパワーを使って膨大な量の数値計算を、高速に行えるようになりました。

所謂デスクトップPCはもちろん、ノートPCもその例外ではありません²。

したがって、コンピュータ、というよりも数値計算に頼る戦法は、今や手計算以上に有効な手段なのです。

1.3 前進差分

三つ目として前進差分というものを説明します。

これは微分をコンピュータ上で実装する、つまり数値計算を使って実装する方法です。

例えば連続な一変数関数 $y=y(x)$ の微分は

\frac{dy(x)}{dx}=\lim_{h \rightarrow 0} \frac{y(x+ h) - y(x)}{h}

と極限を用いて定義されます。

ここで、極限操作「$h$ をだんだん $0$ に限りなく近づけていく」をコンピュータでは厳密に再現できません。

ではどうするのか、という問いに対する答え（の一つ）が前進差分です。

微分の極限操作を、コンピュータでは「$h$ を小さな正の定数にする」ことで代用します。

こうすることで、微分を

\frac{dy(x)}{dx}\approx\frac{y(x+ h) - y(x)}{h}

と「$y(x+ h) - y(x)$ を $h$ で割る」という割り算に置き換えてしまうのです。

これを前進差分といいます。

なお、前進差分以外にも「中心差分」や「後退差分」という方法もありますが、ここでは割愛します。

1.4 常微分方程式へ前進差分を適用

前準備はこれで最後です。

オイラー法の主役となるのが、以下に示す一階常微分方程式

\frac{dy}{dx}=f(x,y)

でしたね。

右辺は単純な数式なので問題ありません。しかし問題は左辺。関数 $y=y(x)$ の一次導関数が入っています。つまり微分が含まれているので、このままではコンピュータで実現することができません。

そこで左辺を前進差分で近似します。

つまり、以下のように

\frac{dy(x)}{dx}\approx\frac{y(x+ h) - y(x)}{h}

とするのです。これより

\frac{y(x+ h) - y(x)}{h}\approx f(x,y)

とすることができます。

さらに $h$ を両辺にかけ、式を整理することで

y(x+h) \approx y(x) + hf(x,y)

となります。

この式の意味するところは、例えば $x=a$ とすると、$y(a)$ が既知の時、 $y(a+h)$ が $y(a) + hf(a,y(a))$ にほぼ等しいということです。

ただし、あくまで「$\approx$（ニアイコール）」ですので両者には誤差があります。

これで前準備は終了です。

2. オイラー法

いよいよオイラー法の説明に入ります。

改めて、オイラー法とは数値計算を使って常微分方程式の特殊解を求めるアルゴリズムの一つです。ただし、求める解は近似解です。

オイラー法で重要なのは、前準備で作り上げた以下の近似式

y(x+h) \approx y(x) + hf(x,y)

です。これからオイラー法のアルゴリズムを作り上げます。

では、最初にそのアルゴリズムを以下に示します。

これを見て理解ができて、もうプログラムが書けるという人は3章までスキップしてください。

プログラムを書いた後で説明を読んだら分かった、ということもよくありますので、敢えてスキップするものありです。

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

1. 微分方程式の初期値 $y_0=y(0)$ と、積分範囲 $0 \leq x \leq a$ を決める。これによりオイラー法の出発点 $(x,y)=(0, y_0)$ が用意できる。

2. 刻み幅 $h$ によって積分範囲を $n$ 分割する。すなわち $x_{k+1}=x_k + h\ (k=0,1,\dots,n-1),\ x_0 = 0,\ x_n = a$ である。

3. $y_k=y(x_k)$ とすると、$(x_k, y_k)$ を近似式 $y_{k+1} = y_k + hf(x_k, y_k)$ に代入して $y_{k+1}$ を求める。

4. 3 の手順を $y_1$ から $y_n$ が求まるまで繰り返す。

2.1 手順1：初期値と積分範囲の指定

まずは微分方程式の初期値 $y_0$ を決めます。

大抵の場合、初期値とは $x=0$ の時の $y$ の値 $y(0)$ を指します。したがって積分範囲も $0$ を開始として問題ないでしょう。

これにより、オイラー法の出発点 $(x,y)=(0,y_0)$ が用意できます。出発点という言葉の持つニュアンスは後で説明します。

2.2 手順2：刻み幅

刻み幅 $h$ というのは、前進差分の説明で登場した $h$ のことです。極限操作の代用とするための、小さな正の値です。

例えば私は $h=10^{-4}$ で試すことが多いです³。

そして、$h$ を使って積分範囲を $n$ 分割します。横一文字の直線を細かく切っていくようなイメージで、「左端の $x_0=0$」から $h$ 単位で「右端の $x_n = a$」まで刻んでいきます。

それを $x_{k+1}=x_k + h\ (k=0,1,\dots,n-1),\ x_0 = 0,\ x_n = a$ と一般表現しています。

すると $x_0, x_1, x_2, \dots, x_n$ というように一つずつ分割点ができます。

2.3 手順3と手順4

手順3、および4は実際に手を動かした方が理解が早いです。

まず $y_k = y(x_k)$ と表記することにします。

手順1で出発点 $(x_0, y_0)=(0,y(0))$ を決めました。

次に、$(x_0, y_0)$ を近似式 $y_{k+1} = y_k + hf(x_k, y_k)$ に代入すれば $y_1$ が求まります。

$y_1$ が求まったので、同様に $(x_1, y_1)$ を近似式に代入すれば $y_2$ が求まります。

このようにして、$y_n = y(x_n)$ まで順番に一個ずつ求めていきます。

これで $y_0,\ y_1,\ \dots,\ y_n$ まで求めることができました。

後は $(x_0,y_0), (x_1, y_1), \dots, (x_n,y_n)$ を二次元平面にプロットすれば近似解の完成です。

このように、初期値 $y_0$ を決め、さらに積分範囲を分割することで出発点を設定でき、またそこから数値計算によって次の点を一つずつ求めることができます⁴。

この一つずつ求めていく処理方法は逐次処理と呼ばれます。

2.5 オイラー法を一言で説明すると？

常微分方程式の特殊解は一変数関数であり、それは二次元平面上でグラフ（一つの曲線）として表現できます。つまり「点の集合」です。

結局、オイラー法とは「特殊解のグラフを形成する点を出発点から順々に求める手法」なのです。

厳密には、準備した $x$ 座標と初期値 $y_0$ を満たすような $y$ 座標を求める方法でしょうか。

$y$ 座標を求めるための数式が $y_{k+1} = y_k + hf(x_k, y_k)$ なのです。

3. オイラー法の実装

それではオイラー法を実装します。

sample1.m

% 手順1,2
y0 = 1; % 初期値
h = 0.01; % 刻み幅
x = 0:h:10; % 積分範囲

y1 = exp(-2*x); % 真の解
y2 =  EulerMethod(@func, x, y0); % 近似解

plot(x, y1,'LineWidth',1.3) % 真の解の描画
hold on
plot(x, y2,'r--', 'LineWidth',1.3); % 近似解の描画
legend('真の解','近似解')
set(gca,'FontSize',13)

grid on

% 微分方程式
function ydot = func(x, y)
  ydot = -2*y;
end

% オイラー法
function Yout = EulerMethod(func, x, Y0)
  % func:微分方程式の関数ハンドル
  % x：積分範囲（hで分割済みのやつ）
  % Y0：初期値
  nrow = length(Y0);
  ncol = length(x);

  Y = zeros(nrow, ncol);
  Y(:,1) = Y0; % 初期値の代入

  h = x(2) - x(1); % 刻み幅

  % 手順3,4
  for k = 1:ncol-1
      Y(:,k+1) = Y(:,k) + h*feval(func, x(k), Y(:,k)); % 前進差分
  end
  Yout = Y;
end

なお、for文の中身を以下のように書くこともできます（@WolfMoon 様のコードを参考にいたしました。）。

Y(:,k+1) = Y(:,k) + h*func(x(k), Y(:,k));

今回は手計算で解ける簡単な微分方程式にしました。初期値 $y_0$ を $1$ として、真の解は $y=\text{exp}(-2x)$ となります。

実行結果は、真の解と近似解を同時に描画したグラフです。

方程式が簡単だったので、ほぼ誤差なく数値計算が行えたようです。しかし方程式が複雑になれば両者はズレますので、それは頭に入れておいてください。

また、EulerMethod() という名前でオイラー法を実装しました。以下、それについて色々と補足します。

3.1 引数

func は微分方程式の関数ハンドルを受け取る引数です。これはfeval() という関数の中で使用します。これがないと近似式の $f(x_k, y_k)$ が実装できません。

x は積分範囲を受け取る引数です。この時点で $h$ によって分割された状態になっています。つまり $x_0, x_1, \dots,\ x_n$ にした状態で渡したのです。

Y0 は初期値を受け取る引数です。なぜ大文字なのかというと、それは後で説明します。

3.2 feval()

手順3と手順4を行っているfor文の中で、一番気になるのは「feval()」でしょう。見たことない人がいてもおかしくない関数です。

eval は evaluate の略で「評価」という意味です。f は関数を表していて、feval()は「関数の評価」という意味を持つ関数となります。

関数の評価というのは、変数に値を代入してその時の関数の値を出すという意味です。

$y=f(x)$ の $x$ に適当な値を代入すれば $y$ が求まりますよね。そのことです。

結局、feval(func, x(k), Y(:,k))というのは $f(x_k, y_k)$ を表していて、for文の中身は近似式を実装しているのです。

なお、関数ハンドルを使えば feval() を使わなくても、以下のように直接的な書き方も可能です。

補足は以上です。

4. オイラー法を高階微分方程式に適用

最後に、私が初期値を大文字にした理由を説明します。結局、単なる私のこだわりでしかありませんが...。

12回では高階（$n$階）微分方程式を解くために、これを連立一階微分方程式に書き換えることを行いました。

実は、この書き換えた方程式にもオイラー法が適用できます。

仮に、以下のような高階微分方程式

a_n x^{(n)}(t) + a_{n-1} x^{(n-1)}(t) + \dots + a_1 \dot{x}(t) + a_0 x(t) = g(t)

を解きたいとします⁵。解は $x=x(t)$ です。

これを連立一階微分方程式に書き換える方法は以下の通りです。

$n$ 階微分方程式を連立一階微分方程式へ書き換える方法

1. $y_1 = x,\ y_2 = \dot{x},\ \dots,\ y_n=x^{(n-1)}$ とおく
2. $\dot{y_1},\ \dot{y_2},\ \dots,\ \dot{y_n}$ を求める

書き換えたとして、ここで $Y=[y_1, y_2, \dots, y_n]^T$ とひとまとめにします。すると、連立一階微分方程式の一般表現は

\dot{Y} = f(t, Y)

となります。これにより、$x=x(t)$ を求める問題が $Y$ を求める問題に変わりました。

しかし、全体としては $Y$ についての一階微分方程式です。$Y$ を $t$ で微分する際、微分演算は各項について適用されるので、各項にオイラー法を適用すれば問題なく数値計算できます。

結果、オイラー法で高階微分方程式を解くことができるのです。

試しに、12回の時に題材としたファン・デル・ポール方程式にオイラー法を適用してみましょう。

sample2.m

% 手順1,2
Y0 = [1;0]; % 初期値
h = 0.01; % 刻み幅
t = 0:h:10; % 積分範囲

Y =  EulerMethod(@func, t, Y0); % 近似解

plot(t, Y, 'LineWidth',1.3); % 近似解の描画
legend('y_1', 'y_2')
set(gca,'FontSize',13)

grid on

% ファン・デル・ポール方程式（書き換えた後）
function Ydot = func(t, Y)
 Ydot = zeros(2,1);
 mu = 2;
 Ydot(1,1) = Y(2);
 Ydot(2,1) = mu*(1-Y(1)^2)*Y(2) - Y(1);
end

% オイラー法
function Yout = EulerMethod(func, t, Y0)
  % func:微分方程式の関数ハンドル
  % t：積分範囲（hで分割済みのやつ）
  % Y0：初期値
  nrow = length(Y0);
  ncol = length(t);

  Y = zeros(nrow, ncol);
  Y(:,1) = Y0; % 初期値の代入

  h = t(2) - t(1); % 刻み幅

  % 手順3,4
  for k = 1:ncol-1
      Y(:,k+1) = Y(:,k) + h*feval(func, t(k), Y(:,k)); % 前進差分
  end
  Yout = Y;
end

真の解が分からないのでそれとの比較はできませんが、計算結果はでてきました。

ode45()コマンドを使った時と、グラフの形がちゃんと似ているはずです。

まとめ

今回はオイラー法を実装しました。

簡単だったでしょうか、それとも難しかったでしょうか。

一度、理屈が分かってしまえば後はどんな言語でもプログラムを作ることができます。

計算速度にこだわりたいなら、より高い効率を実現するコーディングを追求するのもいいでしょう。

終わりに

オイラー法のように理屈が単純な手法は大体、理解しやすい、実装しやすいことと引き換えに、性能が悪いです。性能が良い手法は大体、理屈が難しくて実装も難しい。

世の中、理屈も実装も簡単で、なのに性能が良い手法ってないんですかね？　一つくらいあってもおかしくないと思うのですが。

もしかして「無ければ作ればいい」というやつでしょうか。

よろしければ次回の記事も読んでくださると大変嬉しいです。

※本記事に対する改善点や修正点、またはこんな事が知りたいといったご意見がありましたらぜひご連絡ください。

参考HP

[1] らい・ぶらり, 数値計算を使って常微分方程式を解く～ルンゲクッタ法の解説～, 2020/01/19公開, 2020/03/06更新, 2022/10/4閲覧

[2] MONOist, オイラー法とルンゲ・クッタ法を使って微分方程式を解こう, 2015/1/28公開, 2022/10/4閲覧

[3] マスオ, 高校数学の美しい物語, オイラー法をわかりやすく解説, 2022/1/05更新, 2022/10/4閲覧

1. 近似という言葉からも分かる通り、近似解と真の解、両者の間には必ず誤差が生じます。その理由は割愛しますが、近年の高性能なコンピュータであってもこの誤差を零とすることはできません。 ↩

2. GPU搭載のノートPCが誕生しているくらいです。 ↩

3. ある文献で得た知識なのですが、微分の真値と前進差分（後退差分と中心差分も）との誤差は $h=10^{-4}$ の時が一番小さいのです。 ↩

4. 厳密には、我々が初期値（出発点）という条件を課して、その条件を満たすような近似解をコンピュータが求めてくれるのです。 ↩

5. 3章の実装例から、階数が高くなり、さらに $y=y(x)$ が $x=x(t)$ に変わりました。 ↩