この記事ではラプラス変換について説明します。数学的に深く掘り下げるのではなく、物理学や工学への応用を理解することを重視し、どのように使用するかに焦点を当てます。

前提知識としては高校の微積分、三角関数、指数関数、多項式方程式などで十分です。又複素数も重要ですが、今回の説明ではあえてそこまで触れません。

又おまけにPythonのsympyでのラプラス変換の実装も少し紹介します。

はじめに

ラプラスというと、平成育ちの人だったらやはりポケモンのことだと思うでしょう。私もそうですが、勿論全く別の話です。

（自分がラプラスに乗った画像をFLUX.1 Kontextで生成）

ラプラス変換はよく使われている数学の変換です。ラプラスというのはこの変換方法を考え出したフランスの数学者の名前です。

最近ラプラス変換を使う必要があるので、これは勉強するきっかけとなりました。元々大学時代に勉強したことがありますが、あの時あまり理解ができなく、実際に使う場面も全くなかったから、とっくに忘れたという状態でした。

フーリエ変換を使ったことがあるからラプラス変換も似たものだとわかっていますが、変換の意味はフーリエ変換ほど直感的ではないし、使い所もだいぶ違います。本当に理解できるまで随分大変でした。

せっかくだから自分が理解できたことをここに纏めておきたいです。私は本来専門は物理学なので、数学のガチ説明よりも使い方に重点を置きます。

だからまず解きたい物理学の問題から始め、これを目標として説明していきます。

解きたい運動方程式の問題

ラプラス変換の方法を学ぶ前に、まず問題を設定して、なぜそれを学ぶのかをイメージしやすくします。

ラプラス変換を応用できる問題は様々ですが、最もイメージしやすいのは物体の運動でしょう。これは高校の物理の授業でも見られます。

m\ddot{x}(t) = f(t)

ここで、$f(t)$は質量$m$の物体に作用する力、$ x(t)$は物体の位置であり、両方とも時間$t$の関数です。
ここでの点は時間に関する微分を表します。つまり、$\begin{align} \dot{x}(t) = \frac{d}{dt}x(t) \end{align}$および$\begin{align} \ddot{x}(t) = \frac{d^2}{dt^2}x(t) \end{align}$は時間の関数です。

よく見られる簡単な例は、バネを取り付けた車のシステムで、物体に作用する力はバネの伸びに比例します。

f_{バネ}(t) = - kx(t)

ここで、$k$はバネ定数、$x(t)$はバネの平衡点からの位置です。これで運動方程式はこのようになります。

m\ddot{x}(t) = - kx(t)

このような単純なシステムでは、周期を持った往復運動、つまり高校の授業でおなじみの単振動が発生します。

しかし、ここではより複雑にしてイメージを掴みやすくするため、ダンパー（減衰器）も含むシステムを考えます。

ダンパーは粘性抵抗を利用し、その力は速度に比例します。

f_{ダンパー}(t) = - c\dot{x}(t)

ここで、$c$は粘性減衰係数です。

バネの力とダンパーの力を合わせると、次のようになります。

f(t) = m\ddot{x}(t) = - kx(t) - c\dot{x}(t)

アニメーションを作成すると、以下のような運動になります。前後に繰り返し移動しますが、減衰によって移動距離が次第に短くなり、最終的には静止します。

上の式を整理し、項を移行します。

m\ddot{x}(t) + c\dot{x}(t) + kx(t) = 0

これは微分方程式の形、より具体的には線形常微分方程式と呼ばれるものになります。

このような微分方程式は、一見多項式方程式に似ていますが、実際に解くのが困難です。多項式方程式なら因数分解して解けますが、微分方程式はそれほど単純ではありません。

しかし、微分方程式を多項式方程式を解くのと同じように簡単に解くことを可能にする方法があります。それがラプラス変換です。

したがって、この方程式を解くことを目標としてこれからラプラス変換の紹介を始めます。

定義と計算方法

ラプラス変換は以下のように定義されます。

\mathcal{L}[x(t)] = \int_{0}^{\infty}x(t)e^{-st}dt

つまり、考察対象の値$x$（時間$t$の関数）に$e^{-st}$を乗じたものを、$t=0$から$t=\infty$まで積分したものです。

この$t$は、物理学の問題では通常時間を意味します（一般数学では常に時間である必要はありませんが）。

