Laplace変換による微分方程式の解法

同義語:Laplace変換法solving ODEs by Laplace transform

概要

Laplace変換による微分方程式の解法(solving ODEs by Laplace transform)とは、定数係数の線形初期値問題を、初期値の項を含む変換領域の代数方程式へ移して解く方法である。入力が指数オーダーなら解も変換可能であり、逆変換の一意性と初期値問題の一意性によって得た解を正当化できる。連立一次系にも同じ方法が使える。

$$\newcommand{C}[0]{\mathbb{C}} \newcommand{div}[0]{\mathbin{÷}} \newcommand{N}[0]{\mathbb{N}} \newcommand{Q}[0]{\mathbb{Q}} \newcommand{R}[0]{\mathbb{R}} \newcommand{Z}[0]{\mathbb{Z}} $$

前提知識: Laplace変換, 逆Laplace変換, 線形微分方程式

Laplace変換による解法では、初期値問題を時間側で解く代わりに、未知関数の変換についての代数方程式を解く。微分の変換には初期値の項が自動的に入り、定数係数はそのまま係数として残る。得られた式を部分分数に分ければ、指数関数や振動の和に戻せる。途中で入力が変わる問題でも、入力の式だけを変換して同じ分母を使える点が便利である。
一方、変換を取る前には解が存在し、変換できる増大率をもつことを確かめる必要がある。求めたい未知関数に都合のよい収束を仮定して、その仮定を使って存在を結論するのでは循環する。本記事では、定数係数の初期値問題が指数オーダーの解をもつことを先に示し、その後に変換の計算を正当化する。一般的な存在・一意性と基本解系は線形微分方程式の定理「解の存在と解空間」「定数係数の斉次方程式の基本解系」で扱われている。

問題の設定

定数係数の線形初期値問題

定数 $a_0,\ldots,a_n$、$a_n\ne0$ と、$[0,\infty)$ 上の連続な入力 $g$ に対し、
$$ \sum_{k=0}^na_ky^{(k)}(t)=g(t),\qquad y^{(j)}(0)=b_j\quad(0\le j< n) $$
を満たす $C^n$ 級の関数 $y$ を求める問題を考える。特性多項式は $P(s)=\sum_{k=0}^na_ks^k$ である。

入力 $g$ は指数オーダーであると仮定する。定数係数であることは、変換後に $P(s)Y(s)$ が現れるために本質的である。係数が $t$ に依存すると、未知関数に掛けた積の変換をさらに調べる必要があり、一般には代数方程式一つへは直らない。非線形の $y^2$ を含む場合も、$\mathcal L[y^2]=Y^2$ とはならないので同じ手順をそのまま使えない。
最高階の係数を1にするには方程式全体を $a_n$ で割る。初期値はその操作で変わらないが、入力は $g/a_n$ になる。力学の例で $my''+cy'+ky=f$ を扱うなら、$m$ がこの最高階係数に当たる。変換後の分母から $m$ を落とすと応答の振幅が変わってしまう。以下の数値例は時間と未知量を無次元化した方程式として扱う。
区分的に連続な入力に対しては、解の最高階微分も区分的になる。$y,\ldots,y^{(n-1)}$ は連続で、方程式は入力の連続点で成立する。この拡張では単位ステップ関数の定理「ステップ入力に対する接続条件」を使える。入力にデルタがある場合にはさらに跳びの条件が必要であり、通常の $C^n$ 級の初期値問題とは区別する。

変換してよいことの確認

初期値問題の解は指数オーダー

定義の入力 $g$ が連続で指数オーダーなら、初期値問題の解 $y$ とその導関数 $y',\ldots,y^{(n)}$ は全て指数オーダーである。従って十分右の半平面で微分の変換公式を全ての項に適用できる。

