ラプラス変換で定数係数の線形非同次微分方程式を解く - ケィオスの時系列解析メモランダム

定数係数の線形非同次微分方程式を解くために必要なラプラス変換の知識を整理します．

　定数係数の微分方程式は，ラプラス変換で解かなくても，他に解き方はいくつかあります．

　ここで強調したいことは，定数係数の微分方程式は必ず解に辿りつけるレシピがあるということです．そして，「なぜ，この方法で解けるの？」という疑問を持って，この記事の内容を超えて理解を深めてほしいです (下の動画と下図を参考にしてください)．

www.youtube.com

定数係数の線形非同次微分方程式

　ここで考えるのは，以下のタイプの微分方程式です．

$\displaystyle{ \frac{d^n y(t)}{d t^n}+a_{n-1} \frac{d^{n-1} y(t)}{d t^{n-1}}+ \cdots + a_2\, \frac{d^2 y(t)}{d t^2}+a_1\, \frac{d y(t)}{d t}+a_0\, y(t)= x(t) }$

ここで， $\displaystyle{ \left\{ a _ 0, a _ 1, \cdots, a _ {n-1}\right\} }$ は定数です．工学で扱うシステムでは， $\displaystyle{ x(t) }$ が入力で， $\displaystyle{ y(t) }$ が出力です．

　工学で扱う微分方程式では，通常，以下のことが仮定されています．　

マイナスの時刻には何も起きない

　工学で考える世界では，すべての物語は， $\displaystyle{ t \ge 0 }$ で展開されます． $\displaystyle{ t \lt 0 }$ には，何も起こらないです．わざわざ，マイナスの時刻なんて考えないです．

　このことを関数で表すために，単位ステップ関数

$\displaystyle{ u(t)= \begin{cases}0 & (t \leq 0) \\ 1 & (t > 0)\end{cases} }$

あるいは，

$\displaystyle{ u(t)= \begin{cases}0 & (t<0) \\ 1 & (t \geq 0)\end{cases} }$

が使われます (私は下の図のような考え方をしたいので，一つ目の定義を好みます)． $\displaystyle{ u(t) }$ をかければ， $\displaystyle{ t \lt 0 }$ での値をすべて0にすることができます．

　また，単位インパルス関数 $\displaystyle{ \delta(t) }$ は，

$\displaystyle{ \int_{0}^{\infty} \delta(t) \, d t=1 }$

です． $\displaystyle{ t }$ がマイナスの領域なんて考える必要はありません．ディラックのデルタ関数と違うと疑問を感じた方は以下の記事を参照してください．

デルタ関数，単位インパルス関数は，どこに立っている？ - ケィオスの時系列解析メモランダム

単位インパルス関数については，

$\displaystyle{ \int_{0}^{\infty} \delta(t)\, f(t) \, d t = f(0) }$

あるいは，

$\displaystyle{ \int_{0}^{\infty} \delta(t-a)\, f(t) \, d t = f(a) \quad (a > 0) }$

のように， $\displaystyle{ \delta(t-a) }$ の括弧の中が0になる値 $\displaystyle{ t=a }$ を抜き出すことができます．

また

$\displaystyle{ \begin{aligned} \int_{0}^{0} \delta(t)\, d t &= 0 \\ \int_{0}^{t} \delta(t) \, d\tau &= 1 \quad (t > 0) \end{aligned} }$

です．

微分方程式を解くときに使うラプラス変換

　微分方程式をラプラス変換の公式（対応表）を使って解くとき，導関数のラプラス変換を使いこなす必要があります．ということで，以下の表の対応関係を印象に残してください．表の左列が元の微分方程式に登場する関数（原関数），右列がラプラス変換後の像関数です．