記号$\mathcal{L}$は、通常の変数Lと区別するために特別な書体で書かれたLです。$\mathcal{L}$の後ろの角括弧内は、変換したい関数です。

なお、数学の教科書では関数の変数として$f$が使われることが多いですが、ここでは位置の問題から始めたため変数として$x$を使用しています。混乱を避けるため、この記事では全て変数$x$を使って説明します。

式の中に新しい変数$s$が$e^{-st}$の中にあることがわかります。積分後、時間変数$t$は消え、残る変数は$s$だけになります。これは、ラプラス変換の結果が変数$s$の関数となることを意味するため、通常は次のように書きます。

X(s) = \mathcal{L}[x(t)]

この変数$s$はラプラス変数と呼ばれ、関数$X(s)$は$x(t)$のラプラス変換結果と呼ばれます。

通常、ラプラス変換された関数は大文字で表記され、変換前の元の関数は小文字で表記され、関連性を示します。今回は$s$と$X$です。ただし、これは絶対的な規則ではなく、変数は自由に使用できます。

簡単な変換の例を挙げてみましょう。例えば、次の関数のラプラス変換を試みます。

x(t) = e^{2t}

前述の定義に従って計算すると、以下のようになります。

\begin{align}
X(s) &= \mathcal{L}\left[e^{2t}\right] \\
&= \int_{0}^{\infty}e^{2t}e^{-st}dt \\
&= \int_{0}^{\infty}e^{(2-s)t}dt \\
&= \frac{1}{s-2} \\
\end{align}

したがって、$e^{2t}$のラプラス変換結果は$\frac{1}{s-2}$であることがわかります。

謎の変数 s について

変換式と変換例を見た後、多くの人が最も疑問に思うことは、突然$tの代わりに現れる変数$の代わりに現れる変数$s$についてでしょう。

この$s$は何なのでしょうか？私たちは時間$t$の関数である量を考えているのに、なぜ変換すると未知の変数に置き換わって現れるのでしょうか？

これは、かなり抽象的で解釈が難しい部分です。ラプラス変換の意味を真に理解するのは容易ではなく、想像力が必要です。

そして、率直に言ってしまえば、ラプラス変換を学ぶ主な目的は、その意味を理解することではなく、どのように使うかを理解することですよね。

物理学者やエンジニアにとって数学は、ある現象を説明するための道具であり、時にはその深い由来を完全に理解せずに使用しています。勿論、理解できればそれに越したことはありませんが、応用される数学的技法は非常に多く、すべてを理解するための時間は足りないのですよね。

しかし、そうは言っても、何の説明もなければイメージが湧かず、先に進めないかもしれません。そこで、ここでは物理学的な観点から大まかに説明します。

ラプラス変換は、問題の見方を変えるものです。元々は時間$t$の世界で値$x$を考察していたのを、ラプラス変数$s$の世界で値$X$を考察するように変換します。

世界は違うのです。だから考え方も見方も変わりますね。
~~異世界でラプラスは私たちの知っているラプラスじゃなくて、私も赤髪猫耳かもしれません。なんちゃって。~~

（Qwen-Image-Editで上の画像から編集したもの）

$x$が位置であれば、$X$も位置に関連する何かですが、直接的な意味で解釈できる位置の値ではありません。

$t$が時間であれば、$s$も時間に関連する何かですが、時間そのものではありません。$e^{-st}$の中で掛け合わされているから次元を考察すると、$s$は$t$の逆数の次元を持っていることはわかりますね。

時間の逆数の次元を持つものと言えば、最もよく知られているのは周波数です。

では、$s$は周波数なのでしょうか？実はそういは言い切れません。

イメージを掴むために、先ほど挙げたバネとダンパーが付いた車の問題に戻りましょう。この問題の解は、実際には次のような形になります。

x(t) = A e^{-\alpha t} \sin(\omega t + \theta)

これは、$\sin(\omega t + \theta)$という周期的に繰り返される部分と、$e^{-\alpha t}$という時間とともに減衰する部分から構成されています。グラフに描くと以下のようになります。

ここで、$\omega$は周波数（もっと正確は角周波数）であり、$\alpha$は、放射性崩壊などで見られる減衰定数（又は崩壊定数）と呼ばれるもので、何らかの量の時間に対する減衰（負の場合は増加）の速さを表します。

$\alpha$も$\omega$も時間に関する変化を表す値であり、次元は$s$と同じく時間の逆数です。更に見ていくと、周波数や時間減衰に関連するこれらの値が、$s$に加算又は減算という形で頻繁に現れることがわかります。

