こんにちは、Rafkaです。
前回 → 数値シミュレーション －数理モデル－
今回は、多くの数理モデルに使用される微分方程式についてのさらなる解説と、本シリーズでの共通の表記法を定義します。
TL;DR
- 微分方程式には大きく分けて、常微分方程式と偏微分方程式があります。
- 微分方程式に使用する文字の定義を行います。
- シリーズで使用するプログラムの作成にはC#を使用します。
- 今後のプログラムで使用するクラスも定義していきます。
微分方程式
微分方程式は、未知の関数とその導関数との関係を表した式で、微分方程式を解くことは、未知関数の具体的な形を決めることを意味します。解く際には積分が必要不可欠となるので、解には積分定数が含まれ、積分定数が含まれた状態で表される解を一般解、初期値等の情報を与えることで積分定数の値まで定めることで得られる解を特殊解と言います。コンピュータで数値的に解く方法は、式変形を行って解く方法ではないので、コンピュータによって求まる解は特殊解となります。
微分方程式には大きく分けて、未知関数が本質的に1変数関数である常微分方程式（Ordinary Differential Equation; O.D.E.）と、未知関数が多変数関数であり、その偏導関数を用いて表される偏微分方程式（Partial Differential Equation; P.D.E.）の2種類があります。それぞれについて解説していきます。
常微分方程式（O.D.E.）
\(n\)階の常微分方程式は１階の連立常微分方程式に変形することができますので、常微分方程式は一般に以下のような形式で表すことができると言えます。
$$\frac{d \boldsymbol{x}(t)}{d t} = \boldsymbol{f} \bigl( \boldsymbol{x}(t), t \bigr)$$
ここで\(t\)は変数、\(\boldsymbol{x}(t)\)は未知関数、\(\boldsymbol{f}(\boldsymbol{x}, t)\)は右辺の式をそれぞれ表していて、未知関数と\(\boldsymbol{f}\)はそれぞれ以下のようなベクトルになります。
$$\boldsymbol{x}(t) = [x_1(t), x_2(t), \cdots, x_n(t)]^T$$
$$\boldsymbol{f}\bigl(\boldsymbol{x}(t), t\bigr) = \Bigl[f_1\bigl(\boldsymbol{x}(t), t\bigr), f_2\bigl(\boldsymbol{x}(t), t\bigr), \cdots, f_n\bigl(\boldsymbol{x}(t), t\bigr)\Bigr]^T$$
常微分方程式の具体例として最も有名なのはニュートンの運動方程式でしょう。
ニュートンの運動方程式は、位置を時間に対する未知関数\(x(t)\)とした時の以下のような２階の常微分方程式であると言えます。
$$m \frac{d^2 x(t)}{d t^2} = F$$
また、位置\(x(t)\)、速度\(v(t)\)、加速度\(a(t)\)には以下のような関係があるので、
$$ \frac{d x(t)}{d t} = v(t), \;\; \frac{d v(t)}{d t} = a(t)$$
運動方程式の例を先程の式の表記に置き換えるとすると、\(\boldsymbol{x}(t)\)と\(\boldsymbol{f}\bigl(\boldsymbol{x}(t), t\bigr)\)はそれぞれ以下のようになる。
$$\boldsymbol{x}(t) = [x(t), v(t)]^T$$
$$\boldsymbol{f}\bigl(\boldsymbol{x}(t), t\bigr) = [v (t), F / m]^T$$
この様にして、２階の常微分方程式であるニュートンの運動方程式を１階の連立微分方程式に変換することができたので、同様にして、高次の導関数を１次の導関数の連結で表現することで、高階の常微分方程式を１階の連立常微分方程式に変換できることが分かりますね。
偏微分方程式（P.D.E.）
常微分方程式が、１変数関数の導関数を用いて構成される方程式であったのに対して、偏微分方程式は多変数関数の偏導関数を用いて構成される方程式になります。
例としては、前回の記事で示した以下のナビエ－ストークス方程式や、シュレディンガーの波動方程式があります。
$$\frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla) \boldsymbol{v} =\; – \frac{1}{\rho} \nabla p + \nu \nabla^2 \boldsymbol{v} + \boldsymbol{g}$$
本シリーズの目標である脳シミュレーションに使用する数理モデルは常微分方程式で表されるので、偏微分方程式は紹介程度に留めさせてもらいます。
プログラム
数値計算を行うプログラムも、このシリーズでは示していこうと思います。
使用する言語はC#7.2になります。
数値計算手法は様々なものがありますが、ベクトルや数理モデルの常微分方程式等は変わらないものなので、その辺りのプログラムを見せていきます。
Vector構造体
今まで示してきた式に出てきたように、\(n\)次元ベクトルを頻繁に使用するので、四則演算をオーバーロードした\(n\)次元ベクトル構造体を以下のように定義して使用していきます。
NumericalModelWithODE抽象クラス
常微分方程式で表される数理モデルを表すクラスで、継承して具体的な数理モデルクラスを構築して使用することを想定しています。
シミュレーションを行う上で必要となる、初期値・開始時間・終了時間の３つはモデルの一部であると考え、プロパティとして持たせました。
肝となる常微分方程式の右辺を計算するメソッドは、モデルごとに異なることや、継承先で必ず実装しなければならないことから抽象メソッドとしました。
これによりクラスも抽象クラスとなり、漠然とした数理モデルというインスタンスが作られることも防げる適した構成だと思います。
Solverクラス
NumericalModelWithODE抽象クラスを継承して構築された数理モデルクラスを解くクラスです。
次回、具体的に説明を行いますが、とりあえずオイラー法という数値計算手法を実装してみました。
他の数値計算手法も同様にこのクラスに実装していきます。
使用例
今まで示してきたプログラムを使って、実際に簡単な例について解いてみましょう。
例題に使うのは以下のような連立微分方程式です。
$$\left\{ \begin{align} \frac{d x_1(t)}{d t} &= x_2(t)\\ \frac{d x_2(t)}{d t} &=\; – x_1(t) \end{align} \right.$$
この方程式の解は、\(x_1(t) = \sin(t)\)、\(x_2(t) = \cos(t)\)となります。それは、
$$\begin{align} \frac{d x_1(t)}{d t} = \frac{d \sin(t)}{d t} &= \cos(t) = x_2(t)\\ \frac{d x_2(t)}{d t} = \frac{d \cos(t)}{d t} &= -\sin(t) = -x_1(t) \end{align}$$
であることから確認できます。
さて、次はこの問題をコンピュータで解いても同じ結果が得られるかどうかを確認してみましょう。
解くのに使用するのは、今まで示してきたプログラムになります。
私は今回、以下のようなプログラムを作成してみました。
簡単にプログラムを解説すると、８行目から１２行目で今回の微分方程式の数理モデルクラスを作成しています。そして１５行目で解く際に使用する刻み幅を\(0.01\)にしていて、１６行目で計算を行いながら結果をファイルに出力しています。
得られた計算結果をプロットしてみると、以下のようなグラフが得られました。
赤線が\(x_1(t)\)、青線が\(x_2(t)\)を表しています。
グラフからは確かに、\(x_1(t) = \sin(t)\)で、\(x_2(t) = \cos(t)\)である様子が確認できます。
これでコンピュータを使用して常微分方程式を数値的に解けるということが分かっていただけたと思います。
まとめ
今回は、数理モデルを記述するための微分方程式について説明し、最後には実際に数値計算によって解く様子をお見せしました。
次回は、今回微分方程式を解くために使用したオイラー法の解説を行います。
次回 → 数値シミュレーション －オイラー法－