︠273d6969-986d-4c08-8ed4-e605489fcf8di︠ %html desolve_in_japanese system:sage

Hiroshi TAKEOTO (take@pwv.co.jp)

エンジニアにとって身近な微分方程式をsageを使って解いてみます。

RC回路

RC直列回路の放電

以下のような抵抗（R）とコンデンサー（C）の回路で、

とし、$t=0$で$V_C$は$V_0=1$で、すぐに0にした場合のVc電圧の変化$V(t)$を考えてみます。

キルヒホッフの法則から回路を一周したときの電圧降下は０となります。 $$ \begin{array}{l l l} V_C + V_R = 0, V_C = V(t) & \cdots & (0) \end{array} $$ となり、コンデンサーからの電流が$i(t)=C\frac{dV(t)}{dt}$、$V_R(t)=Ri(t)$であることから、

$$ \begin{array}{l l l} V(t) + RC\frac{dV(t)}{dt} = 0 & \cdots & (1) \end{array} $$ となります。

これをsageを使って解くと以下のようになります。微分方程式の解法にはdesolve関数を使います。 

desolve(微分方程式, [求める関数と変数のリスト], [初期値のリスト])
︡7f532d18-dab5-4b26-bf6e-58a4510c97fb︡{"html": "desolve_in_japanese\r\nsystem:sage\r\n\r\n\r\n

Hiroshi TAKEOTO (take@pwv.co.jp)

\r\n

\u30a8\u30f3\u30b8\u30cb\u30a2\u306b\u3068\u3063\u3066\u8eab\u8fd1\u306a\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092sage\u3092\u4f7f\u3063\u3066\u89e3\u3044\u3066\u307f\u307e\u3059\u3002\r\n

\r\n\r\n

RC\u56de\u8def

\r\n

RC\u76f4\u5217\u56de\u8def\u306e\u653e\u96fb

\r\n\u4ee5\u4e0b\u306e\u3088\u3046\u306a\u62b5\u6297\uff08R\uff09\u3068\u30b3\u30f3\u30c7\u30f3\u30b5\u30fc\uff08C\uff09\u306e\u56de\u8def\u3067\u3001\r\n

\r\n\r\n\r\n

\r\n\u3068\u3057\u3001$t=0$\u3067$V_C$\u306f$V_0=1$\u3067\u3001\u3059\u3050\u306b0\u306b\u3057\u305f\u5834\u5408\u306eVc\u96fb\u5727\u306e\u5909\u5316$V(t)$\u3092\u8003\u3048\u3066\u307f\u307e\u3059\u3002\r\n\r\n

\r\n\u30ad\u30eb\u30d2\u30db\u30c3\u30d5\u306e\u6cd5\u5247\u304b\u3089\u56de\u8def\u3092\u4e00\u5468\u3057\u305f\u3068\u304d\u306e\u96fb\u5727\u964d\u4e0b\u306f\uff10\u3068\u306a\u308a\u307e\u3059\u3002\r\n$$\r\n\\begin{array}{l l l}\r\nV_C + V_R = 0, V_C = V(t) & \\cdots & (0)\r\n\\end{array}\r\n$$\r\n\u3068\u306a\u308a\u3001\u30b3\u30f3\u30c7\u30f3\u30b5\u30fc\u304b\u3089\u306e\u96fb\u6d41\u304c$i(t)=C\\frac{dV(t)}{dt}$\u3001$V_R(t)=Ri(t)$\u3067\u3042\u308b\u3053\u3068\u304b\u3089\u3001\r\n

\r\n$$\r\n\\begin{array}{l l l}\r\nV(t) + RC\\frac{dV(t)}{dt} = 0 & \\cdots & (1)\r\n\\end{array}\r\n$$\r\n\u3068\u306a\u308a\u307e\u3059\u3002\r\n\r\n

\r\n\u3053\u308c\u3092sage\u3092\u4f7f\u3063\u3066\u89e3\u304f\u3068\u4ee5\u4e0b\u306e\u3088\u3046\u306b\u306a\u308a\u307e\u3059\u3002\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u6cd5\u306b\u306fdesolve\u95a2\u6570\u3092\u4f7f\u3044\u307e\u3059\u3002\r\n

\r\n\r\n
\r\ndesolve(\u5fae\u5206\u65b9\u7a0b\u5f0f, [\u6c42\u3081\u308b\u95a2\u6570\u3068\u5909\u6570\u306e\u30ea\u30b9\u30c8], [\u521d\u671f\u5024\u306e\u30ea\u30b9\u30c8])\r\n
\r\n"}︡ ︠335130ab-76c9-4c90-85ce-afb86b776360︠ # 使用する変数t, C, Rと関数Vを定義します var('t C R') V = function('V',t) # 微分方程式を定義 de = V + C*R*diff(V,t) == 0 # t=0のV(t)を1として、微分方程式を解く des = desolve(de,[V,t],[0,1]);view(des) ︡62f2546a-8e80-4b98-bcf3-974b16d84479︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}e^{-\\frac{t}{C R}}"}︡ ︠3931ad38-5375-4ba5-82d7-b1cdaef8e589︠ # 解をR=1, C=1のVの変化をプロットする f(t,R,C) = des plot(f(t,1,1),[t,0,5]) ︡85f234f8-ed23-46b3-8e7d-327dd70d76fd︡{"html": ""}︡ ︠13fc6ac9-e604-4f47-9ff9-6bbfe8d9734fi︠ %html