したがって、ラプラス変数$s$の世界では、$s$は時間的な変化の性質（周期的な変化を表す場合は周波数、減衰や増加を表す場合は減衰定数）を示していると解釈できます。

こんな感じの説明ですが、恐らくこれだけではまだ理解が難しいと思います。まだ本番の説明はれしていないから、イメージは湧かないでしょう。この記事全体を読み終えて全体像を掴んだ後で、もう一度読み返せば、より明確に理解できるかもしれません。

変換表

ラプラス変換は積分を実行することで行えますが、多少面倒です。一般に考察する問題の関数形は大体同じようなものが多いため、事前に計算された変換結果の表を参照して比較する方法が使えます。勿論、変換表を暗記する必要はありません。必要な時に表を見て変換すればよいのです。

ラプラス変換表は検索したらすぐ出るのでここでは載せませんが、リンクを貼っておきます。

線形性

次に、今回の主な目標である物体の運動問題を解くために必要な、ラプラス変換の重要な性質を紹介します。

ラプラス変換には線形性と呼ばれる性質があります。つまり、ある関数を定数倍したもののラプラス変換結果は、その関数のラプラス変換結果を定数倍したものです。

\mathcal{L}[ax(t)] = a\mathcal{L}[x(t)] = aX(s)

ここで、$a$は定数です。更に、2つの関数を足し合わせたもののラプラス変換は、それら2つの関数のラプラス変換の和に等しくなります。

\mathcal{L}[x_1(t)+x_2(t)] = \mathcal{L}[x_1(t)] + \mathcal{L}[x_2(t)] = X_1(s)+X_2(s)

勿論、各関数が定数倍されていても。

\begin{align}
      \mathcal{L}[a_1x_1(t)+a_2x_2(t)] &= a_1\mathcal{L}[x_1(t)] + a_2\mathcal{L}[x_2(t)] \\
      &= a_1X_1(s)+a_2X_2(s)
\end{align}

又は、関数が何個あっても良いので、一般的な形で書くと以下のようになります。

\begin{align}
\mathcal{L}[a_1x_1(t) + a_2x_2(t) + \cdots + a_nx_n(t)] & = a_1\mathcal{L}[x_1(t)] + a_2\mathcal{L}[x_2(t)] + \cdots + a_n\mathcal{L}[x_n(t)] \\
& = a_1X_1(s) + a_2X_2(s) + \cdots + a_nX_n(s)
\end{align}

使用例。

\begin{align}
      \mathcal{L}[2t + 3e^{-\alpha t}] &= 2\mathcal{L}[t] + 3\mathcal{L}[e^{-\alpha t}] \\
      &= \frac{2}{s^2} + \frac{3}{s+\alpha}
\end{align}

微分

ラプラス変換にはもう一つ非常に重要な性質があり、それは任意の関数の微分をラプラス変換すると次のようになることです。

\mathcal{L}[\dot{x}(t)] = sX(s) - x(0)

これは、元の時間$t$の世界での微分が、ラプラス変数の世界では$s$を掛けることになることを意味します。ただし、$x(0)$は初期時刻における$x$の値です。単なる定数なので、それほど複雑にはなりません。

更に、二次微分を求めてみると、次のようになります。

\mathcal{L}[\ddot{x}(t)] = s^2X(s) - sx(0) - \dot{x}(0)

そして、これを$n$次微分のラプラス変換の公式として書くこともできます。

\begin{align}
      \mathcal{L}\left[\frac{d^n}{dt^n}x(t)\right] = s^nX(s) - s^{n-1}x(0) - \cdots - \frac{d^{n-1}}{dt^{n-1}}x(0)
\end{align}

多くの場合、すべての初期条件がゼロである問題を考えることが多く、その場合は次のようになります。

\mathcal{L}\left[\frac{d^n}{dt^n}x(t)\right] = s^nX(s)

したがって、微分の形をした関数は、ラプラス変換後には単に$s$を掛けるだけであると簡単に考えることができます。これは非常に重要な点です。

これでなぜラプラス変換が微分方程式を解くのに役立つのか、そのイメージが湧いてきたのではないでしょうか。

それはつまり、ラプラス変換後、微分方程式が多項式方程式になるからです。

積分

微分と反対の操作が積分です。したがって、微分のラプラス変換が$s$を掛けることであるのに対し、積分のラプラス変換は$s$で割ることと等しくなります。

