Pythonで微分方程式とラプラス変換

よびのり先生のYouTube動画で、物理工学の教授二人と話しているのがあり、その中でこう語られていました。

• 理学部では主にフーリエ変換をやるが、工学部では主にラプラス変換をやる
• 耳も目もフーリエ変換器。光コンピュータでも高速フーリエ変換をする
• 制御論はラプラス変換の世界
• 物事を思う通りに動かすには制御が必要で、それはラプラス変換の世界である。
• 実時間応答・フィードバックをしようとしたら、ラプラス変換なしでは語れない。
• ラプラス変換は、微分方程式を解くための道具でとどまらない。
• 物事を制御しようとしたらラプラス変換をして式で動的な応答を追ってみて、「ここをこうすれば安定する」などを、実験を色々と変えたりしなくてもわかる極めて強力なツールである。
• すくなくとも工学部の学生は機械でも電気でも、ラプラス変換で飯を食っているようなものなので、ラプラス変換がわからないとそもそも話にならない。
• 建築でも大きな建造物は制御論で成り立っている。
• ものづくり全体がラプラス変換・制御論で成り立っている。制御論を学ばないと、ちゃんとしたものは作れない。
• ラプラス変換は若いうちに学ぶと深い部分に入っていける

・・・・なんかすごそう！

ということで、ラプラス変換について理解を深めていきたいと思います。

ラプラス変換は微分方程式を簡単にするとはどういうことなのか。まずは微分方程式のややこしさを理解していきます。

微分方程式

微分方程式とは、$\frac{dy}{dx} = y$を満たす関数を見つけるための方程式。この解は$ y = ce^x 　(c: 任意数) $です。

では問題を解いていきます。

例1

$$ y' = 2x+3, y(0) = 0(初期条件) $$

解:

$ y = x^2 + 3x + c　(一般解) $
$ y(0) = 0 $より
$ y = x^2 + 3x （特殊解）$

例2 :

$$ y^{\prime\prime} = e^x + 1, y(0) = 0, y'(0) = 1 $$

解：

$ y' = e^x + x + c_1 $
$ y = e^x + \frac{1}{2} x^2 + c_1 x + c_2 $
$ y(0) = 0 $より、$ 1 + c_2 = 0 $　$\therefore c_2 = -1$
$ y'(0) = 1 $より、$ 1 + c_1 = 1 $　$\therefore c_1 = 0$
以上より
$ y = e^x + \frac{1}{2} x^2 - 1 $

ここまでの微分方程式はかなり単純な微分方程式でした。

ただし、実際の現象をモデル化して微分方程式（＝数理モデル）を得て、それを解いて得た解を見てまた微分方程式を修正して、というのを繰り返して得られた微分方程式は、かなり複雑で解くのがかなり難しいものばかり。

まずは簡単な人工モデルを見ていきます。

例：人工モデル 「個体数の増加率はその時の個体数に比例する」（マルサスの法則）

つまり、

$$ y' = ay, y(0) = y_0 $$

これを解くと、

$ y = ce^{at} $
$ y(0) = y_0 $より、$ c = y_0 $
$ \therefore y = y_0 e^{at} $

この関数は指数関数的増加をするが、実際には資源に限りがあるので、微分方程式を以下のように修正する必要がある。

$$ y' = ay - by^2 $$

これを解くと、最初は指数関数的増加をしますが、あるところで一定になっていきます。

人口モデルはシンプルでしたが、たとえば人工衛星についての微分方程式はかなり複雑です。

こちらによれば、以下の常微分方程式（未知関数が本質的にただ一つの変数を持つ微分方程式）を解くと、木星の周りを回る人工衛星の位置と速度がわかります。

てすと

$$ \frac{d^2x}{dt²} = - \frac{GM(x + tV_j)}{((x + tV_j)² + y²)^{\frac{3}{2}}} $$

$$ \frac{d^2y}{dt²} = - \frac{GMy}{((x + tV_j)² + y²)^{\frac{3}{2}}} $$

$$ \frac{dx}{dt} = v_x $$

$$ \frac{dy}{dt} = v_y $$

$G$…万有引力定数、$ M $…木星の質量、$V_j$…木星の平均公転速度

ラプラス変換

ラプラス変換とは、座標変換の一種で、「微分方程式が簡単になる」という大きなメリットがあります。

微分方程式を解くのは大抵の場合難しいですが、微分方程式をラプラス変換して$s$世界に関する解を求めた後、それをラプラス逆変換して$t$世界に戻せば、微分方程式を解くことができるというわけです。

では、ラプラス変換とは何か。ラプラス変換は「関数は指数関数$e^{st}$のまとまりで表せる」という大胆な思想です。ちなみにフーリエ変換は「関数は三角関数のまとまりで表せる」という思想により、$t$（時間）空間から$\omega$（周波数）空間に移すもので、テイラー展開は「関数はその関数の色々な次数関数のまとまりで表せる」という思想でした。

微分方程式のひとつである運動方程式を考えてみます。

