11回目のMATLAB記事です。
今回と次回に分けて、「微分方程式」を扱いたいと思います。
微分方程式は、この世に存在する現象（ここでは自然現象や物理現象のこと）を理解、および記述するために、人間が生み出した方法の「一つ」です1。
微分方程式は物理学をはじめ、多くの分野で活躍しています。
例えば、力学においては乗り物の動き、投げたボールの動き等、これらの動きを全て微分方程式で記述します。
他にも、微分方程式の応用例は枚挙にいとまがありません。
本記事は前編として、まず微分方程式の基本事項を簡単に説明します。
その後、微分方程式をMATLABで定義する方法と、MATLABでの解き方を説明します。
微分方程式
微分方程式とは、ある関数 $y$ とその導関数 $y^{'}$ を含む方程式のことです。
例えば、ある方程式が $x$ の関数 $y=y(x)$ と、その導関数 $\frac{dy}{dx}$ を含んでいる時、その方程式は微分方程式です。
微分方程式の例
\begin{equation}
\frac{dy}{dx} = -3y
\end{equation}
方程式の主旨は、その等式（方程式）をみたす解を求めることです。
微分方程式における解とは「関数 $y$」のことです。
微分方程式をみたす関数を求めることを、微分方程式を解くといいます。
ちなみに、上の例における微分方程式の解は以下の通りです。
\begin{equation}
y=y(x)=C\exp{(-3x)}\ (C \in \mathbb{R})
\end{equation}
ただし、$C$ は任意の定数です。
一般解と特殊解
微分方程式の解には、一般解と特殊解の二種類があります。
解の種類
- 一般解：微分方程式の全ての解のこと。任意定数を含む。
- 特殊解：一般解のうち、初期条件2を満たす解。
例えば、ある二次方程式を満たす解は、複素数の範囲に重複を含めて二個存在します。
微分方程式の場合、解は無数に存在します。
なぜでしょうか。
微分方程式には導関数が含まれるので、解を求めるには積分が必要で、その際、任意定数 $C$ が出てきます3。
先ほどの微分方程式の例
\begin{equation}
\frac{dy}{dx} = -3y
\end{equation}
においても、解となる関数へ任意定数が含まれてましたね。
そこで、微分方程式の全ての解を総称して「一般解」と呼ぶのです。
さて、解が無数に存在しては非常に扱いにくいので、任意定数の値を具体的に決める条件を与えます。
上記の微分方程式の一般解である
\begin{equation}
y=y(x)=C\exp{(-3x)}
\end{equation}
へ含まれる $C$ を決めるには、どのような条件が必要なのでしょうか。
それは「$y(0)$ の値」です。
これを初期条件、あるいは初期値といいます。
このように、ある初期値を満たす解（関数）$y$ を求める問題を初期値問題といいます。
微分方程式を解く、というのは基本的に初期値問題のことを指します。
微分方程式の階数
最大で $n$ 次の導関数が登場するような微分方程式を $n$ 階微分方程式といいます。
例えば $y^{(2)} = -2y$ は二階微分方程式です。
線形微分方程式
$y^{2}, yy^{'}$ といった項を含まない以下のような微分方程式を線形微分方程式といいます。
\begin{equation}
P_n(x)y^{(n)}+P_{n-1}(x)y^{(n-1)}+\dots+P_1(x)y^{'}+P_0(x)y=Q(x)
\end{equation}
$P_i (x)$ がどのような関数であっても、線形微分方程式です。
注目しているのは、$y$ やその導関数の次数なのです。
同次、非同次
線形微分方程式は $Q(x)$ の値によって以下のように分けることできます。
-
同次微分方程式：$Q(x)=0$ のもの。
-
非同次微分方程式：$Q(x) \neq 0$ のもの。
このうち、線形同次微分方程式には重ね合わせの原理という性質が成り立ちます。
重ね合わせの原理
-
ある線形同次微分方程式の解が $y_1$ であれば、そのスカラ倍 $cy_1$ も解である。
-
ある線形同次微分方程式の解が $y_1, y_2$ であれば、それらの和 $y_1 + y_2$ も解である。
これが成り立つのは非常に嬉しいことです。
元々、微分方程式の目的は現象の記述、および理解することです。
重ね合わせの原理が成り立つと非常に解析がしやすくなります。
さらに嬉しいことに、物理現象の多くが線形同次微分方程式で表すことができます。
自由落下運動、バネの運動、電気回路の回路方程式がその代表例です。
常微分方程式と偏微分方程式
これまで、一変数関数とその導関数を含む微分方程式について考えました。
しかし、実際には多変数関数とその偏導関数を含む方程式を考える場面もあります。
ここで、前者を常微分方程式、後者を偏微分方程式といいます。
-
常微分方程式：一変数関数 $y(x)$ と、その導関数 $y^{'}, y^{(2)},y^{(3)}\dots$ を含む微分方程式のこと。
-
偏微分方程式：多変数関数 $y(x_1, x_2, x_3,\dots,x_n)$ と、その偏導関数 $y_{x_1},y_{x_2},\dots,y_{x_1 x_1}, y_{x_2 x_2},y_{x_1 x_2},\dots$ を含む微分方程式のこと。
普段、微分方程式というと大概は常微分方程式のことを言います。
私の記事でも微分方程式とは常微分方程式を指すことにします。
ところで、偏微分方程式で表す現象の一つが熱です。
熱は座標と時間を変数に持つ多変数関数で表します。
興味のある人は「熱伝導方程式」を調べてみてください。
数値的に解く
さて、微分方程式には大きな問題があります。
それは基本的に手計算で解けないということです。
この世には膨大な個数の微分方程式が存在します。
その中で、手計算で解けるものは氷山の一角といってよくて、まさに奇跡なのです。
微分方程式の元々の目的は、現象を表現し、さらに理解することです。
しかし、せっかく微分方程式で現象を表せたのに、解が分からなければ先に進めません（解のことを真の解といいます）。
これをうけて、人間が編み出したのがコンピュータで解を求めることです。
ただし、コンピュータで求めた解は本物ではなく近似的に求めたもので、これを近似解といいます。
プログラミングで微分方程式を解く場合、導かれるのは近似解です。
繰り返しになりますが、手計算で解ける微分方程式は氷山の一角です。
したがって、自然と、あるいは否が応でも近似解に頼る場面が多くなります。
一階微分方程式の実装
それでは、MATLABで微分方程式を実装してみましょう。
MATLABでは、関数の形式で実装します。
題材は以下の一階微分方程式です。
\begin{equation}
\frac{dy}{dx} = -3y
\end{equation}
これをMATLABでは次のように実装します。
function dy_dx = func(x, y)
dy_dx = -3*y;
end
これは x, y という引数を受け取って、微分 dydx の値を返す関数となります。
x, y はそれぞれ $y=y(x)$ の $x$ と $y$ のことです。
今回は方程式中に $x$ の項が現れていないので結局、引数の x を使いませんでした。
しかし、使わないと実装できない微分方程式も存在します。
一階微分方程式を解く
次に微分方程式を解いてみます4。
MATLABで微分方程式を解くためのコマンドはいくつか存在します。
ここでは、私が普段使っている ode45() というコマンドを使います。
x = 0:0.1:5; % 積分範囲
y0 = 3; % 初期値
[x,y] = ode45(@func, x, y0);
plot(x,y,'b-')
hold on
grid on
function dy_dx = func(x, y)
dy_dx = -3*y;
end
まず一行目で積分範囲を指定します。
数学では不定積分という範囲を指定しない積分もありますが、プログラミングでは基本的に範囲を指定しなければなりません。
二行目では初期値を $y(0)=1$ と指定します。
ちなみにこの時の任意定数の値は $C=1$ になります。
三行目では ode45() コマンドを使って微分方程式を解きます。
引数は左から「微分方程式を表す関数」「積分範囲」「初期値」です。
関数のところにある@マークを忘れないでくださいね。
そして戻り値を受け取る [x,y] には、左から [指定した積分範囲, 解となる関数の積分範囲における値] が格納されます。
もちろん、初期条件をちゃんと満たしています。
最後に plot() によって解を描画してみると次のようなグラフになります。
これは特殊解 $y=\text{exp}(-3x)$ のグラフです。
長くなりますので今回はここで区切りたいと思います。
まとめ
今回の前編では、微分方程式の基本事項を簡単に説明した後、MATLABで一階微分方程式とそれを解くプログラムを実装しました。
次回の後編では、二階以上の微分方程式を解いてみたいと思います。
終わりに
私が初めてMATLABにふれたのは、大学のとある講義でした。
当時はMATLAB完全初心者な上、今ほどプログラミングの知識も経験もありませんでした。
その講義における最初の課題が、実は微分方程式を解くことでした。
同級生や先輩に手伝ってもらい、四苦八苦してなんとか完成させたのを今でも覚えています。
よろしければ次回の記事も読んでくださると大変嬉しいです。
※本記事に対する改善点や修正点、またはこんな事が知りたいといったご意見がありましたらぜひご連絡ください。
参考HP