\begin{align}
\mathcal{L}\left[\int_{0}^{t}x(\tau)d\tau\right] = \frac{1}{s} X(s)
\end{align}

物体の運動問題への応用

ラプラス変換の方法と必要な性質を理解した後、先ほど挙げたバネとダンパーに繋がれた車の問題に戻りましょう。

まず、ラプラス変換を適用します。

\begin{align}
m\ddot{x}(t) + c\dot{x}(t) + kx(t) &= 0 \\
\mathcal{L}[m\ddot{x}(t) + c\dot{x}(t) + kx(t)] &= 0
\end{align}

線形性の性質を用いると、

m\mathcal{L}[\ddot{x}(t)] + c\mathcal{L}[\dot{x}(t)] + k\mathcal{L}[x(t)] = 0

次に、微分のラプラス変換の公式を用いると、

\begin{align}
m[s^2X(s) - sx(0) - \dot{x}(0)] + c[sX(s) - x(0)] + kX(s) &= 0 \\
(ms^2 + cs + k)X(s) - (ms+c)x(0) - m\dot{x}(0) &= 0
\end{align}

最終的に、

\begin{align}
X(s) &= \frac{mx(0)s+cx(0) + m\dot{x}(0)}{ms^2 + cs + k}
\end{align}

ここで、$x(0)$と$\dot{x}(0)$は初期位置と初速度です。分かりやすくするために、$x_0$と$v_0$と書き直すと、

\begin{align}
      X(s) &= \frac{mx_0s+(cx_0 + mv_0)}{ms^2 + cs + k}
\end{align}

これで、位置のラプラス変換後の表現である$X(s)$が得られました。しかし、これが最終目標ではありません。私たちが本当に必要としているのは、時間の関数としての位置$x(t)$ですよね。そのためには、逆ラプラス変換を行う必要があります。

逆ラプラス変換

逆ラプラス変換は次のように書けます。

x(t) = \mathcal{L}^{-1}[X(s)]

逆ラプラス変換にも計算の定義がありますが、複素数も絡んでやや複雑なので、ここでは詳細は省きます。逆ラプラス変換によく使われる方法は、ラプラス変換の時と同じく変換表と照らし合わせる方法です。ただ逆に見るだけ。

さて、バネとダンパーに繋がれた車の問題を再考します。解くためには、変換表に合わせて形を整え、逆変換できるようにする必要があります。

\begin{align}
      X(s) &= \frac{mx_0s+(cx_0 + mv_0)}{ms^2 + cs + k} \\
      &= x_0 \frac{s+(c/m + v_0/x_0)}{s^2 + (c/m)s + k/m} \\
      &= x_0 \frac{(s+c/2m) - c/2m + (c/m + v_0/x_0)}{(s+c/2m)^2 - (c/2m)^2 + k/m} \\
      &= x_0 \frac{s+c/2m}{(s+c/2m)^2 + \sqrt{k/m - (c/2m)^2}^2} \\
      &+ x_0 \frac{c/2m + v_0/x_0}{\sqrt{k/m - (c/2m)^2}} \frac{\sqrt{k/m - (c/2m)^2}}{(s+c/2m)^2 + \sqrt{k/m - (c/2m)^2}^2}
\end{align}

次のように置きます。

\alpha = c/2m

\omega = \sqrt{k/m - (c/2m)^2}

そして、次のように書き直せます。

\begin{align}
      X(s) &= x_0 \frac{s+\alpha}{(s+\alpha)^2 + \omega^2} + x_0 \frac{\alpha + v_0/x_0}{\omega} \frac{\omega}{(s+\alpha)^2 + \omega^2}
\end{align}

ここで変換表を見ると、

\begin{align}
      \mathcal{L}^{-1}\left[ \frac{s+\alpha}{(s+\alpha)^2+\omega^2} \right] = e^{-\alpha t} \cos(\omega t) \\
      \mathcal{L}^{-1}\left[ \frac{\omega}{(s+\alpha)^2+\omega^2} \right] = e^{-\alpha t} \sin(\omega t)
\end{align}

とあります。したがって、逆ラプラス変換を行い、次の結果を得ます。

\begin{align}
      x(t) &= x_0e^{-\alpha t} \cos(\omega t) + x_0\frac{\alpha + v_0/x_0}{\omega}e^{-\alpha t} \sin(\omega t) \\
      &= x_0e^{-\alpha t} \left(\cos(\omega t) + \frac{\alpha + v_0/x_0}{\omega}\sin(\omega t) \right)
