当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

【Excel】微分方程式の数値解の求め方

【Excel】常微分方程式の数値解

微分方程式の数値解を計算するための初等的なメソッドの一つに、オイラー法があります。オイラー法とは、ある関数 $y=f(x)$ について次のような形をした微分方程式を解くために用いられる手法です。
 \[\frac{dy}{dx}=F(x,y)\]
右辺は $x$ と $y$ の 2 変数関数(もちろん 1 変数でも構わない)なので、かなり幅広い種類の微分方程式を含むことになります。$x$ から少しだけ離れた $x+\Delta x$ における $y$ の値をテイラー展開の最初の 2 項をとって
 \[f(x+\Delta x)=f(x)+f’ (x)\Delta x\]
と近似します。$f’ (x)=F(x,y)$ なので
 \[f(x+\Delta x)=f(x)+F(x,y)\Delta x\]
と書き直せます。初期値 $(x_0,y_0)$ を与えて、そこから少しだけずれた点の $f(x_1)$ を求め、さらにその値をもとに $f(x_2),f(x_3),\cdots$ を計算していけば、$x$ と $y$ のデータが得られます。

例として、微分方程式
 \[\frac{dy}{dx}=x(1-y)\]
の数値解を求めてみましょう。オイラー法では
 \[f(x+\Delta x)=f(x)+x(1-y)\Delta x\]
で近似計算することになります。まず下図のようなシートを作ります。
 
Excel微分方程式ワークシート作成

$\Delta x=0.1$ に設定します。$x$ の値を 0.1 刻みで-4 から 4 まで入力してください。また $x=0$ における $y=f(x)$ の値、すなわち $f(0)$ を 2 としておきます。この $x=0$ の値を基準に順次オイラー法で値を入れるのですが、$x$ が増加する方向では $\Delta x$ が正、すなわち + 0.1 ですが、$x$ が減少する方向では $\Delta x$ は負、すなわち -0.1 となることに注意して、数式を書き分ける必要があります。
 
Excel微分方程式0の前後で数式変更

セル C43 には先ほど記入しておいたセル F3 の初期条件を絶対参照して

=$F$3

と入力します。セル C42 には

=C43+B43*(1-C43)*(-0.1)

という数式を入れてセル C3 までコピーします。セル C44 には

=C43+B43*(1-C43)*0.1

と入力し、右下隅をダブルクリックして行末まで数値を埋めます。入力されたデータ範囲を選択して [挿入] 、[散布図] 、[散布図] の順に選択すると次のようなグラフが表示されるはずです。
 
微分方程式の解曲線グラフ1

$(0,f(0))$ を頂点として正負の方向に減衰する関数です。しかし初期値のとりかたによって、グラフの形(解曲線)は大きく変化します。今度は $f(0)=0$ としてみましょう。
 
微分方程式の解曲線グラフ2

今度は $(0,f(0))$ が谷底となるグラフです。
次は初期値を $f(0)=1$ としてみると…
 
Excelでプロットした微分方程式の解曲線グラフ3

$y=1$ という直線になってしまいます (いわゆる特解です)。
この微分方程式の理論解は初期値を用いて
 \[y=1+(f(0)-1)\exp\left(-\frac{x^2}{2}\right)\]
と表せて、$f(0)=1$ を境に解曲線を 2 つの群に分けることが知られています。変数分離で簡単に解けるので、気になる人は試してみてください。
 
≫ Pythonのサイトでもオイラー法による数値解を扱っています

エクセルや数学に関するコメントをお寄せください