トップ > ライブラリ > 数値解析 > ラプラス変換

数値解析：ラプラス変換

　現在工事中です。
　ラプラス変換は、工学の各分野で偏微分方程式を解くために利用される汎用性のある、また応用性のある解析的な応用数学技術のひとつである。 純数学的な知識として身につける必要はまったくなく、単純に応用数学の世界で有効に利用できるテクニックのひとつと考え、 応用・利用できる程度の知識を身につければよい。
　ラプラス変換は、いわゆる数学的に実空間を写像空間に変換し、数学的な取り扱いを簡便にする方法の一つで、写像空間での処理の後、 逆変換することで、実空間での現象に戻すという一連の操作である。ある種の偏微分方程式の解析解を導出することも可能である。
　以下、ラプラス変換の基礎知識を紹介するとともに、工学への応用例を簡単に紹介する。

先頭に戻る

【注】数式表示にはMathJaxを利用しています。IE8以下では表示が遅くなる可能性があります。FireFox などIE8以外のブラウザを利用下さい。

ラプラス変換とは

　工学の分野ではしばしばラプラス変換を使って、非定常の伝熱方程式（Fourierの熱伝導式）や非定常の拡散方程式（Fickの拡散方程式） を解析的に、非常にスマートな形で解くことをしている。
　条件によっては解析的に解くことができない場合もあるが、定数係数を持つ、非定常１次元の熱伝導方程式（または拡散方程式）では 解析解が、バイブルといわれる、Carslaw-Jaegerの図書3)やCrankの図書4)に数多く紹介されている。

　筆者が学生の頃、コンピュータがまだ世の中に出現していない頃の古い文献で、ラプラス変換を用い偏微分方程式の解析解が記載されていると、スマートな導出方法に憧れたものです。 これが起因となって、独学でラプラス変換の手法を身につけました。

　さて、ラプラス変換の数学的な定義はつぎのとおりです。
　関数\(f(t)\ \)が\(t>0\ \)で定義されているとき、複素数\(s\ \)について、積分 \[ \begin{align*} F(s)=\int^\infty_0 e^{-st}f(t)dt \equiv \lim_{T \to \infty} \int^T_0 e^{-st}f(t)dt \tag{1} \end{align*} \] が存在すれば、これを\(f(t)\ \)のラプラス変換といい、 \[ \begin{align*} \mathcal{L} f(t)=F(s) \end{align*} \] で表す。記号\(\mathcal{L}\ \)は実変数\(t\ \)の関数\(f(t)\)に対して複素変数\(s\ \)の関数\(F(s)\)を対応させる変換を表します。
　ラプラス変換は実変数\(t\ \)の空間から、複素数\(s\ \)の空間に写像したことになる。一方、ラプラスの逆変換は複素数空間からもとの実変数空間 に写像することをいう。

例題１：ラプラス変換

　\(f(t)=c,\quad t>0\)のラプラス変換は \[ \begin{align*} \mathcal{L}c=\frac{c}{s}, \quad \Re(s)>0 \tag{2} \end{align*} \] である。ここで\(c\ \)は定数である。また\( \Re(s)\)は\(s\ \)の実部である。
【解答】　\(f(t)=c \)を定義式(1)に代入すると \[ \begin{align*} \mathcal{L}c & =\int^\infty_0 e^{-st}cdt=\Big[ \frac{c}{-s}e^{-st}\Big]^\infty_0 \cr &=\frac{c}{-s} \Big[e^{-s \times \infty}-e^0 \Big] \cr &=\frac{c}{-s} [0-1] = \frac{c}{s} \end{align*} \] ここで\(s\ \)の実数は正であるから、\( \lim_{t \to \infty}e^{-st}=0\)を利用した。

　いくつかの基本的な関数のラプラス変換を表１にまとめた。

表１：基本的な関数のラプラス変換
\[f(t),t>0\]\[F(s)=\mathcal{L}f(t)\]
\[1\]\[\frac{1}{s}\]
\[t\]\[\frac{1}{s^2}\]
\[\frac{t^n}{n!}\]\[\frac{1}{s^{n+1}}\]
\[e^{at}\]\[\frac{1}{s-a}\]
\[\sin \omega t\]\[\frac{\omega}{s^2+\omega^2}\]
\[\cos \omega t\]\[\frac{s}{s^2+\omega^2}\]