\end{align}

ここで三角関数の性質（ここに参照）を用いると、次の形に書くことができます。

x(t) = A e^{-\alpha t} \sin(\omega t + \theta)

ただし、

\begin{align}
      A &= x_0 \sqrt{1 + \left( \frac{\alpha + v_0/x_0}{\omega} \right)^2} \\
      \theta &=  \text{atan2}\left(\omega, \alpha + \frac{v_0}{x_0}\right)
\end{align}

得られた結果は、上記で示したグラフとアニメーションの通りとなります。$e^{-\alpha t}$という指数関数的減衰部分と、$\sin(\omega t + \theta)$という周期的振動部分があり、したがって、減衰振動となり、最終的には停止します。

ただし、ここで$\omega = \sqrt{k/m - (c/2m)^2}$は、バネ定数が非常に小さい$k \le c^2/4m$の場合、0又は負の値になる可能性があります。この場合、臨界減衰又は過減衰が発生し、周期的な振動は起こりません。

更に、ダンパーがない場合、つまり$c = 0$の場合、$\alpha = 0$となり、指数関数部分が消え、

x(t) = A \sin(\omega t + \theta)

となります。これは永遠に振動し続ける単振動です。

この例から、ラプラス変換を用いて物理学の問題のような微分方程式を解くことができることがわかります。ラプラス変換を行い、ラプラス変数の世界で形を整えたり様々な操作を行い、終わったら逆ラプラス変換で戻します。

しかし、問題が複雑になればなるほど、逆ラプラス変換は難しくなり、常に戻せるとは限りません。

ただし、逆変換が常に必須というわけではありません。ラプラス変数の世界での値は、私たちが慣れ親しんだ時間の世界ほど直感的に理解しやすいものではありませんが、そのままの形で解釈できる場面もあります。これはフーリエ変換での周波数解析に似ています。

最終値の定理

ラプラス変換のもう一つの重要でよく使われる性質として、最終値の定理があります。次のように説明されます。

x(\infty) = \lim_{s \to 0}sX(s)

つまり、時間$t$が無限大に経過したときの値は、$s$を 0 に置き換えるだけで簡単に求めることができ、逆ラプラス変換を行う必要はありません。

ただし、条件として、$x_0$の値が時間の経過とともにある値に収束しなければなりません。発散する値の場合は、そもそも値を求めることができません。

使用例として、上記で説明したバネとダンパーに繋がれた車の問題に適用し、非常に長い時間が経過した後の車の位置を求めてみましょう。

\begin{align}
      x(\infty) &= \lim_{s \to 0}s\frac{mx_0s+(cx_0 + mv_0)}{ms^2 + cs + k} \\
      &= \lim_{s \to 0}\frac{mx_0s^2+(cx_0 + mv_0)s}{ms^2 + cs + k} \\
      &= 0
\end{align}

初期位置がどこであれ、初速度がどれほどであれ、最終的に車は位置 0 で停止することがわかります。

更に理解を深めるために、車に一定の力$f_0$を加える例を追加します。これは、バネとダンパーからの力に加えて、車に一定の補助力を加えることを意味し、運動方程式は次のようになります。

m\ddot{x}(t) + c\dot{x}(t) + kx(t) - f_0 = 0

これをラプラス変換します。

\begin{align}
m\mathcal{L}[\ddot{x}(t)] + c\mathcal{L}[\dot{x}(t)] + k\mathcal{L}[x(t)] - \mathcal{L}[f_0] &= 0 \\
m[s^2X(s) - sx(0) - \dot{x}(0)] + c[sX(s) - x(0)] + kX(s) - \frac{f_0}{s} &= 0 \\
(ms^2 + cs + k)X(s) - (ms+c)x(0) - m\dot{x}(0) - \frac{f_0}{s} &= 0
\end{align}

そして最終的に、

\begin{align}
X(s) &= \frac{mx(0)s+cx(0) + m\dot{x}(0) + f_0/s}{ms^2 + cs + k}
\end{align}

最終値の定理を用いて無限時間後の値を求めます。

\begin{align}
      x(\infty) &= \lim_{s \to 0}sX(s) \\
      &= \lim_{s \to 0}s\frac{mx_0s+cx_0 + mv_0 + f_0/s}{ms^2 + cs + k} \\
      &= \lim_{s \to 0}\frac{mx_0s^2+cx_0s + mv_0s + f_0}{ms^2 + cs + k} \\
      &= \frac{f_0}{k}