$x=(y,y',\ldots,y^{(n-1)})^T$ とすれば、定数の同伴行列 $A$ とベクトル $b=(0,\ldots,0,1/a_n)^T$ によって $x'=Ax+bg$ になる。線形微分方程式の定理「定数変化法の公式(Duhamel の公式)」から、解は
$$ x(t)=e^{At}x(0)+\int_0^t e^{A(t-u)}b g(u)\,du. $$
最大行和ノルムを使い $K=\|A\|$ とする。行列指数の級数と劣乗法性から $\|e^{At}\|\le\sum_{m=0}^\infty K^mt^m/m!=e^{Kt}$。$\lvert g(u)\rvert\le Me^{au}$ を全域で成り立つように取り、$c=\max(K,a,0)$ とすれば
$$ \|x(t)\|\le e^{Kt}\|x(0)\|+M\|b\|\int_0^te^{K(t-u)+au}\,du \le (\|x(0)\|+M\|b\|t)e^{ct}. $$
任意の $\varepsilon>0$ について $t\le e^{\varepsilon t}/\varepsilon$ なので、右辺は $C_\varepsilon e^{(c+\varepsilon)t}$ 以下となる。各成分、すなわち $y,\ldots,y^{(n-1)}$ は指数オーダーである。方程式を $y^{(n)}=(g-\sum_{k< n}a_ky^{(k)})/a_n$ と解けば、有限個の指数オーダー関数の線形結合なので $y^{(n)}$ も指数オーダーになる。連続性とともに微分公式の条件が満たされた。$\square$

この評価は鋭い増大率を求めるためのものではない。安定な方程式でも $\|A\|$ は正であり、評価の上界は増大することがある。それでも変換の存在を保証するには十分である。実際の長時間の振る舞いは分母の根を調べて判断する。ノルムによる粗い上界と、根から得る厳密な減衰率を同じものだと解釈しない。
行列の成分に単位がある場合、異なる成分をそのまま最大ノルムで比較するのは物理量としての意味をもたない。存在の議論では各変数を固定した単位で無次元化して数値として扱える。厳密な物理量の評価が目的なら、各成分のスケールを決めてからノルムを定義する。解法の代数式自体は、同じ単位で全ての項をそろえた方程式に適用すればよい。

変換後の方程式と手順の正しさ

初期値を含む代数方程式

$G=\mathcal L[g]$、$Y=\mathcal L[y]$ とし、
$$ R(s)=\sum_{k=1}^na_k\sum_{j=0}^{k-1}s^{k-1-j}b_j $$
とおく。十分右の半平面で
$$ P(s)Y(s)=G(s)+R(s),\qquad Y(s)=\frac{G(s)+R(s)}{P(s)} $$
が成り立つ。右辺の逆変換として連続な指数オーダー関数 $z$ を得られたなら、$z$ は初期値問題の解 $y$ に等しい。

prop-lto-orderによって全ての必要な変換が存在する。Laplace変換の公式「微分の公式」により、各 $k\ge1$ について
$$ \mathcal L[y^{(k)}]=s^kY-\sum_{j=0}^{k-1}s^{k-1-j}b_j. $$
これを $a_k$ 倍して足すと、$Y$ の係数は $\sum_{k=0}^na_ks^k=P$、初期値の項は $-R$ になる。方程式の両辺を変換すれば $PY-R=G$。多項式 $P$ の根は有限個なので、共通の収束半平面を必要ならさらに右へ狭め、その中で $P\ne0$ とできる。そこで割れば $Y=(G+R)/P$ を得る。
候補 $z$ の変換がこの右辺なら、その半平面の十分大きな実数上で $\mathcal L[z]=\mathcal L[y]$。両者は連続で指数オーダーなので、逆Laplace変換の定理「Lerchの一意性定理」により $z=y$ が全ての時刻で成立する。解 $y$ の存在は変換を取る前に確保しているので、この議論は逆変換の存在を仮定して解の存在を証明する循環にはなっていない。$\square$

この正当化では、候補を毎回 $n$ 回微分して方程式へ代入する必要はない。ただし具体的な計算の符号や係数を確かめるには、代入と初期値の確認は有効である。逆変換の表を誤って使った候補は、変換が本当に $Y$ になることを満たさない。方程式の残差をゼロにする確認と、初期値を満たす確認は別に行う。斉次解を足した誤った候補は方程式を満たしても初期値を満たさないことがある。
初期値の多項式 $R$ の最高次数は $n-1$ である。2階なら $R=a_2(sb_0+b_1)+a_1b_0$、3階なら $R=a_3(s^2b_0+sb_1+b_2)+a_2(sb_0+b_1)+a_1b_0$ である。最高階の微分だけに初期値を入れるのでなく、各微分項が生む分を全て足す。右辺へ移すと符号は正になる。

方程式の項変換適用条件・根拠
$a_0y$$a_0Y$変換の線形性
$a_1y'$$a_1(sY-b_0)$微分の公式
$a_2y''$$a_2(s^2Y-sb_0-b_1)$微分の公式
$a_ny^{(n)}$$a_n(s^nY-\sum_{j< n}s^{n-1-j}b_j)$微分の公式
$u(t-a)g(t-a)$$e^{-as}G(s)$$a\ge0$、時間推移の公式

解法の計算の基本は LebDQ26 §6.2 にある。定数変化法と比較すると、時間側で入力を積分する操作が、変換側では $G/P$ という積・商にまとまる。初期値がゼロなら $R=0$ なので、入力の変換に $1/P$ を掛けるだけでよい。この $1/P$ の逆変換がインパルス応答であり、入力との畳み込みで解を表せる。

計算の手順

  1. 方程式を定数係数の標準形へ直し、入力と初期値を全て書く。変換可能な入力であるか確かめる。
  2. 微分の変換を各項に適用し、$P(s)$ と $R(s)$ を作る。$R$ は初期値の項だけ、$G$ は入力だけから来る。
  3. $Y=(G+R)/P$ を通分して約分し、単根・重根・複素共役根を判定する。遅延があれば因子 $e^{-as}$ を外に保つ。
  4. 逆変換し、共役な項を実数の正弦・余弦へまとめる。代表時刻の値を計算する。
  5. 初期値と方程式への代入を確かめる。ステップ入力なら切替点の左右の接続も確認する。
    分母の根が重なると $t^ke^{at}$ が現れる。共振も同じ仕組みで説明できる。入力の変換にある極が $P$ の根と重なると、$Y$ により高位の極が生まれ、時間側の多項式次数が増える。これは線形微分方程式の命題「指数関数を右辺にもつ方程式の解」と同じ結果を変換側から読むものであり、共振の一般定理を別に作る必要はない。
    初期値を含む分子と入力を含む分子を最初から分けておくと、異なる初期値や入力で同じ系を調べる際に再計算を減らせる。例えば $y''+3y'+2y=g$ の分母はいつも $(s+1)(s+2)$ であり、初期値だけ変えれば $R$ だけが変わる。入力が複数のパルスの和なら、線形性で各入力の応答を足す。

例題

非ゼロの初期速度

$y''+3y'+2y=e^{-3t}$、$y(0)=0$、$y'(0)=1$ を解く。答えは $y(t)=\tfrac32e^{-t}-2e^{-2t}+\tfrac12e^{-3t}$ であり、$y(1)\approx0.306042$。thm-lto-algebraと部分分数を使う。

途中の計算を開く

$G=1/(s+3)$、$R=1$。$Y=(1+1/(s+3))/((s+1)(s+2))=(s+4)/((s+1)(s+2)(s+3))$。単根の係数は $3/2,-2,1/2$。$y(0)=0$、$y'(0)=-3/2+4-3/2=1$。

共振による時間因子

$y''+4y=\cos2t$、$y(0)=y'(0)=0$ の答えは $y(t)=\tfrac14t\sin2t$。$y(1)\approx0.227324$。重なった極を含む逆変換を使う。

途中の計算を開く

$Y=s/(s^2+4)^2$。$\mathcal L[t\sin2t]=4s/(s^2+4)^2$ なので係数は $1/4$。$y''+4y=\cos2t$ を微分して確かめられる。

3階の初期値多項式

$y'''+6y''+11y'+6y=0$、$y(0)=1$、$y'(0)=y''(0)=0$ の答えは $y(t)=3e^{-t}-3e^{-2t}+e^{-3t}$。$y(1)\approx0.747420$。

途中の計算を開く

$P=(s+1)(s+2)(s+3)$、$R=s^2+6s+11$。単根の展開係数は3、−3、1。初期値は $3-3+1=1$、$-3+6-3=0$、$3-12+9=0$ になる。

各例の入力と初期値は違うが、手順は同じである。1題目では入力の極 $-3$ が新しい指数を加える。2題目では入力の極と特性根が一致し、時間因子 $t$ が加わる。3題目は入力がゼロでも分子に初期値が残ることを示す。初期値がゼロだから斉次方程式の解もゼロになる場合と、非ゼロの初期値による自由応答を混同しない。

連立方程式

連立1階の変換

定数行列 $A$、連続で指数オーダーのベクトル入力 $b(t)$ に対し、$x'=Ax+b(t)$、$x(0)=x_0$ の変換は
$$ (sI-A)X(s)=x_0+B(s),\qquad X(s)=(sI-A)^{-1}(x_0+B(s)) $$
である。十分右の半平面で逆行列が存在する。

成分ごとに微分の公式を使えば $sX-x_0=AX+B$ となる。指数オーダーの評価はprop-lto-orderの行列の議論をそのまま適用できる。$\det(sI-A)$ は最高次係数1の多項式で根は有限個だから、その全ての実部より大きな半平面で逆行列が存在する。そこで左から掛ければ式を得る。$\square$

三角形の連立系

$x'=-x+z$、$z'=-2z$、$x(0)=0$、$z(0)=1$ を解く。$z=e^{-2t}$、$x=e^{-t}-e^{-2t}$ であり、$t=1$ では $(x,z)\approx(0.232544,0.135335)$。

途中の計算を開く

$(s+2)Z=1$、$(s+1)X=Z$。従って $X=1/((s+1)(s+2))=1/(s+1)-1/(s+2)$。両成分に逆変換を使えば答えを得る。

連立系で逆行列を使う場合、行列の掛ける順序を保つ。スカラーの $G/P$ の感覚で、行列を左右に自由に移してはいけない。2次の行列なら行列式と余因子で各成分を有理関数に直し、個別に逆変換できる。大きな系では全ての部分分数を作るより、行列の指数関数の計算法を使って $e^{At}$ を求める方が簡単な場合もある。
連立系の軌道は、各成分の式を同じ $t$ で評価して描く。$x$ と $z$ の別々の時刻の値を組み合わせても軌道にはならない。上の例では $z$ は単調に減衰し、$x$ は一度増えてからゼロへ戻る。$x$ の増加は不安定さを意味せず、初期に残っている $z$ の成分が $x$ を駆動することによる。

使える条件と反例

外す条件反例成り立たなくなること
定数係数$y'+ty=0$$\mathcal L[ty]=-Y'$ が現れ、代数方程式一つにはならない
線形性$y'+y^2=0$$\mathcal L[y^2]=Y^2$ という計算が使えない
入力の指数オーダー$y'=e^{t^2}$、$y(0)=0$解は存在するが全ての実数 $s$ で変換が発散する
初期値の取り込み$y'+2y=0$、$y(0)=1$ を $sY+2Y=0$ とする誤って $Y=0$ を得る

第1行では $sY-y(0)-Y'=0$ になる。これは $s$ を変数とする微分方程式なので、計算が単なる代数式にならない。Laplace変換そのものが使えないと断定しているのではなく、本記事の定数係数の手順の結論が変わることを示している。変数係数には別の工夫が必要である。
非線形の反例を確かめるため、単純な関数 $y=e^{-t}$ を考えると $\mathcal L[y^2]=1/(s+2)$、$Y^2=1/(s+1)^2$。$s=1$ ではそれぞれ $1/3$ と $1/4$ で一致しない。したがって非線形の項を変換後の未知量の同じ冪に置き換える操作は正しくない。なお $y'+y^2=0$ 自体は変数分離など別の方法で解ける。
増大する入力の例では、解は $y(t)=\int_0^te^{u^2}\,du$ であり、連続な入力なので通常の存在・一意性に問題はない。しかし $t\ge2$ で $y(t)\ge\int_{t-1}^te^{u^2}\,du\ge e^{(t-1)^2}$。任意の実数 $s$ に対して $e^{-st}y(t)$ は無限遠で増大するので変換が存在しない。解の存在と、Laplace変換でその解を扱えることは同じ条件ではない。
最後の行の正しい式は $(s+2)Y-1=0$ であり、$Y=1/(s+2)$、$y=e^{-2t}$ になる。誤ったゼロ関数は方程式だけを満たすが、初期値を満たさない。この反例は、変換後の右辺が入力だけでなく初期値も含むことを具体的に示す。

よくある誤りと補足

誤り正しくは理由
2階の微分で $sb_0$ を落とす$s^2Y-sb_0-b_1$微分の公式
初期値を変換したあとに別で合わせる初めから $R$ に入れるthm-lto-algebra
重根を単根として展開$t^ke^{at}$ を使う重根の展開公式
分母の根だけで振幅も決まる分子と初期値も必要例題の違い
ステップ入力で解を跳ばせる通常は低階導関数を連続に接続接続条件

減衰があるかどうかは、初期値だけでなく入力にも依存する。特性根の実部が負なら自由応答は減衰するが、一定入力なら解は通常ゼロでない定常値に近づく。指数的に増大する入力を与えれば、安定な分母でも応答が増大する。分母だけで判断できるのは系の自由応答やインパルス応答の性質であり、全ての入力に対する出力の値ではない。
時間側での定数変化法、変換側での $Y=(G+R)/P$、入力とインパルス応答の畳み込みは、同じ初期値問題を表す三つの形である。閉じた式を得るには部分分数が便利で、任意の入力を一つの式で表すには畳み込みが便利である。定常値だけが必要なら、条件を確認したうえで最終値定理を使えば、全ての逆変換を計算する必要がないこともある。

初期値を再び読み取る

求めた $Y$ から初期値を確かめるには、$s$ が大きいときの展開も使える。$y$ と必要な導関数が変換可能で原点で右連続なら、導関数の変換と 初期値定理と最終値定理 の初期値定理から
$$ \lim_{s\to\infty}sY(s)=y(0),\qquad \lim_{s\to\infty}s\{sY(s)-y(0)\}=y'(0) $$
を得る。この二つの式は初期値の項を落とした誤りを検出する助けになる。ただし、途中で跳びを持つ入力がある場合には、適用する導関数が原点で右連続であることを確認する。
例えば $y''+3y'+2y=e^{-3t}$、$y(0)=0$、$y'(0)=1$ では
$$ Y(s)=\frac{1+1/(s+3)}{(s+1)(s+2)}. $$
分子が $s\to\infty$ で $1$ に近づき、分母の最高次は $s^2$ なので、$sY\to0$ と $s^2Y\to1$ が得られる。これは指定した二つの初期値に一致する。入力の変換 $1/(s+3)$ だけを分子に残した誤った式では、$s^2Y\to0$ となり初速度を失う。
同じ確認を高階の問題にも適用できる。前の初期値を順に引いた
$$ s^jY(s)-\sum_{r=0}^{j-1}s^{j-1-r}y^{(r)}(0) $$
は $y^{(j)}$ の変換である。従ってこれに $s$ を掛けて大きな $s$ の極限を取れば $y^{(j)}(0)$ が得られる。形式的な漸近展開を使う場合も、背後にある導関数の変換可能性がこの読み取りを正当化している。

自由応答と零状態応答

式 $Y=(G+R)/P$ は二つに分けられる。$R/P$ が初期状態から来る自由応答、$G/P$ が入力から来る零状態応答である。入力を零にすれば前者だけが残り、初期値をすべて零にすれば後者だけが残る。線形性により両者の和が一般の解になる。
この分け方は、同じ入力を異なる初期状態から与える比較に便利である。二つの解に同じ入力 $g$ を与えれば、差は $g$ の項を含まない斉次方程式の解になる。その差が減衰するか増えるかは、入力でなく斉次部分の根による。入力が有界でも、自由応答に右半平面の根があれば解は増大しうる。
また入力の変換が特性多項式の根と同じ位置に極を持つと、積 $G/P$ の極の次数が上がる。それに対応して時間領域では $t$ の因子が現れる。これが共振の例で見た構造であり、線形微分方程式 の例「強制振動と共振」と同じ現象である。変換法は微分方程式の新しい種類の解を作るのではなく、既存の解の構造を極とその次数で整理している。
最終的な式を確かめる際には、初期値、微分方程式、入力の切替位置の三つを別々に見る。初期値が正しいだけでは、右辺の係数や周波数を誤った式を排除できない。微分方程式が正しいだけでも、異なる初期値の斉次解を足した可能性が残る。変換式との一致と時間領域での確認を組み合わせることで、計算の全体を確かめられる。
入力と初期状態を分けた表示では、同じ系の複数の問題を比較するときに、どの項が変化したかも追える。特性多項式を変えずに入力だけを変える計算と、係数そのものを変える計算を区別すると、分母を共用できる範囲が明確になる。

図で見る対応

図: 時間領域と変換領域の対応 図: 時間領域と変換領域の対応
左は強制入力 $e^{-3t}$ に対する二階系の解、右は連立系の軌道である。右図の横軸は $x$、縦軸は $z$ を表す。

関連項目

参考文献

Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する