$$
\begin{equation}
\left\{ \,
\begin{aligned}
& v = \frac{dx}{dt} \\
& m\frac{dv}{dt} = F　(運動方程式) \\
& \frac{dv}{dt} = -\frac{k}{m}x　(バネの周期的な動きを表す微分方程式) \\
\end{aligned}
\right.
\end{equation}
$$

これより、

$$
\begin{equation}
\frac{d^2x}{dt^2} = -\frac{k}{m}x \\
\end{equation}
$$

となりますが、２回微分すると元に戻って$-\frac{k}{m}$という係数が前につく関数が$x$ということになります。２回微分すると元に戻る関数を考えると、たとえば

$$
\begin{equation}
\begin{aligned}
& sin(t) \\
& cos(t) \\
& Ae^{st}　(A=定数, s=\pm{i \sqrt{\frac{k}{m}}})
\end{aligned}
\end{equation}
$$

などがあります。なので、$x$は指数関数とか三角関数だろうと予想できますが、今のは計算ではなく、代入法です。頭が良い人ならわかるかもしれませんが、計算でわかるようなものではなさそうです。

そこで関数を指数関数のまとまりで表すことで、その関数の微分は指数関数$e^{st}$の$s$がおりてきた一次関数に、2階微分は指数関数$e^{st}$の$s$が２回おりてきて二次関数になり、一次関数や二次関数なら誰でも解くことができるようになります。つまり、「ラプラス変換により微分方程式を$s$の$n$次関数で表すことができる」ということです。

またオイラーの公式$ e^{i\theta} = cos\theta + i sin\theta $より、指数関数は複素数を使うと、三角関数を含むようになります。

関数を指数関数のまとまりで表すことで表すことができるということは、

$$ f(t) = \int F(s) e^{st} ds 　(原関数)$$

と表せますが、この$F(s)$を計算する方法のことをラプラス変換といいます。

$F(s)$だけを取り出すには$e^{st}$を消す必要があるので、ラプラス変換は以下のような数式で表せます。

$$ F(s) = \mathcal{L}[f(t)] = \int^{\infty}_0 f(t)e^{-st} dt 　(像関数)$$

$f(t)$が$F(s)$になってます。つまり、$t$の関数が$s$の関数に変換されています。

たとえば水力発電を例に取ってみると、川の流れの速さを$x$、水車の回る速度が$y$、モーターで発生する電圧が$z$、得られる電力が$w$とする時、普通$x$と$y$の間、$y$と$z$の間、$z$と$w$の間は微分方程式によって繋がれていますが、これが$x$から$w$を得る際に大変であり、ラプラス変換を使えば簡単に得られるということです。

$f(t) = 1$をラプラス変換すると、$F(s) = \frac{1}{s}$
$f(t) = t$をラプラス変換すると、$F(s) = \frac{1}{s^2}$
$f(t) = \frac{t^2}{2}$をラプラス変換すると、$F(s) = \frac{1}{s^3}$

となります。つまり、$f(t)$を$t$で微分すると、$ F(s)$を$s$倍することになっています。なので、

$f(t) = sin(t)$をラプラス変換すると、$F(s) = \frac{1}{s^2+1}$であれば
$f(t) = cos$のラプラス変換は、$F(s) = \frac{s}{s^2+1}$となります。

時間領域で解析しようと思っていたけれど微分方程式が難しいので、微分法則が簡単になる$F(s)$だけ見れば良い、ということができます。

$f(t) = x(t)$をラプラス変換すると、$F(s) = X(s)$とすると、
$f(t) = \dot{x}(t)$をラプラス変換すると、$F(s) = sX(s)$になるため、$x(t)$の部分積分を頑張らなくても、$x(t)$のラプラス変換を求めていればそれを$s$倍すれば良いだけ。つまり、

$$ \mathcal{L}[\dot{x}] = s\mathcal{L}[x] $$

が成り立ちます。

例題1

$$ \ddot{x} - 3\dot{x} + 2x = 3t^2 $$

この式をラプラス変換すると、

$ s^2X - 3sX + 2X = \frac{6}{s^3} $
$ (s^2 - 3s + 2)X = \frac{6}{s^3} $
$ X = \frac{1}{(s^2 - 3s + 2)} * \frac{6}{s^3} $

と求められます。そしてもし

$$ \ddot{x} - 3\dot{x} + 2x = f(t) $$

とするなら、そのラプラス変換による解は

$ X(s) = \frac{1}{(s^2 - 3s + 2)} * F(s) $

と表せます。$\frac{1}{(s^2 - 3s + 2)}$の部分は伝達関数といい、$G(s)$として

$$ X(s) = G(s)* F(s) $$

と表すこともあります。このように入出力の関係が見やすくなることがわかります。

そしてこれを時間領域$t$に逆変換すれば微分方程式を簡単に解くことができます。しかしラプラス変換のすごさはそういうことではなく、$ G(s) $がわかれば、システムの安定性がだいたいわかるということにあります。

つまり伝達関数$ G(s) $は制御システムに対する入力と出力の関係を表すものであり、入出力システムの挙動や安定性を評価することができる指標として考えることもできるということです。

このような入出力の関係を表す関数には入力信号の種類に関係なく、入力の設定部やセンサによる検出部、必要な制御入力を作る調節部などの制御要素に固有の関数となる事が要求されます。

ではシステムの安定性とは何かといえば、それは「入力なしの時に出力が収束するかどうか」と言えます。

Pythonでマルサスの法則を見てみる

マルサスの法則をPythonの記号計算ライブラリsympyを使って見てみます。

まずはライブラリのインポートと設定、変数や関数の定義をしていきます。

import sympy as sym
from sympy.plotting import plot
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
sym.init_printing(use_unicode=True)
%matplotlib inline

x,y = sym.symbols("x y")
f = sym.Function('f')

マルサスの法則は$ y' = ay, y(0) = y_0 $でした。これを実装していきます。

a = 0.0078 # 日本の人工増加率

eq = sym.Eq(f(x).diff(x,1) - a*f(x))
eq

$ - 0.0078 f{\left(x \right)} + \frac{d}{d x} f{\left(x \right)} = 0 $

# 一般解
ans = sym.dsolve(eq)
ans

$f{\left(x \right)} = C_{1} e^{0.0078 x}$

# 特殊解
f_0 = 123731000 # 日本の人口
ans = sym.dsolve(eq, ics={f(0):f_0})
ans

可視化してみる

plot(ans.rhs, (x, 9000, 12000))