　これら関数以外にも、「数学公式II」（岩波全書）1)に多数の変換式がリストアップされている。

先頭に戻る

ラプラス変換の性質

　ラプラス変換を一般の関数\(f(t)\)に適用する際に、便利な性質を以下に示す。偏微分方程式を解析的に解く際に有用な関係式である。

　いま、\(c_1,c_2,a\ \)を定数とし、\(\mathcal{L}f(t)=F(s)\)とすれば、次の関係が成立する。 また\(f'(t)、f^n(t)\)は、１階微分およびｎ階微分を示す。 \[ \begin{align*} \mathcal{L}[c_1f_1(t)+c_2f_2(t)] & =c_1\mathcal{L}f_1(t)+c_2\mathcal{L}f_2(t) \tag{3} \cr \mathcal{L}f(at) & =\frac{1}{a}F(\frac{s}{a}) \tag{4} \cr \mathcal{L}[e^{at}f(t)] & =F(s-a) \tag{5} \cr \mathcal{L}f'(t) & =sF(s)-f(+0) \tag{6} \cr \mathcal{L}f^{(n)} & =s^nF(s)-s^{n-1}f(+0)-s^{n-2}f'(+0)-\cdots \cr & -sf^{n-2}(+0)-f^{n-1}(+0) \tag{7} \cr \mathcal{L}\Big[\int^t_0f(t)dt \Big] & = \frac{F(s)}{s} \tag{8} \cr \mathcal{L}[tf(t)] & =-\frac{d}{ds}F(s) \tag{9} \cr \mathcal{L}\Big[ \frac{f(t)}{t} \Big] & =\int^\infty_s F(s)ds \tag{10} \cr \end{align*} \] 　この中で、特に(6)式の一時微分が、変換後の関数\(F(s)\)と変換前の初期値\(f(+0)\)の線形結合 で表されることが特徴的である。この性質により、非定常問題の偏微分方程式が、常微分方程式に変換できることとなる。

先頭に戻る

ラプラス逆変換

　ラプラス変換とは逆に、複素数空間の関数\(F(s)\)から実空間の関数\(f(t)\)を求めることを逆変換といい、 \[ \begin{align*} \mathcal{L}^{-1}F(s)=f(t) \tag{11} \end{align*} \] と表す。
　先に示した表１のような変換表、「数学公式II」1)などの図書による変換表を元に逆変換することができる。

例題２：ラプラス逆変換

次のラプラス逆変換 \[ \mathcal{L}^{-1} \frac{s+2}{(s-1)^2s^3} \] を求めよ。
【解答】　例題の有理関数を次のように部分分数に分解すると便利である。 \[ \begin{align*} \frac{s+2}{(s-1)^2s^3}=\frac{3}{(s-1)^2}-\frac{8}{s-1}+\frac{8}{s}+\frac{5}{s^2}+\frac{2}{s^3} \cr \end{align*} \] 表１の変換表を用いると、 \[ \begin{align*} \mathcal{L}^{-1} \frac{s+2}{(s-1)^2s^3} & =3\mathcal{L}^{-1}\frac{1}{(s-1)^2} -8\mathcal{L}^{-1} \frac{1}{s-1}+ 8\mathcal{L}^{-1} \frac{1}{s} +5\mathcal{L}^{-1} \frac{1}{s^2}+2\mathcal{L}^{-1} \frac{1}{s^3} \cr & = 3te^t-8et+8+5t+t^2 \cr & = (3t-8)e^t+(t^2+5t+8) \end{align*} \] となる。

先頭に戻る

留数定理

　ここでは複素変数についての数学的な定理を述べる。詳細は専門の図書を参照されたい。

　複素関数\(F(s)\)のラプラス逆変換は、k個のすべての極\(s_n\ (n=1,2,\cdots,k)\ \)の留数の総和で表される。 \[ \begin{align*} \mathcal{L}^{-1}F(s) = \sum^k_{n=1} {\rm Res}[e^{st}F(s),s_n] \tag{12} \end{align*} \] ここで、\({\rm Res}[f(s),a] \)は複素関数\(f(s)\)の点\(a\)における留数を示す。すなわち、複素関数\(F(s)\)のラプラス逆変換 は\(e^{st}F(s)\)の留数の総和で表される。これを留数定理という。