\end{align}

これは、ある力を加えたときに得られるべきバネの伸びの理論上の長さに一致します。引き始めにはまず振動があるが、最終的には$f_0/k$の位置で静止します。

ただし注意点として、ここでは$c=0$を代入しても$f_0/k$が得られますが、これは、減衰力がなければシステムは永遠に振動し続け、いかなる値にも収束しないという事実に反します。したがって、収束しない場合、この最終値の定理は成り立たず、使用できないことを示しています。よって使用前には、値が収束することを確認する必要があります。

Pythonのsympyでラプラス変換

最後にPythonでラプラス変換と逆ラプラス変換をする便利な使い方を紹介したいと思います。

よく使われるのはsympyというライブラリです。sympyは様々な数式処理をすることができます。詳しい使い方は割愛しますが、ここでは特にラプラス変換に関する機能だけ紹介します。

まずはラプラス変換の結果を求める例です。

import sympy as sp

a,s,t = sp.symbols('a,s,t') # 変数のシンボルの定義

x = sp.sin(a*t) # 変換したい関数。
X = sp.laplace_transform(x,t,s)[0] # 変換実行
print(X)

laplace_transform関数からの返り値は、変換結果関数の他に変数条件の情報も含まれたtupleとなりますが、ここでは必要ないので[0]を付けて変換結果関数Xだけを取ります。

結果。

a/(a**2 + s**2)

そして逆ラプラス変換も簡単にできます。

inv = sp.inverse_laplace_transform(X,s,t)
print(inv)

a*sin(t*sqrt(a**2))*Heaviside(t)/sqrt(a**2)

ここではHeaviside(t)が現れます。これはヘヴィサイドの階段関数のことです。ステップ関数とも呼ばれます。つまり$t \ge 0$以外の値は全部0となります。

元々ラプラス変換では$t=0$から$\infty$までしか部分計算に使わないのでたとえ元の関数が$t<0$に値があっても意味がありません。だから逆ラプラス変換の結果も厳密には$t<0$の値があるはずがないでしょう。

ただし実際に物理学でラプラス変換を使うのは$t \ge 0$しか興味ない場合が殆どなので、これを意識する必要があまりないと思いますね。

latexを使ったらlatexコードで表示することができます。

print(sp.latex(inv))

結果

\frac{a \sin{\left(t \sqrt{a^{2}} \right)} \theta\left(t\right)}{\sqrt{a^{2}}}

ここでヘヴィサイドの階段関数は$\theta$として表示されます。ところでたとえ$\theta(t)$がなくてもまだ元の$sin(at)$には戻っていませんね。$\sqrt{a^{2}}$が$a$にしたら元に戻るはずですが、なぜそうしないでしょう？

それは$a$が複素数で虚数が含まれる可能性があるからです。今回は複素数なしで説明するつもりなのでこれについては割愛しますが、このようなガチ数学で考えると慎重に取り組む必要があります。今回の記事ではそんな厳密な条件の説明をだいぶ省けています。

シンボル定義の時にちゃんとpositive=True（正の実数のみ）などの条件を付けることによって、思った通りの簡単な結果が得られます。

a,s,t = sp.symbols('a,s,t',positive=True)
x = sp.sin(a*t)
X = sp.laplace_transform(x,t,s)[0]
inv = sp.inverse_laplace_transform(X,s,t)
print(inv) # sin(a*t)

バネとダンパーが付いた車の問題でも試してみましょう。

s,t,m,c,k = sp.symbols('s,t,m,c,k',positive=True) # 正の実数
x_0,v_0 = sp.symbols('x_0,v_0',real=True) # 実数
X = (m*x_0*s + c*x_0 + m*v_0)/(m*s**2 + c*s + k)
inv = sp.inverse_laplace_transform(X,s,t)
print(sp.latex(inv))

\begin{align}
\frac{2 m \left(- \frac{c x_{0}}{2 m} + \frac{c x_{0} + m v_{0}}{m}\right) e^{- \frac{c t}{2 m}} \sin{\left(\frac{t \sqrt{- c^{2} + 4 k m}}{2 m} \right)}}{\sqrt{- c^{2} + 4 k m}} + x_{0} e^{- \frac{c t}{2 m}} \cos{\left(\frac{t \sqrt{- \frac{c^{2}}{4} + k m}}{m} \right)}
\end{align}

