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

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

オイラー法

オイラー法とは、常微分方程式の数値解を求めるアルゴリズムの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の近似値を求めることができます。ここでは、yxについての変化率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から増加し続けて、符号も反転し続けることを示しています。このような状態を数値不安定といいます。この状態になると、変数の値が発散して、エラーが起きます。

まとめ

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

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

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です