　留数の性質として、点\(a\ \)が\(f(z)\)の１位の極であれば、点\(a\ \)における留数は \[ \begin{align*} {\rm Res}[f(z),a]=\lim_{z \to a}(z-a)f(z) \tag{13} \end{align*} \] である。また、点\(a\ \)が\(f(z)\)の\(k\ \)位の極であれば、点\(a\ \)における留数は \[ \begin{align*} {\rm Res}[f(z),a]=\frac{1}{(k-1)!}\lim_{z \to a}\frac{d^{k-1}}{dz^{k-1}} [(z-a)^kf(z)] \tag{14} \end{align*} \] である。ここで、点\(a\ \)が\(f(z)\)のk位の極とは、\((z-a)^kf(z)\)は点\(a\ \)で正則である。 あるいは、\(f(z)\)がつぎの形 \[ f(z)=\frac{g(z)}{(z-a)^k} \tag{15} \] に書けるとき、点\(a\ \)は\(f(z)\)のk位の極であるという。

例題３：留数の例

複素関数\(f(z)\)が次式 \[ f(z)=\frac{e^z}{(z-1)(z+3)^2} \] のとき、\({\rm Res}[f(z),-3]\)を求めよ。
【解答】　-3は\(f(z)\)の2位の極である。(14)式から、 \[ \begin{align*} {\rm Res}[f(z),a] & =\frac{1}{(2-1)!}\lim_{z \to -3}\frac{d^{2-1}}{dz^{2-1}} [(z+3)^kf(z)] \cr & = \lim_{z \to -3} \frac{d}{dz}\frac{e^z}{(z-1)} \cr & = \lim_{z \ to -3} \frac{(z-1)e^z-e^z}{(z-1)^2} =-\frac{5}{16}e^{-3} \cr \end{align*} \]

先頭に戻る

Heavisideの展開定理

　ラプラスの逆変換するとき、複素関数\(F(z)\)の分子分母が、次のようにが有理関数\(P(z)/Q(z)\)と表され、以下のように部分分数分解できると、逆変換が容易にできる。

　いま、\(P(z),Q(z)\)は\(z\ \)の実数係数多項式で、\(Q(z)\)の次数＞\(P(z)\)の次数とする。また \[ \begin{align*} Q(z)= (z-z_1)^{n_1}(z-z_2)^{n_2}\cdots (z-z_k)^{n_k} \tag{16} \end{align*} \] と因数分解できるとする。\(z_1,z_2,\cdots,z_k\)は相異なる複素数とする。
　先の留数のところで述べたように、\(z_1\)は\(n_1\)位の極、\(z_2\)は\(n_2\)位の極、･･･に相当する。
　定数\(a_{ij}, 1 \ge i \ge k, 1 \ge j \ge n_j\)をうまく選べば \[ \begin{align*} \frac{P(z)}{Q(z)}= \sum^k_{i=1} \sum^{n_i}_{j=1} \frac{a_{ij}}{(z-z_i)^j} \tag{17} \end{align*} \] と部分分数和で表すことができる。

　ここで、定数\(a_{ij}\ \)は \[ \begin{align*} a_{ij} = \frac{1}{(n_i-j)!} \lim_{z \to z_i} \frac {d^{n_i-j}}{dz^{n_i-j}}\Big\{(z-z_i)^{n_i} \frac{P(z)}{Q(z)}\Big\} \tag{18} \end{align*} \] で計算できる。以上がHeavisideの展開定理といわれる。