書き方は違いますが、整理したら上述の解と一致していますね。

最後に力$f_0$を加えた場合も解いてみましょう。

s,t,m,c,k = sp.symbols('s,t,m,c,k',positive=True)
x_0,v_0,f_0 = sp.symbols('x_0,v_0,f_0',real=True)
X = (m*x_0*s + c*x_0 + m*v_0 + f_0/s)/(m*s**2 + c*s + k)
inv = sp.inverse_laplace_transform(X,s,t)
print(sp.latex(inv))

\begin{align}
\frac{f_{0}}{k} + \frac{2 m \left(- \frac{c \left(- f_{0} m + k m x_{0}\right)}{2 k m^{2}} + \frac{- c f_{0} + c k x_{0} + k m v_{0}}{k m}\right) e^{- \frac{c t}{2 m}} \sin{\left(\frac{t \sqrt{- c^{2} + 4 k m}}{2 m} \right)}}{\sqrt{- c^{2} + 4 k m}} + \frac{\left(- f_{0} m + k m x_{0}\right) e^{- \frac{c t}{2 m}} \cos{\left(\frac{t \sqrt{- \frac{c^{2}}{4} + k m}}{m} \right)}}{k m}
\end{align}

.subsメソッドで値を代入することができます。

sub_lis = [
    (x_0,0),
    (v_0,0),
    (f_0,200),
    (k,20),
    (m,500),
    (c,40),
]
xt = inv.subs(sub_lis)
print(sp.latex(xt))

\begin{align}
10 - \frac{5 \sqrt{6} e^{- \frac{t}{25}} \sin{\left(\frac{2 \sqrt{6} t}{25} \right)}}{6} - 10 e^{- \frac{t}{25}} \cos{\left(\frac{2 \sqrt{6} t}{25} \right)}
\end{align}

そしてグラフをプロットする機能もあります。

sp.plot(xt,(t,0,100))

ただしその中身はmatplotlibを使っています。もっと自由に描きたい場合lambdify関数を使ってnumpy配列で計算できるようにして直接matplotlibでプロットした方がいいかもしれません。

xfunc = sp.lambdify(t,xt)

import numpy as np
import matplotlib.pyplot as plt

t_graph = np.linspace(0,100,1001)
plt.figure(figsize=[6,4],dpi=100,facecolor='#ffeedd')
plt.ylabel(r'$x(t)$',size=14)
plt.xlabel(r'$t$',size=14)
plt.plot(t_graph,xfunc(t_graph),'r')
plt.grid(ls=':')
plt.show()

なお、このグラフを見たら最終値の定理で説明した通り、$f_0/k = 200/20 = 10$に収束するとわかりますね。

纏め

ラプラス変換は、物理問題の解決を容易にするために使用できる
ラプラス変換は、時間 t の世界からラプラス変数 s の世界への変換である
時間 t の世界での微分は、ラプラス変数の世界では s を掛けることに対応する
時間 t の世界での積分は、ラプラス変数の世界では s で割ることに対応する
これによって本来微分方程式は多項式方程式として解く問題となって、楽になる
ラプラス変数 s の次元は時間 t の逆数であり、周波数または減衰定数と見なすことができる
ラプラス変数 s の世界から、慣れ親しんだ時間 t の世界の値に、逆ラプラス変換を行うことで戻すことができる
どのように解釈するかを理解していれば、ラプラス変数の世界の値を戻さずに解析できる場合もある
逆ラプラス変換は、直接計算するよりも、変換表に合うように式を整理して行うことが多い
収束するシステムの無限時間後の値のみに関心がある場合は、最終値の定理を用いて逆ラプラス変換せずに求めることができる

ラプラス変換についてはまだ説明しきれない詳細が山程あります。ここで説明したのは初歩的な部分だけで、ラプラス変換が何のためにあり、どのように使うのかのイメージを持ってもらうためです。興味があれば、さらに詳しく学んでいくことができます。

もっと読む

制御工学分野におけるラプラス変換の活躍はこの記事に書いてあります。

又、ラプラス変換は基本的に連続時間に使うものですが、ディジタル信号処理などでは「離散時間のラプラス変換」と呼ばれるZ変換が使われます。このZ変換についても記事を書きました。

参考

ラプラス変換の基本

動画

sympy