原関数 $\displaystyle{y(t)}$	像関数 $\displaystyle{Y(s)}$
$\displaystyle{y(t)}$	$\displaystyle{Y(s)}$
$\displaystyle{\frac{d y(t)}{d t}}$	$\displaystyle{s Y(s) - y(0)}$
$\displaystyle{\frac{d^ 2 y(t)}{d t^ 2}}$	$\displaystyle{s^ 2 Y(s) - s y(0) - \frac{d y(0)}{d t}}$
$\displaystyle{\frac{d^ 3 y(t)}{d t^ 3}}$	$\displaystyle{s^ 3 Y(s) - s^ 2 y(0) - s \frac{d y(0)}{d t} - \frac{d^ 2 y(0)}{d t^ 2}}$
$\displaystyle{\frac{d^ 4 y(t)}{d t^ 4}}$	$\displaystyle{s^ 4 Y(s) - s^ 3 y(0) - s^ 2 \frac{d y(0)}{d t} - s \frac{d^ 2 y(0)}{d t^ 2} - \frac{d^ 3 y(0)}{d t^ 3} }$
$\displaystyle{\frac{d^ n y(t)}{d t^ n}}$	$\displaystyle{s^ n Y(s) - s^ {n-1} y(0) - s^ {n-2} \frac{d y(0)}{d t} - \cdots - s \frac{d^ {n-2} y(0)}{d t^ {n-2}} - \frac{d^ {n-1} y(0)}{d t^ {n-1}} }$

一度，ラプラス変換の定義に従って計算してみると，公式のルールが印象に残ると思います (部分積分が使えます)．

ラプラス変換は項ごとに分けて考えて良い

　ラプラス変換は線形な演算なので，項ごとに分けて公式を当てはめてよいことも覚えておいてください．ラプラス変換も逆変換も，足し算と引き算で分けられた項ごとに，考えればよいのです．定数倍は，そのまま定数倍です．

　例えば，

$\displaystyle{ \begin{aligned} \frac{d^2 y(t)}{dt^2} + 3 \frac{d y(t)}{dt} + 2 y(t) &= 0 \\ \end{aligned} }$

の両辺をラプラス変換すると，

$\displaystyle{ \begin{aligned} \left(s^2 Y(s) - s y(0) - \frac{d y(0)}{dt} \right) + 3 \Big(s Y(s) - y(0) \Big) + 2 Y(s) &= 0 \end{aligned} }$

です．これは，

$\displaystyle{ \begin{aligned} \frac{d^2 y(t)}{dt^2} + 3 \frac{d y(t)}{dt} + 2 y(t) &= 0 \\ \end{aligned} }$

に含まれる $\displaystyle{ \frac{d^ 2 y(t)}{dt^ 2} }$ ， $\displaystyle{ \frac{d y(t)}{dt} }$ ， $\displaystyle{ y(t) }$ のそれぞれにラプラス変換の公式を当てはめただけです．係数の３と２は，そのままくっつけておいてください．

初期値が全部0なら $\displaystyle{s^ n}$ をかけるだけ

　　 $\displaystyle{t = 0}$ での初期値が全部ゼロ，

$\displaystyle{ y(0) = 0,\ \frac{d y(0)}{d t} = 0, \ \cdots }$

であれば，ラプラス変換の公式はもっと簡単になります．

原関数 $\displaystyle{y(t)}$	初期値がすべて0のときの像関数 $\displaystyle{Y(s)}$
$\displaystyle{y(t)}$	$\displaystyle{Y(s)}$
$\displaystyle{\frac{d y(t)}{d t}}$	$\displaystyle{s Y(s)}$
$\displaystyle{\frac{d^ 2 y(t)}{d t^ 2}}$	$\displaystyle{s^ 2 Y(s)}$
$\displaystyle{\frac{d^ 3 y(t)}{d t^ 3}}$	$\displaystyle{s^ 3 Y(s)}$
$\displaystyle{\frac{d^ n y(t)}{d t^ n}}$	$\displaystyle{s^ n Y(s)}$

$\displaystyle{ n }$ 階の導関数なら， $\displaystyle{ s }$ の $\displaystyle{ n }$ 乗をかけるだけです．楽ちんです．

その他のラプラス変換の公式は，wikipediaを参考にしてください．

ラプラス変換 - Wikipedia

