計算科学をはじめよう！ニュートンの冷却法則②

前々回の記事では、ニュートンの冷却法則の解析解を導出しました。今回の記事では、いよいよ数値計算によって解を求めます！

オイラー法

オイラー法とは、常微分方程式の数値解を求めるアルゴリズムの1つです。対象とする微分方程式を差分方程式に書き換えて、微分方程式の数値解を求める標準的な方法です。

オイラーのアルゴリズム

オイラーのアルゴリズムを理解するために、ここでは次の1階微分方程式について考えます。

(1) $\begin{align*} \cfrac{dy}{dx} = f(x,y) \end{align*}$

たとえば、 $x=x_0$ において $y$ の値は $y_0$ であるとします。式(1)から $(x_0,y_0)$ における $y$ の傾き（どのように変化するか）がわかりますので、xの変化 $\Delta x$ が十分に小さいとき、 $x_1=x_0+\Delta x$ における $y$ の近似値を求めることができます。ここでは、 $y$ の $x$ についての変化率 $f(x,y)$ が $x_0$ から $x_1$ までの区間で一定であるとします。この仮定により、 $x_1$ における $y$ の近似値は次のように与えられます。

(2) $\begin{align*} y(x_1) = y(x_0) + \Delta y \approx y(x_0) + f(x_0,y_0) \Delta x \end{align*}$

これを一般化するれば、任意の点 $x_{n+1}=x_n+\Delta x$ における $y$ の近似値は漸化式

(3) $\begin{align*} y_{n+1} = y_n + f(x_n,y_n)\Delta x \end{align*}$

で計算することができます。この方法をオイラー法といいます。

ニュートンの冷却法則の数値解

では次にニュートンの冷却法則の数値解をオイラー法を使って求めてみます。前々回の記事で書きましたが、この法則を表す微分方程式と解析解は以下のとおりです。

(4) $\begin{align*} \cfrac{dT}{dt} = -\gamma (T-T_s) \end{align*}$

(5) $\begin{align*} T(t) = T_s + (T_0-T_s)e^{-\gamma t} \end{align*}$

数値解を求めたい(4)式にオイラー法を適用すると、次のようになります。

(6) $\begin{align*} T_{n+1} = T_n - \gamma ( T_n-T_s)\Delta t \end{align*}$

この式を数値的に解けば、(5)式の解が得られるはずです。今回は、以前の記事で扱ったProcessingで数値計算とグラフ作成をしました。以下に解析解と数値計算の結果を比較して示します。各パラメータは、 $T_0=1.0,T_s=0.0,\gamma=1.0,\Delta t=0.05$ としました。

グラフの青い実線が解析解(5)式、赤いプロットが数値計算の結果になります。オイラー法による数値計算の結果が解析解とよく一致しています。

数値計算の安定性

次にパラメータ $\Delta t$ を変化させていったときの結果を示します。

【 $\Delta t=1.0$ の結果】

【 $\Delta t=2.0$ の結果】

【 $\Delta t=2.5$ の結果】

なぜこのような結果になるか、オイラー法のアルゴリズムから考察してみます。

まず、(6)式に今回設定したパラメータを代入します。 $T_s=0.0,\gamma=1.0$ なので、

(7) $\begin{align*} T_{n+1} &= T_n - T_n \cdot \Delta t \\ &= (1-\Delta t) T_n \end{align*}$

となります。次にそれぞれのΔtの場合について考えていきます。

Δt=1.0のとき

このとき、(7)式は次のようになります。

(8) $\begin{align*} T_{n+1} = 0 \end{align*}$

つまり、初期の $T_0$ 以外はすべてゼロになるということです。先ほどグラフで示したように、Δtを1.0まで大きくすると、解析解の曲線を再現できない、ということです。Δtを大きくしすぎると、数値計算では誤差が大きくなってしまいますので、条件の設定には注意が必要です。

Δt=2.0のとき

このとき、(7)式は次のようになります。

(9) $\begin{align*} T_{n+1} = -T_n \end{align*}$

これは $T_{n+1}$ の絶対値は初期状態 $T_0$ から変わらず、符号が反転し続けることを示しています。グラフの赤いプロットで示したように、数値計算の解は振動し続けます。

Δt=2.5のとき

このとき、(7)式は次のようになります。

(10) $\begin{align*} T_{n+1} = -1.5T_n \end{align*}$

これは $T_{n+1}$ の絶対値は初期状態 $T_0$ から増加し続けて、符号も反転し続けることを示しています。このような状態を数値不安定といいます。この状態になると、変数の値が発散して、エラーが起きます。

まとめ

オイラー法によってニュートンの冷却法則の微分方程式を数値計算で解きました。数値計算を使えば、複雑な微分方程式を解くこともできます。ただし、今回の記事で扱ったように、パラメータの設定によっては誤差が大きくなったり、振動や発散が起きてしまいます。数値計算で得られた結果が妥当なものであるか、確認する必要があります。

今後は、他の物理現象に数値計算を適用した例や、オイラー法以外の数値計算のアルゴリズムを紹介していきます！