　特に\(Q(z)\)の極の\(z_1,z_2,\cdots,z_k\)がすべて1位のとき、すなわち(16)式が \[ \begin{align*} Q(z)= (z-z_1)(z-z_2)\cdots (z-z_k) \tag{16'} \end{align*} \] で表されると、留数定理の(12),(13)式から \[ \begin{align*} \mathcal{L}^{-1}F(s) = \mathcal{L}^{-1}\frac{P(s)}{Q(s)} = \sum_{n=1}^k\frac{P(z_n)}{Q'(z_n)}e^{z_nt} \tag{19} \end{align*} \] と逆変換される。ここで\(Q'(s)\)は\(dQ(s)/ds\)を意味する。

先頭に戻る

微分方程式の解析解

　ラプラス変換を使うとある種の微分方程式（定数係数をもつ微分方程式など）の解析解を求めることが可能となる。

常微分方程式

　ラプラス変換を利用して常微分方程式を解くことを、例を通して考える。

例題４：常微分方程式の解法

次の常微分方程式 \[ \begin{align*} \frac{d^2y}{dt^2}+3\frac{dy}{dt}+2y=0 \tag{20} \end{align*} \] を、初期条件\(y(0)=y_0,y'(0)=v_0\ \)のもとで解け。ただし\(t \ge 0\ \)とする。
【解答】　\(\mathcal{L}y(t)=Y(s)\ \)とする。両辺をラプラス変換し、(6),(7)式を適用すると \[ \begin{align*} & \mathcal{L}y''+3\mathcal{L}y'+2\mathcal{L}y = 0 \tag{21a} \cr & \{s^2Y(s)-sy(0)-y'(0)\}+3\{sY(s)-y(0)\}+2Y(s) = 0 \tag{21b} \cr \end{align*} \] 初期条件を代入し、\(Y(s)\)について解くと \[ \begin{align*} Y(s)=\frac{y_0s+3y_0+v_0}{s^2+3s+2}=\frac{2y_0+v_0}{s+1}-\frac{y_0+v_0}{s+2} \tag{22} \end{align*} \] となり、部分分数に分解される。これを、逆変換すると微分方程式の解として \[ \begin{align*} y(t)=\mathcal{L}^{-1}Y(s)=(2y_0+v_0)e^{-t}-(y_0+v_0)e^{-2t} \tag{23} \end{align*} \] が得られる。

先頭に戻る

偏微分方程式の解析例

　次は古い文献などで偏微分方程式の解析解が結果として導出されていることがあり、その導出までの経緯は大概省略されていたりする。 ここでは、典型的な偏微分方程式で解析解の導出方法をいくつか例を交えて紹介する。

　また工学的な熱伝導の問題や拡散の問題については、バイブルといわれるCarslaw and Jaeger3)の図書、Crank4)の図書 に、偏微分方程式の解析解が数多く紹介されている。ただし、これら図書では途中の導出を省略していることが多く、導出した結果として得られる式だけが記載されていることが多い。

先頭に戻る

半無限固体中の非定常伝熱（直交座標系）

　まずはじめに、直交座標系の非定常問題を取り上げる。

例第５：半無限固体中の非定常伝熱5)

　\(z=0\ \)から\(z=\infty\)まである半無限の固体が最初\(T=T_i\)に保たれ、時間\(t \gt 0\)で \(z=0\)の表面温度が\(T_w\)となりそのまま維持されるとき、固体中の温度分布の時間変化を求めよ。
　z方向の熱伝導の基礎式(5-1)式、および初期条件(5-2)式、境界条件(5-3)式は以下のとおりとする。熱伝導度\(\lambda\)は定数とす。 \[ \begin{align*} \frac{\partial T}{\partial t}=\lambda \frac{\partial^2 T}{\partial z^2} \tag{5-1} \end{align*} \] ここで、\(\lambda\)は固体の熱伝導度、\(T\)は温度、\(t\)は時間、\(z\)は表面からの距離とする。 \[ \begin{align*} {\rm I.C.} \quad t=0, z \ge 0 : T=T_i \tag{5-2} \end{align*} \] \[ \begin{align*} \left. \eqalign{ {\rm B.C.} \quad t>0,z=0 : & \quad T=T_w \cr z=\infty : & \quad T=T_i \cr } \right\} \tag{5-3} \end{align*} \] ここで、\(T_i\)は固体の初期温度、\(T_w\)は固体の表面温度で、ともに\(t,z\)には依存しない定数とする。

　例題５の解答例は、こちらラプラス変換・例題５に詳述しています。

先頭に戻る

粒子充填層の滞留時間分布（直交座標系）

　次の例題６も、直交座標系で非定常問題であるが、こちらはラプラス変換を利用した解法であり、 式導出が難しい。

例題６：触媒充填層の混合拡散モデルL1342),6)