注意してほしいことは， $\displaystyle{ t \lt 0 }$ のマイナスの時刻なんて考えないので，どの関数にも単位ステップ関数 $\displaystyle{ u(t) }$ がかけてあることです．

インパルス応答：単位インパルスを入力

$\displaystyle{ \frac{d^n y(t)}{d t^n}+a_{n-1} \frac{d^{n-1} y(t)}{d t^{n-1}}+ \cdots + a_2\, \frac{d^2 y(t)}{d t^2}+a_1\, \frac{d y(t)}{d t}+a_0\, y(t)= \delta (t) }$

の解をインパルス応答と呼びます．この式は，単位インパルス $\displaystyle{ \delta (t) }$ を入力したときの出力 $\displaystyle{ y (t) }$ を記述しています．

　何でわざわざ，単位インパルス $\displaystyle{ \delta (t) }$ を入力するの？と疑問に感じるかもしれません．これは，単位インパルスの応答であるインパルス応答がわかれば，任意の入力 $\displaystyle{ x (t) }$ に対する出力がわかるからです．

インパルス応答を求める際は，

　 $\displaystyle{ \delta (t) }$ のラプラス変換は１，

ということを忘れないでください (ラプラス変換の計算をしてみれば当たり前です)．

このことと，以下の対応表を使えば，インパルス応答を求めることができます．

像関数 $\displaystyle{Y(s)}$	原関数 $\displaystyle{y(t)}$
$\displaystyle{\frac{1}{s+\alpha}}$	$\displaystyle{\mathrm{e}^ {-\alpha t} \cdot u(t)}$
$\displaystyle{\frac{\omega}{(s+\alpha)^ 2+\omega^ 2}}$	$\displaystyle{\mathrm{e}^ {-\alpha t} \sin (\omega t) \cdot u(t)}$
$\displaystyle{\frac{s+\alpha}{(s+\alpha)^ 2+\omega^ 2}}$	$\displaystyle{\mathrm{e}^ {-\alpha t} \cos (\omega t) \cdot u(t)}$
$\displaystyle{\frac{1}{(s+\alpha)^ {n}}}$	$\displaystyle{\frac{t^ {n-1}}{(n - 1)!} \mathrm{e}^ {-\alpha t} \cdot u(t)}$
$\displaystyle{\frac{\omega}{s^ 2+\omega^ 2}}$	$\displaystyle{\sin (\omega t) \cdot u(t)}$
$\displaystyle{\frac{s}{s^ 2+\omega^ 2}}$	$\displaystyle{\cos (\omega t) \cdot u(t)}$
$\displaystyle{\frac{1}{s^ {n}}}$	$\displaystyle{\frac{t^ {n-1}}{(n - 1)!} \cdot u(t)}$

一般の入力 $\displaystyle{x(t)}$ での解は畳み込みで求められる

$\displaystyle{ \frac{d^n y(t)}{d t^n}+a_{n-1} \frac{d^{n-1} y(t)}{d t^{n-1}}+ \cdots + a_2\, \frac{d^2 y(t)}{d t^2}+a_1\, \frac{d y(t)}{d t}+a_0\, y(t)= x (t) }$

の解は，

$\displaystyle{ \frac{d^n y(t)}{d t^n}+a_{n-1} \frac{d^{n-1} y(t)}{d t^{n-1}}+ \cdots + a_2\, \frac{d^2 y(t)}{d t^2}+a_1\, \frac{d y(t)}{d t}+a_0\, y(t)= \delta (t) }$

の解であるインパルス応答を $\displaystyle{ h(t) }$ とすれば，次の畳み込み積分で求めることができます．

$\displaystyle{ y(t)=\int_{0}^{t} x(\tau) h(t-\tau) \, d \tau }$

　もちろん，

$\displaystyle{ \frac{d^n y(t)}{d t^n}+a_{n-1} \frac{d^{n-1} y(t)}{d t^{n-1}}+ \cdots + a_2\, \frac{d^2 y(t)}{d t^2}+a_1\, \frac{d y(t)}{d t}+a_0\, y(t)= x (t) }$