RC直列回路の充電

つぎに、t=0から電圧V0を掛けて充電した場合のVcの変化を見てみましょう。 以下のような回路で、スイッチを1から2に入れたときの変化になります。

$$ \begin{array}{l l l} V(t) + RC\frac{dV(t)}{dt} = V_0 u(t) & \cdots & (2) \end{array} $$ ここで、$u(t)$は、単位ステップ関数、 $$ u(t)=\left\{ \begin{array}{l l} 0, & 0 \lt 0 \\ 1, & 1 \ge 0 \\ \end{array} \right. $$ です。

式(1)では右辺が０であり、このような方程式は斉次方程式と呼ばれ、式(2)のように右辺が０でない方程式は非斉次方程式と呼ばれます。 非斉次方程式の解は、 

非斉次方程式の一般解　＝　斉次方程式の一般解　＋　非斉次方程式の一つの解（特解）
から求めることができます。

$t=0$で$V(t)=0$であるとすると、$A=-V_0$となることが分かります。 従って、微分方程式(2)の求める解は、 $$ \begin{array}{l l l} V(t) = V_0(1-e^{\frac{-t}{RC}}) & \cdots & (3) \end{array} $$ となります。 ︡76c2ea5f-1936-4ebf-8ba7-19491a65fa9c︡{"html": "\r\n

RC\u76f4\u5217\u56de\u8def\u306e\u5145\u96fb

\r\n\u3064\u304e\u306b\u3001t=0\u304b\u3089\u96fb\u5727V0\u3092\u639b\u3051\u3066\u5145\u96fb\u3057\u305f\u5834\u5408\u306eVc\u306e\u5909\u5316\u3092\u898b\u3066\u307f\u307e\u3057\u3087\u3046\u3002\r\n\u4ee5\u4e0b\u306e\u3088\u3046\u306a\u56de\u8def\u3067\u3001\u30b9\u30a4\u30c3\u30c1\u30921\u304b\u30892\u306b\u5165\u308c\u305f\u3068\u304d\u306e\u5909\u5316\u306b\u306a\u308a\u307e\u3059\u3002\r\n

\r\n\r\n

\r\n$$\r\n\\begin{array}{l l l}\r\nV(t) + RC\\frac{dV(t)}{dt} = V_0 u(t) & \\cdots & (2)\r\n\\end{array}\r\n$$\r\n\u3053\u3053\u3067\u3001$u(t)$\u306f\u3001\u5358\u4f4d\u30b9\u30c6\u30c3\u30d7\u95a2\u6570\u3001\r\n$$\r\nu(t)=\\left\\{\r\n \\begin{array}{l l}\r\n 0, & 0 \\lt 0 \\\\\r\n 1, & 1 \\ge 0 \\\\\r\n \\end{array}\r\n\\right.\r\n$$\r\n\u3067\u3059\u3002\r\n\r\n

\r\n\u5f0f(1)\u3067\u306f\u53f3\u8fba\u304c\uff10\u3067\u3042\u308a\u3001\u3053\u306e\u3088\u3046\u306a\u65b9\u7a0b\u5f0f\u306f\u6589\u6b21\u65b9\u7a0b\u5f0f\u3068\u547c\u3070\u308c\u3001\u5f0f(2)\u306e\u3088\u3046\u306b\u53f3\u8fba\u304c\uff10\u3067\u306a\u3044\u65b9\u7a0b\u5f0f\u306f\u975e\u6589\u6b21\u65b9\u7a0b\u5f0f\u3068\u547c\u3070\u308c\u307e\u3059\u3002\r\n\r\n\u975e\u6589\u6b21\u65b9\u7a0b\u5f0f\u306e\u89e3\u306f\u3001\r\n

\r\n\u975e\u6589\u6b21\u65b9\u7a0b\u5f0f\u306e\u4e00\u822c\u89e3\u3000\uff1d\u3000\u6589\u6b21\u65b9\u7a0b\u5f0f\u306e\u4e00\u822c\u89e3\u3000\uff0b\u3000\u975e\u6589\u6b21\u65b9\u7a0b\u5f0f\u306e\u4e00\u3064\u306e\u89e3\uff08\u7279\u89e3\uff09\r\n
\r\n\u304b\u3089\u6c42\u3081\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u3002\r\n

\r\n\r\n\r\n\r\n$t=0$\u3067$V(t)=0$\u3067\u3042\u308b\u3068\u3059\u308b\u3068\u3001$A=-V_0$\u3068\u306a\u308b\u3053\u3068\u304c\u5206\u304b\u308a\u307e\u3059\u3002\r\n\u5f93\u3063\u3066\u3001\u5fae\u5206\u65b9\u7a0b\u5f0f(2)\u306e\u6c42\u3081\u308b\u89e3\u306f\u3001\r\n\r\n$$\r\n\\begin{array}{l l l}\r\nV(t) = V_0(1-e^{\\frac{-t}{RC}}) & \\cdots & (3)\r\n\\end{array}\r\n$$\r\n\r\n\u3068\u306a\u308a\u307e\u3059\u3002\r\n\r\n"}︡ ︠8765b41e-7687-4298-9c95-45327b3a790f︠ v(t, R, C) = (1 - des) plot(v(t, 1, 1),[t,0,5]) ︡e698d0ce-e5ac-43cf-b5ba-e0cc88e221f9︡{"html": ""}︡ ︠f162d47b-f91e-4395-836a-6eaddccb7cfdi︠ %html

ラブラス変換を使った微分方程式の解法

ラプラス変換を使用して微分方程式の解を求めることができます。ただし、求めることができるのは一般解ではなく「初期値問題」における特殊解です。

ラブラス変換とは

ラプラス変換とは、関数$f(t)$に$e^{-st}$ を掛け、t について0 から無限大まで積分したものです。その積分はs の関数$F(s)$になるが、これをもとの関数$f(t)$のラプラス変換とよび、$\mathcal{L}(f(t))$と表します。 $$ F(s) = \mathcal{L}(f(t)) = \int_{0}^{\infty}e^{-st}f(t)dt $$ ここで、$s$は$t$に無関係な変数であり、$t$は実変数である。逆に$f(t)$は$F(s)$の逆変換とよび、$\mathcal{L}^{-1}(F(s))$と表します。 ︡4043a69b-4775-45ae-9df8-cec5879bb0b0︡{"html": "\r\n

\u30e9\u30d6\u30e9\u30b9\u5909\u63db\u3092\u4f7f\u3063\u305f\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u6cd5

\r\n\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u4f7f\u7528\u3057\u3066\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u3092\u6c42\u3081\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u3002\u305f\u3060\u3057\u3001\u6c42\u3081\u308b\u3053\u3068\u304c\u3067\u304d\u308b\u306e\u306f\u4e00\u822c\u89e3\u3067\u306f\u306a\u304f\u300c\u521d\u671f\u5024\u554f\u984c\u300d\u306b\u304a\u3051\u308b\u7279\u6b8a\u89e3\u3067\u3059\u3002\r\n\r\n

\u30e9\u30d6\u30e9\u30b9\u5909\u63db\u3068\u306f

\r\n\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3068\u306f\u3001\u95a2\u6570$f(t)$\u306b$e^{-st}$ \u3092\u639b\u3051\u3001t \u306b\u3064\u3044\u30660 \u304b\u3089\u7121\u9650\u5927\u307e\u3067\u7a4d\u5206\u3057\u305f\u3082\u306e\u3067\u3059\u3002\u305d\u306e\u7a4d\u5206\u306fs \u306e\u95a2\u6570$F(s)$\u306b\u306a\u308b\u304c\u3001\u3053\u308c\u3092\u3082\u3068\u306e\u95a2\u6570$f(t)$\u306e\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3068\u3088\u3073\u3001$\\mathcal{L}(f(t))$\u3068\u8868\u3057\u307e\u3059\u3002\r\n$$\r\nF(s) = \\mathcal{L}(f(t)) = \\int_{0}^{\\infty}e^{-st}f(t)dt\r\n$$\r\n\r\n\u3053\u3053\u3067\u3001$s$\u306f$t$\u306b\u7121\u95a2\u4fc2\u306a\u5909\u6570\u3067\u3042\u308a\u3001$t$\u306f\u5b9f\u5909\u6570\u3067\u3042\u308b\u3002\u9006\u306b$f(t)$\u306f$F(s)$\u306e\u9006\u5909\u63db\u3068\u3088\u3073\u3001$\\mathcal{L}^{-1}(F(s))$\u3068\u8868\u3057\u307e\u3059\u3002\r\n"}︡ ︠06803588-9e36-4e44-bcdc-368fd49d1853︠ var('t s') f = function('f', t) # ラプラス変換のsageでの表現 view(laplace(f, t, s)) ︡b3944c04-9aa1-468e-9ad1-806f979343c9︡{"html": "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathcal{L}\\left(f\\left(t\\right), t, s\\right)"}︡ ︠844e1482-2346-4c5b-80af-3266300ca17ai︠ %html

一般的なラプラス変換の例

主な関数のラプラス変換をした結果を以下に示します。
元の関数ラプラス変換結果
$\delta$1
$1$$\frac{1}{s}$
$t$$\frac{1}{s^{2}}$
$t^2$$\frac{1}{s^{3}}$
$t^{n}$$s^{{(-n - 1)}} \Gamma\left(n + 1\right)$
$\cos\left(\omega t\right)$$\frac{s}{{(\omega^{2}+s^{2})}}$
$\sin\left(\omega t\right)$$\frac{\omega}{{(\omega^{2} + s^{2})}}$
$e^{a t}$$-\frac{1}{{(a - s)}}$
$t e^{a t}$$\frac{1}{{(a - s)}^{2}}$
︡14d18dbd-c109-4928-9838-0a1bbf78cdf6︡{"html": "\r\n

\u4e00\u822c\u7684\u306a\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u306e\u4f8b

\r\n\u4e3b\u306a\u95a2\u6570\u306e\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u3057\u305f\u7d50\u679c\u3092\u4ee5\u4e0b\u306b\u793a\u3057\u307e\u3059\u3002\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n
\u5143\u306e\u95a2\u6570\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u7d50\u679c
$\\delta$1
$1$$\\frac{1}{s}$
$t$$\\frac{1}{s^{2}}$
$t^2$$\\frac{1}{s^{3}}$
$t^{n}$$s^{{(-n - 1)}} \\Gamma\\left(n + 1\\right)$
$\\cos\\left(\\omega t\\right)$$\\frac{s}{{(\\omega^{2}+s^{2})}}$
$\\sin\\left(\\omega t\\right)$$\\frac{\\omega}{{(\\omega^{2} + s^{2})}}$
$e^{a t}$$-\\frac{1}{{(a - s)}}$
$t e^{a t}$$\\frac{1}{{(a - s)}^{2}}$
\r\n\r\n"}︡ ︠2b6fdca7-1650-41ad-bd19-cead4c2c84fc︠ # 変数の宣言と制約条件n > 0をセット var('n omega a') assume(n>0) # 主な関数のラプラス変換の結果 print latex(dirac_delta), laplace(dirac_delta(t), t, s) print 1, laplace(1, t, s) print t, laplace(t, t, s) print t^2, laplace(t^2, t, s) print latex(t^n), laplace(t^n, t, s) print latex(cos(omega*t)), laplace(cos(omega*t), t, s) print latex(sin(omega*t)), laplace(sin(omega*t), t, s) print latex(exp(a*t)), laplace(exp(a*t), t, s) print latex(t*exp(a*t)), laplace(t*exp(a*t), t, s) ︡15132aa3-4178-4dc4-a93f-8dabf334945e︡{"stdout": "\\delta 1\n1 1/s\nt s^(-2)\nt^2 2/s^3\nt^{n} s^(-n - 1)*gamma(n + 1)\n\\cos\\left(\\omega t\\right) s/(omega^2 + s^2)\n\\sin\\left(\\omega t\\right) omega/(omega^2 + s^2)\ne^{a t} -1/(a - s)\nt e^{a t} (a - s)^(-2)"}︡ ︠9d0f43a2-3c71-45e2-90e9-428e4d35d845i︠ %html

ラプラス変換の性質

ラプラス変換の性質で、特質すべきはラプラス変換の微分によって$s$が掛けられる点です。 $$ \mathcal{L}(f') = s\mathcal{L}(f) + f(0) $$ これをsageを使って表現すると以下のようになります。 ︡7533ba21-21a8-4ca9-8c7e-fefaacb30a87︡{"html": "\r\n

\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u306e\u6027\u8cea

\r\n\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u306e\u6027\u8cea\u3067\u3001\u7279\u8cea\u3059\u3079\u304d\u306f\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u306e\u5fae\u5206\u306b\u3088\u3063\u3066$s$\u304c\u639b\u3051\u3089\u308c\u308b\u70b9\u3067\u3059\u3002\r\n$$\r\n\\mathcal{L}(f') = s\\mathcal{L}(f) + f(0)\r\n$$\r\n\r\n\u3053\u308c\u3092sage\u3092\u4f7f\u3063\u3066\u8868\u73fe\u3059\u308b\u3068\u4ee5\u4e0b\u306e\u3088\u3046\u306b\u306a\u308a\u307e\u3059\u3002\r\n"}︡ ︠15b34879-a4af-456a-bff2-037f20e47831︠ f = function('f', t) laplace(diff(f,t), t, s) ︡4946c743-8a2e-4cfc-ad43-31ca5fa01804︡{"stdout": "s*laplace(f(t), t, s) - f(0)"}︡ ︠61c68922-ce48-4b78-9f1e-1a4c0519dde7i︠ %html

微分方程式の解法

ラブラス変換を使った微分方程式の解法は、以下の手順で行います。
  1. 微分方程式にラプラス変換を施し、微分方程式を変数t領域から変数s領域へ移します
  2. 得られた方程式は演算子$s$の代数方程式となるので、代数計算により解を求める
  3. 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換する
それでは、sageで微分方程式(1)をラプラス変換を使って解くと、以下のようになります。 ︡8a0f6583-b332-4456-8534-35f60da0a2eb︡{"html": "\r\n

\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u6cd5

\r\n\u30e9\u30d6\u30e9\u30b9\u5909\u63db\u3092\u4f7f\u3063\u305f\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u6cd5\u306f\u3001\u4ee5\u4e0b\u306e\u624b\u9806\u3067\u884c\u3044\u307e\u3059\u3002\r\n
    \r\n
  1. \u5fae\u5206\u65b9\u7a0b\u5f0f\u306b\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u65bd\u3057\u3001\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092\u5909\u6570t\u9818\u57df\u304b\u3089\u5909\u6570s\u9818\u57df\u3078\u79fb\u3057\u307e\u3059
    2. \r\n
  2. \u5f97\u3089\u308c\u305f\u65b9\u7a0b\u5f0f\u306f\u6f14\u7b97\u5b50$s$\u306e\u4ee3\u6570\u65b9\u7a0b\u5f0f\u3068\u306a\u308b\u306e\u3067\u3001\u4ee3\u6570\u8a08\u7b97\u306b\u3088\u308a\u89e3\u3092\u6c42\u3081\u308b
    3. \r\n
  3. \u6c42\u3081\u305fs \u95a2\u6570\u306e\u89e3\u3092\u9006\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u7528\u3044\u3066t \u306e\u95a2\u6570\u306b\u5909\u63db\u3059\u308b
    4. \r\n
\r\n\r\n\u305d\u308c\u3067\u306f\u3001sage\u3067\u5fae\u5206\u65b9\u7a0b\u5f0f(1)\u3092\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u4f7f\u3063\u3066\u89e3\u304f\u3068\u3001\u4ee5\u4e0b\u306e\u3088\u3046\u306b\u306a\u308a\u307e\u3059\u3002\r\n"}︡ ︠f0870190-3ee8-43d4-a858-6ddb23364c3f︠ # deをラプラス変換する l1 = laplace(de, t, s); l1 ︡7e52b92e-99c9-4e8b-a528-c3e254d50a7a︡{"stdout": "(s*laplace(V(t), t, s) - V(0))*C*R + laplace(V(t), t, s) == 0"}︡ ︠516269fd-411f-4ffc-b8cc-814bb61169d8︠ # この方程式をlapace(V(t), t, s)について解くと s1 = solve(l1, laplace(V(t), t, s)); show(s1) # 解の右辺に初期値V(0)=1を代入する s11 = s1[0].rhs().subs_expr(V(0) == 1); s11 ︡32e9475c-29fe-430e-b6d3-f3041c342cb7︡{"html": "
\\left[\\mathcal{L}\\left(V\\left(t\\right), t, s\\right) = \\frac{C R V\\left(0\\right)}{{(C R s + 1)}}\\right]
"}︡{"stdout": "C*R/(C*R*s + 1)"}︡ ︠4236e8d1-e229-45e4-b2b2-2e960e051b73︠ # 得られた解C*R/(C*R*s + 1)を逆ラプラス変換する inverse_laplace(s11, s, t) ︡00a3c993-60e2-41f9-a2e8-a42ee7e0eba3︡{"stdout": "e^(-t/(C*R))"}︡ ︠c2be3e64-e7bf-41e7-9856-40011fff3a0ai︠ %html これで、ラプラス変換からも$e^{\frac{-t}{RC}}$の解を得ることができました。

次に微分方程式(2)に対しても同様にやってみます。 ︡a1f9f592-d0cb-435a-be17-fe8917110eb0︡{"html": "\r\n\u3053\u308c\u3067\u3001\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u304b\u3089\u3082$e^{\\frac{-t}{RC}}$\u306e\u89e3\u3092\u5f97\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3057\u305f\u3002

\r\n\u6b21\u306b\u5fae\u5206\u65b9\u7a0b\u5f0f(2)\u306b\u5bfe\u3057\u3066\u3082\u540c\u69d8\u306b\u3084\u3063\u3066\u307f\u307e\u3059\u3002\r\n"}︡ ︠0241e9fb-3c1f-488e-b148-f0278f87dcb4︠ # 微分方程式(2)を定義する de2 = V + C*R*diff(V,t) == unit_step(t) l2 = laplace(de2, t, s); l2 # ラプラス変換した後、 s2 = solve(l2, laplace(V(t), t, s)); show(s2) ︡5743bd8f-9a49-4727-8ac9-df67aa3307fc︡{"html": "

\\left[\\mathcal{L}\\left(V\\left(t\\right), t, s\\right) = \\frac{{(C R s V\\left(0\\right) + 1)}}{{(C R s^{2} + s)}}\\right]
"}︡ ︠46f75cac-a178-426e-8229-22e7c6a7fc7c︠ # 得られた解に初期値V(0)=0を代入し s22 = s2[0].rhs().subs_expr(V(0) == 0); show(s22) # 逆ラプラス変換する inverse_laplace(s22, s, t) ︡153b89ea-9009-4895-9f51-e3039300e439︡{"html": "
\\frac{1}{C R s^{2} + s}
"}︡{"stdout": "-e^(-t/(C*R)) + 1"}︡ ︠aa040267-a819-4277-add6-f7e3cca884bei︠ %html 微分方程式(1)と同様微分方程式(2)も期待した解を得ることができました。

maximaを使った微分方程式の解法

sageでもある程度の微分方程式は解くことができますが、初期値の柔軟性という点では、maximaに一日の長があるようです。

次に以下の微分方程式をmaximaを使って解いてみます。 $$ x^2\frac{dy}{dx}+3xy = \frac{sin(x)}{x}, y(\pi) = 0 $$ maximaを使った微分方程式の解法は以下の手順で行います。

_sage_()関数は、別の処理系（maxima, scilab, R等）の結果をsageの表現に変換するために使用します。 ︡e248f1ca-873d-4dd8-b153-37ac8eecf178︡{"html": "\r\n\u5fae\u5206\u65b9\u7a0b\u5f0f(1)\u3068\u540c\u69d8\u5fae\u5206\u65b9\u7a0b\u5f0f(2)\u3082\u671f\u5f85\u3057\u305f\u89e3\u3092\u5f97\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3057\u305f\u3002\r\n\r\n\r\n\r\n

maxima\u3092\u4f7f\u3063\u305f\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u6cd5

\r\nsage\u3067\u3082\u3042\u308b\u7a0b\u5ea6\u306e\u5fae\u5206\u65b9\u7a0b\u5f0f\u306f\u89e3\u304f\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u304c\u3001\u521d\u671f\u5024\u306e\u67d4\u8edf\u6027\u3068\u3044\u3046\u70b9\u3067\u306f\u3001maxima\u306b\u4e00\u65e5\u306e\u9577\u304c\u3042\u308b\u3088\u3046\u3067\u3059\u3002

\r\n\r\n\u6b21\u306b\u4ee5\u4e0b\u306e\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092maxima\u3092\u4f7f\u3063\u3066\u89e3\u3044\u3066\u307f\u307e\u3059\u3002\r\n$$\r\nx^2\\frac{dy}{dx}+3xy = \\frac{sin(x)}{x}, y(\\pi) = 0\r\n$$\r\n\r\nmaxima\u3092\u4f7f\u3063\u305f\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u6cd5\u306f\u4ee5\u4e0b\u306e\u624b\u9806\u3067\u884c\u3044\u307e\u3059\u3002\r\n

\r\n\r\n_sage_()\u95a2\u6570\u306f\u3001\u5225\u306e\u51e6\u7406\u7cfb\uff08maxima, scilab, R\u7b49\uff09\u306e\u7d50\u679c\u3092sage\u306e\u8868\u73fe\u306b\u5909\u63db\u3059\u308b\u305f\u3081\u306b\u4f7f\u7528\u3057\u307e\u3059\u3002\r\n"}︡ ︠9b2bba5f-b82b-404d-90ca-94a93903ea3e︠ # 変数x, yを定義 var('x y') # maximaを使って微分方程式を定義：diffの前の'に注意 deq = maxima("x^2*'diff(y,x) + 3*x*y = sin(x)/x") ︡10e07fff-6a23-401a-b883-12d1f9c9cdb3︡︡ ︠8f82e2dc-30eb-4619-8f31-e89baea1aaab︠ # 微分方程式を解く d2 = deq.ode2(y,x).ic1(x=pi,y=0) # 解を表示 show(deq.ode2(y,x)) # maxima の結果をsageの表現に変換するには_sage()_を使い、rhs()を使ってyの関数のみを取り出すことができる d3 = d2._sage_(); d3.rhs() ︡0a563c46-adde-4382-83fb-cbe5d5c47598︡{"html": "
y={{{\\it c}-\\cos x}\\over{x^3}}
"}︡{"stdout": "-(cos(x) + 1)/x^3"}︡ ︠22589209-fb2e-47a5-b4bb-c0dcde97eb4b︠ # ic1を使って初期値を指定して微分方程式を解く d4 = deq.ode2(y,x).ic1(x=pi,y=0); d4 ︡3591cf90-5797-486c-8ca0-0f0638dce6f1︡{"stdout": "y=-(cos(x)+1)/x^3"}︡ ︠8c0b52a8-abd1-47c3-9975-8445f7f2db29︠ # 結果をプロットする plot(d4._sage_().rhs(), [x, 0, pi]) ︡53a24dd4-1888-427c-90d9-20b67fa53ff8︡{"html": ""}︡ ︠4455ecde-63ed-4549-a471-c35613f566f4i︠ %html このようにsageを使って微分方程式を様々な手法で解くことにより、解の検証も可能です。

是非、この機会にsageを使って微分方程式を解いてみてください。 ︡fa25e81c-0c45-4980-aa73-cf11293aad38︡{"html": "\r\n\u3053\u306e\u3088\u3046\u306bsage\u3092\u4f7f\u3063\u3066\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092\u69d8\u3005\u306a\u624b\u6cd5\u3067\u89e3\u304f\u3053\u3068\u306b\u3088\u308a\u3001\u89e3\u306e\u691c\u8a3c\u3082\u53ef\u80fd\u3067\u3059\u3002

\r\n\r\n\u662f\u975e\u3001\u3053\u306e\u6a5f\u4f1a\u306bsage\u3092\u4f7f\u3063\u3066\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092\u89e3\u3044\u3066\u307f\u3066\u304f\u3060\u3055\u3044\u3002\r\n"}︡