　流体が\(x\)方向に平均流速\(u\)で、軸方向混合拡散係数\(E_a\)で流れている。ある位置\(x\)ある時間\(t\)での成分濃度\(c\)の物質収支 式は次式(6-1)で表される。ここで\(\psi(c)\)は反応項を示す。
　初期条件(6-2)式、境界条件(6-3)式とするとき、このモデルの滞留時間分布関数\(E(\phi)\)を求めることを考える。ただし、 反応項\(psi(c)=0\)とし、\(\delta(x)\)はデルタ関数（ステップ応答を求めるため、デルタ関数を初期値とす）、\(L\)は反応器長さとする \[ \begin{align*} \frac{\partial c}{\partial t}+u\frac{\partial c}{\partial x}= E_a \frac{\partial^2c}{\partial x^2}+\psi(c) \tag{6-1} \cr \end{align*} \] \[ \begin{align*} \left. \eqalign{ {\rm I.C.} \quad t<0 : & \quad c=0 \cr t=0 : & \quad c=\delta(x) \cr } \right\} \tag{6-2} \end{align*} \] \[ \begin{align*} \left. \eqalign{ {\rm B.C.} \quad x=0 : & -E_a\Bigl( \frac{dc}{dx}\Bigr)_{x=0+}=u(c-c_{x=0+}) \cr x=L : & -E_a\Bigl( \frac{dc}{dx}\Bigr)=0 \cr } \right\} \tag{6-3} \end{align*} \]

　求める滞留時間分布関数\(E(\phi)\)は次式となる。 \[ \begin{align*} E(\phi)=2 \sum^\infty_{n=1}\frac{(-1)^{n+1}\mu_n^2 \exp({\rm Pe}/2)}{({\rm Pe}/2)^2+{\rm Pe}+\mu_n^2} \exp \Bigl\{ \Bigl( -\frac{({\rm Pe}/2)^2+\mu_n^2}{{\rm Pe}} \Bigr) \phi \Bigr\} \tag{6-4} \end{align*} \] ここで、\(\phi=tu/L\)で表される無次元化された時間、\({\rm Pe}=uL/E_a\)で定義されるPecret数、\(\mu_n\)は次式のn番目の正根である。 \[ \begin{align*} \tan \mu_n & =\frac{{\rm Pe} \mu_n}{\mu_n^2-({\rm Pe}/2)^2} \tag{6-5} \end{align*} \] 以上、ラプラス変換を用いて証明せよ。

　例題６の解答例は、こちらラプラス変換・例題６に詳述しています。

先頭に戻る

着目成分の非定常拡散方程式（円筒座標系）

　次の例題７は、円筒座標系の非定常問題である。ラプラス変換を利用した解法であり、式導出が難しい。 また、円筒座標系で頻出するベッセル関数を取り扱う。

例７：円柱からの拡散方程式の解法4)

　時刻ゼロで均一の成分濃度\(C_2\)に保たれた、無限長さをもった円柱を考える。その表面から着目成分が蒸発し、時間とともに成分濃度\(C\)が変化する。 円柱自体の形状の時間変化はないものとし、円柱内の濃度分布の時間変化を求める。
　このときの基礎式は、円筒座標系を用い、次の(7-1)式で表される。 \[ \begin{align*} \frac{\partial C}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\Bigl( r D \frac{\partial C}{\partial r} \Bigr) \tag{7-1} \end{align*} \] ここで、\(D\ \)は成分の拡散係数で時間・位置には依存しないとする。
　初期条件は、次のように、任意の半径位置で一定の値\(C_2\ \)をもつ。 \[ \begin{align*} {\rm I.C.} \quad t=0,r \ge 0 : \quad C=C_2 \tag{7-2} \end{align*} \] 境界条件は、次のように円柱表面での蒸発を考慮し、また円柱中心の対称性から、 \[ \begin{align*} \left. \eqalign{ {\rm B.C.1} \quad r=a : \quad & -D\frac{\partial C}{\partial r}=\alpha(C_s-C_0) \cr {\rm B.C.2} \quad r=0 : \quad & \frac{\partial C}{\partial r}=0 \cr } \right\} \tag{7-3} \end{align*} \] である。ここで、\(\alpha\ \)は表面蒸発の物質移動係数、\(C_s,C_0\ \)はそれぞれ表面濃度と平衡濃度である。また\(a\ \)は円柱の半径である。
　これら基礎式を初期条件、境界条件にしたがって(7-1)式を解析的に解け。
　さらに、表面からの蒸発量の積算値の時間変化\(M_t\ \)と、無限時間放置したときの蒸発量の積算値\(M_\infty\ \)との比\(M_t/M_\infty\ \)を表す式を導出せよ。

　例題７の解答例は、こちらラプラス変換・例題７に詳述しています。

先頭に戻る

Literature Cited

　本ページを作成するにあたって、以下の書籍を参考にした。

先頭に戻る