の両辺をラプラス変換しても，解を求められます．

ちょっと問題を解いてみる

例題１

$\displaystyle{ \frac{d y(t)}{dt} + 2 \, y(t) = \delta (t) }$

で， $\displaystyle{ y(0) = 0 }$ となる場合を解いてみます (初期値ゼロなら楽ちんです)．

$\displaystyle{ \frac{d y(t)}{dt} + 2 \, y(t) = \delta (t) }$

の両辺をラプラス変換すると，

$\displaystyle{ \begin{aligned} s Y(s) + 2 Y(s) &= 1 \\ (s + 2 ) Y(s) &= 1 \\ Y(s) &= \frac{1}{s + 2} \end{aligned} }$

となります．

　ラプラス変換の対応表で，

像関数 $\displaystyle{Y(s)}$	原関数 $\displaystyle{y(t)}$
$\displaystyle{\frac{1}{s+\alpha}}$	$\displaystyle{\mathrm{e}^ {-\alpha t} \cdot u(t)}$

となっていますので，

$\displaystyle{ \begin{aligned} Y(s) &= \frac{1}{s + 2} \end{aligned} }$

では， $\displaystyle{ \alpha = 2 }$ になってます．したがって，逆変換をすれば答えが，

$\displaystyle{ \begin{aligned} y(t) &= \mathrm{e}^ {-2 t} \cdot u(t) \end{aligned} }$

と求められます． $\displaystyle{ u(t) }$ がくっついている形を見慣れていない人もいるかもしれませんが，これは単位ステップ関数です．別の書き方をすれば，

$\displaystyle{ y(t)= \begin{cases}0 & (t \leq 0) \\ \mathrm{e}^ {-2 t} & (t > 0)\end{cases} }$

となります．この場合，

$\displaystyle{ u(t)= \begin{cases}0 & (t \leq 0) \\ 1 & (t > 0)\end{cases} }$

の定義を採用しました．

例題２

$\displaystyle{ \frac{d y(t)}{dt} + 2 \, y(t) = \mathrm{e}^ {-3 t} }$

で， $\displaystyle{ y(0) = 0 }$ となる場合を解いてみます．

　　まずは，畳み込み積分で計算してみます．例題１でインパルス応答が，

$\displaystyle{ \begin{aligned} h(t) &= \mathrm{e}^ {-2 t} \cdot u(t) \end{aligned} }$

と求まっているので，例題２の答えは，

$\displaystyle{ \begin{aligned} y(t) &= \int_{0}^{t} x(\tau) h(t-\tau) d \tau \\ &= \int_{0}^{t} \mathrm{e}^ {-3 \tau} \mathrm{e}^ {-2 (t - \tau)} \cdot u(t)\, d \tau \\ &= \mathrm{e}^ {-2 t} \int_{0}^{t} \mathrm{e}^ {-3 \tau} \mathrm{e}^ {2 \tau} \, d \tau \\ &= \mathrm{e}^ {-2 t} \int_{0}^{t} \mathrm{e}^ {-\tau}\, d \tau \\ &= \mathrm{e}^ {-2 t} \left[ - \mathrm{e}^ {-\tau} \right]_{0}^{t} \\ &= - \mathrm{e}^ {-3 t} + \mathrm{e}^ {-2 t} \quad (t \geq 0) \end{aligned} }$

となります．

　別の解き方として，

$\displaystyle{ \frac{d y(t)}{dt} + 2 \, y(t) = \mathrm{e}^ {-3 t} }$

の両辺をラプラス変換してみます．

$\displaystyle{ \begin{aligned} s \, Y(s) + 2 Y(s) &= \frac{1}{s+3} \\ (s + 2) Y(s) &= \frac{1}{s+3} \\ Y(s) &= \frac{1}{(s + 2)(s+3)} \\ Y(s) &= \frac{1}{s + 2} - \frac{1}{s+3} \\ y(t) &= \mathrm{e}^ {-2 t} - \mathrm{e}^ {-3 t} \quad (t \geq 0) \end{aligned} }$