Return to SymPy で理工系の数学C

SymPy で常微分方程式

必要なパッケージの import と設定

In [1]:
# 1文字変数の Symbol の定義が省略できる
from sympy.abc import *
# a Python library for symbolic mathematics
from sympy import *

# SymPy Plotting Backends (SPB): グラフを描く際に利用
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

変数分離法

例題

$$ \frac{dy}{dx} = -2 x\, y $$

  • y を（変数ではなく）「関数」として定義するには y = Function('y')
  • 方程式・等式は Eq(左辺, 右辺)
  • 微分であることを表すには Derivative()（diff(y(x), x) だと実際に微分してしまう）
In [2]:
# y を x の関数として宣言
y = Function('y')

# 解くべき微分方程式
eq = Eq(Derivative(y(x), x), -2*x*y(x))
eq
Out[2]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = -2 x y{\left(x \right)}$

SymPy で微分方程式を解くには，以下のように dsolve() 関数を使います。

In [3]:
dsolve(eq, y(x))
Out[3]:
$\displaystyle y{\left(x \right)} = C_{1} e^{-x^{2}}$

上の出力で $C_1$ は積分定数です。

ちなみに，以下のように，解くべき方程式を右辺がゼロになるように

$$\frac{dy}{dx} + 2 x\, y =0$$

として，左辺のみを dsolve() の引数にしても，同じ答えがでます。

In [4]:
dsolve(Derivative(y(x), x)+2*x*y(x), y(x))
Out[4]:
$\displaystyle y{\left(x \right)} = C_{1} e^{-x^{2}}$

マルサスの人口モデル

$$\frac{dN}{dt} = \gamma \, N$$

を初期条件 $t = t_0$ で $N(t_0) = N_0$ として解く。

In [5]:
# 関数および変数の宣言
N = Function('N')
var('gamma')

eqm = Eq(Derivative(N(t), t), gamma*N(t))
eqm
Out[5]:
$\displaystyle \frac{d}{d t} N{\left(t \right)} = \gamma N{\left(t \right)}$
In [6]:
# 初期値変数の宣言
var('t0, N0')

# 初期条件は ics={} で以下のように...
dsolve(eqm, N(t), ics={N(t).subs(t,t0):N0})
Out[6]:
$\displaystyle N{\left(t \right)} = N_{0} e^{\gamma t} e^{-\gamma t_{0}}$

参考：米国の人口とマルサスモデル

上記によれば，1790年から1930年のアメリカの人口は以下のようになっている（人口の単位は百万人）。

In [7]:
var('usa')
# 行列 Matrix を使ってみる。
usa = Matrix([
# 西暦, 人口（百万人）
  [1790,  3.9], 
  [1800,  5.3],
  [1810,  7.2],
  [1820,  9.6],
  [1830, 12.9],
  [1840, 17.1],
  [1850, 23.2],
  [1860, 31.4],
  [1870, 38.6],
  [1880, 50.2],
  [1890, 62.9],
  [1900, 76.0],
  [1910, 92.0],
  [1920,106.5],
  [1930,123.2]
])

$t_0 = 1790$ （年）とすると，$N_0 = 3.9$ 。

Python の要素を取り出すインデックスは 0 ゼロはじまりであるから…

In [8]:
t0 = usa[0, 0] # 1行目の1列目
N0 = usa[0, 1] # 1行目の2列目
t0, N0
Out[8]:
(1790, 3.90000000000000)

残りのパラメーター $\gamma$ は，別の時刻 $t_1$ における $N(t_1)$ から求めてみる。

\begin{eqnarray}
N(t_1) &=& N_0 e^{\gamma\, (t_1 -t_0)} \\
\log \frac{N(t_1)}{N_0} &=& \gamma\, (t_1 -t_0) \\
\therefore\ \ \gamma &=& \frac{1}{t_1 -t_0} \log \frac{N(t_1)}{N_0}
\end{eqnarray}

たとえば，$t_1 = 1830$ とすると（Python のインデックスは 0 ゼロはじまりだから）…

In [9]:
t1 = usa[4, 0]  # 5行目の1列目 
N1 = usa[4, 1]  # 5行目の2列目 

1/(t1-t0) * log(N1/N0)
Out[9]:
$\displaystyle 0.0299062689558006$

したがって，マルサスの人口モデルを $N_m(t)$ とすると

In [10]:
def Nm(t):
    t0 = usa[0, 0]
    N0 = usa[0, 1]
    gamma = 1/(t1-t0) * log(N1/N0)
    return N0 * exp(gamma*(t-t0))

このようにして求められた $N_m(t)$ を，人口データ usa と共グラフにしてみます。

In [11]:
# SymPy Plotting Backends (SPB) で陽関数を描く

p1 = plot(Nm(t), (t, 1790, 1940), "マルサスモデル", 
          legend = True, 
          xlabel = "年", ylabel = "人口（単位：百万人）");

SPB で点，$x$ 座標 $y$ 座標の数値データをグラフにする際の書式は

plot_list([x1, ..., xn], [y1, ..., yn])

ですので，

  • 行列 usa の 1 列目 .col(0) をリストにした list(usa.col(0)) を$x$ 座標，
  • 2 列目 list(usa.col(1)) を$y$ 座標,
  • 数値データを点 is_poinst = True

でグラフにします。.col() を使いたいがために，usa を行列としたのでした。

In [12]:
p2 = plot_list(list(usa.col(0)), list(usa.col(1)), 
               "アメリカの人口", 
               # 線で繋がす点で, 凡例を表示
               is_point = True, legend = True, 
               xlabel = "年", ylabel = "人口（単位：百万人）");

SPB で上記2つのグラフを重ねて表示させる例。

In [13]:
p3 = p1 + p2
p3.show();

マルサス・モデルは $\gamma > 0$ の場合には人口の際限ない増加（指数関数的増加）を予測する。しかし，実際にはいろいろな要因により，このような無制限な増加は続かない。

上のグラフでも，1870年あたりから，マルサス・モデルのグラフは実際のアメリカの人口データからずれはじめ，1900年以降は全くあっていない。

ヴェアフルストによる修正人口モデル

$$
\frac{dN}{dt} = \gamma N \left(1 -\frac{N}{N_{\rm max}}\right)
$$

ヴェアフルストの修正人口モデルの方程式（ロジスティック方程式）を eqv とし，この微分方程式を dsolve() で解きます。

In [14]:
N = Function('N')
var('gamma Nmax')

eqv = Eq(Derivative(N(t), t), gamma*N(t)*(1-N(t)/Nmax))
eqv
Out[14]:
$\displaystyle \frac{d}{d t} N{\left(t \right)} = \gamma \left(1 -\frac{N{\left(t \right)}}{Nmax}\right) N{\left(t \right)}$
In [15]:
var('t0 N0')

ansv = dsolve(eqv, N(t), ics = {N(t).subs(t,t0):N0})
ansv.simplify()
Out[15]:
$\displaystyle N{\left(t \right)} = \frac{N_{0} Nmax e^{\gamma t}}{N_{0} e^{\gamma t} -\left(N_{0} -Nmax\right) e^{\gamma t_{0}}}$

$\displaystyle n(t) \equiv \frac{N(t)}{N_0}, \ nmax \equiv \frac{N_{max}}{N_0}$ とおくと，上の解は以下のような「規格化された」解として書けます。

$$n(t) = \frac{nmax}{1 -(1 -nmax)\, e^{\gamma\, (t_0 -t)}}$$

参考：米国の人口とヴェアフルストモデル

初期条件を $t_0 = 1790$（年）のとき $N(t_0) = N_0 = 3.9$（百万人）とします。

In [16]:
t0 = usa[0, 0]
N0 = usa[0, 1]

規格化された解 $n(t)$ では $n(t_0) = 1$ です。$n(t)$ の2つのパラメータ $nmax$ と $\gamma$ は以下のように，2つの年のデータから求めることにします。

まず，関数 $n(t)$ を定義します。

In [17]:
var('nmax')
def n(t):
    t0 = usa[0, 0]
    return nmax/(1 -(1-nmax)*exp(gamma*(t0-t)))
n(t)
Out[17]:
$\displaystyle \frac{nmax}{-\left(1 -nmax\right) e^{\gamma \left(1790 -t\right)} + 1}$

例として，t1 = 1850（年）と t2 = 1910（年）のデータを使ってみます。

In [18]:
t1 = usa[6, 0]
N1 = usa[6, 1]
t1, N1
Out[18]:
(1850, 23.2000000000000)
In [19]:
t2 = usa[12, 0]
N2 = usa[12, 1]
t2, N2
Out[19]:
(1910, 92.0000000000000)

この2組のデータを使って，未知変数 $nmax$ と $\gamma$ についての連立方程式をつくります。

In [20]:
Eq(N1/N0, n(t1))
Out[20]:
$\displaystyle 5.94871794871795 = \frac{nmax}{-\left(1 -nmax\right) e^{-60 \gamma} + 1}$
In [21]:
Eq(N2/N0, n(t2))
Out[21]:
$\displaystyle 23.5897435897436 = \frac{nmax}{-\left(1 -nmax\right) e^{-120 \gamma} + 1}$

簡単な代数方程式にするため，$e^{-60 \gamma} \equiv T$ とおきます。

In [22]:
var('T')
eqv1 = Eq(N1/N0, n(t1).subs(exp(-60*gamma), T))
eqv1
Out[22]:
$\displaystyle 5.94871794871795 = \frac{nmax}{-T \left(1 -nmax\right) + 1}$
In [23]:
eqv2 = Eq(N2/N0, n(t2).subs(exp(-60*gamma), T))
eqv2
Out[23]:
$\displaystyle 23.5897435897436 = \frac{nmax}{-T^{2} \cdot \left(1 -nmax\right) + 1}$

eqv1, eqv2 を連立させて $nmax$, $T$ について解きます。

In [24]:
ans = solve([eqv1, eqv2], [nmax, T])
ans
Out[24]:
[(49.9687003440904, 0.151115116017121)]

一般には，連立方程式から複数組の解が得られる可能性がありますが，今回は1組のみです。

ans[0] として1組目の解を表示します。

In [25]:
ans[0]
Out[25]:
(49.9687003440904, 0.151115116017121)

こうして得られた $nmax$, $T$ の値を元の $N_{max}$, $\gamma$ の値になおし，
ヴェアフルストモデルの解 $N(t)$ に入れた関数を Nv(t) として定義します。

In [26]:
print('Nmax = nmax*N0 =', ans[0][0]*N0)
print('gamma = -log(T)/60 =', -log(ans[0][1])/60)

Nmax = nmax*N0 = 194.877931341953
gamma = -log(T)/60 = 0.0314952229147497
In [27]:
nmax0, T0 = ans[0]

def Nv(t):
    t0 = usa[0, 0]
    N0 = usa[0, 1]
    Nmax = nmax0 * N0
    gamma = -1/60 * log(T0)
    return N0*Nmax/(N0-(N0-Nmax)*exp(gamma*(t0-t)))

グラフを描きます。

In [28]:
p4 = plot(Nv(t), (t, 1790, 1940), "ヴェアフルストモデル", 
          legend = True, 
          xlabel = "年", ylabel = "人口（単位：百万人）");

3つのグラフを一緒に表示します。

In [29]:
p5 = p3+p4
p5.show();

このグラフをみると，（$N_{max}$ と $\gamma$ をうまくあわせることによって）ヴェアフルストモデルのほうが，この年代の米国の人口をよくあらわしていることがわかる。

1階線形微分方程式と積分因子法

SymPy では，1階線形微分方程式は，特に積分因子法を coding しなくても，dsolve() 関数で解くことができます。

例題

$$ \frac{dy}{dx} + \frac{y}{x} = \frac{\sin x}{x} $$

SymPy での三角関数 $\sin x, \cos x, \tan x$ は sin(x), cos(x), tan(x) と書きます。

In [30]:
y = Function('y')

eq = Eq(Derivative(y(x), x)+y(x)/x, sin(x)/x)
eq
Out[30]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} + \frac{y{\left(x \right)}}{x} = \frac{\sin{\left(x \right)}}{x}$
In [31]:
dsolve(eq, y(x))
Out[31]:
$\displaystyle y{\left(x \right)} = \frac{C_{1} -\cos{\left(x \right)}}{x}$

参考：あえて積分因子法の公式で…

この例題をあえて積分因子法の公式に従ってやってみる。

一般的な1階線形微分方程式（非同次方程式の格好しているもの）を

$$\frac{dy}{dx} + P(x) y = Q(x)$$

と書くと，積分因子 $g(x)$ は

$$g(x) = \exp \left\{\int^x P(x’)\,dx’ \right\}$$

解は

$$y = \frac{1}{g(x)} \left\{\int^x g(x’) Q(x’) \, dx’ + C \right\}$$

この例題では

$$P(x) = \frac{1}{x}, \ Q(x) = \frac{\sin x}{x}$$

In [32]:
def P(x):
    return 1/x

def Q(x):
    return sin(x)/x

g = exp(integrate(P(x), x))
g
Out[32]:
$\displaystyle x$
In [33]:
y = 1/g * (integrate(g * Q(x), x) + C)
y
Out[33]:
$\displaystyle \frac{C -\cos{\left(x \right)}}{x}$

最も簡単な定数係数2階微分方程式

最も簡単な定数係数2階線形微分方程式 $y” + y = 0$ と $y” -y = 0$ を解く。

In [34]:
y = Function('y')

eq1 = Derivative(y(x), x, 2) + y(x)
eq1
Out[34]:
$\displaystyle y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)}$
In [35]:
dsolve(eq1, y(x))
Out[35]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(x \right)} + C_{2} \cos{\left(x \right)}$
In [36]:
eq2 = Derivative(y(x), x, 2) - y(x)
eq2
Out[36]:
$\displaystyle -y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)}$
In [37]:
dsolve(eq2, y(x))
Out[37]:
$\displaystyle y{\left(x \right)} = C_{1} e^{-x} + C_{2} e^{x}$

参考：脊髄反射によらずに解く

$y” + y = 0$ を脊髄反射によらずに解く。両辺に $2 y’$ をかけて整理すると

$$\left(\left(y’\right)^2 + y^2\right)’ = 0$$

微分してゼロということは，かっこの中身は定数であるということなので，

$$ \left(y’\right)^2 + y^2 = \mbox{const.} \equiv a^2$$

とおける。つまり

$$\frac{dy}{dx} = \pm \sqrt{a^2 -y^2}$$

を解けばよい。まずは $+$ の式から…

$\displaystyle \frac{dy}{dx} = + \sqrt{a^2 – y^2}$ を解きます。

注意すべきことは，$\displaystyle \sqrt{a^2}$ です。SymPy がどのような答えを返すか見てみると…

In [38]:
var('a')
sqrt(a**2)
Out[38]:
$\displaystyle \sqrt{a^{2}}$

$\displaystyle \sqrt{a^2} = |a|$ ですが SymPy はそのまま返します。では $a > 0$ と宣言すると…

In [39]:
var('a', positive = True)
sqrt(a**2)
Out[39]:
$\displaystyle a$

参考までに，$a < 0$ と宣言すれば…

In [40]:
var('a', negative = True)
sqrt(a**2)
Out[40]:
$\displaystyle – a$

以下では簡単のために，$a > 0$ として話を進めます。

In [41]:
# 一般性を失うことなく第1の積分定数の平方根 a を正にできる
var('a', positive = True)

eq1 = Eq(Derivative(y(x), x), sqrt(a**2 - y(x)**2))
eq1
Out[41]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = \sqrt{a^{2} -y^{2}{\left(x \right)}}$
In [42]:
dsolve(eq1, y(x))
Out[42]:
$\displaystyle y{\left(x \right)} = a \sin{\left(C_{1} + x \right)}$

次に $\displaystyle \frac{dy}{dx} = -\sqrt{a^2 -y^2}$ を解きます。

In [43]:
eq2 = Eq(Derivative(y(x), x), -sqrt(a**2 - y(x)**2))
eq2
Out[43]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = -\sqrt{a^{2} -y^{2}{\left(x \right)}}$
In [44]:
dsolve(eq2, y(x))
Out[44]:
$\displaystyle y{\left(x \right)} = a \sin{\left(C_{1} -x \right)}$

$y(x) = a \sin\left(C_1 -x\right) = -a \sin(x -C_1)$
ということで 2 つの解をまとめると，積分定数あらためて $\pm a \Rightarrow C_2, \pm C_1 \Rightarrow C_1$ とおいて

\begin{eqnarray}
y &=& C_2 \sin(x + C_1) \\
&=& C_2 \left(\sin x \cos C_1 + \cos x \sin C_1 \right) \\
&=& (C_2 \cos C_1) \sin x + (C_2 \sin C_1) \cos x \\
&\equiv& A \cos x + B \sin x\\
\mbox{ここで} \ \ A &\equiv& C_2 \sin C_1 \\
B &\equiv& C_2 \cos C_1
\end{eqnarray}

のように書けるということがわかる。

最も簡単な定数係数2階微分方程式：続き

$y” + K y = 0$ あるいは移項して $y” = -K y$ を解く。

$K > 0$ の場合

$K > 0$ は SymPy では var('K', positive = True)

In [45]:
y = Function('y')
var('x')
var('K', positive = True)

eq1 = Eq(Derivative(y(x), x, 2), -K*y(x))
eq1
Out[45]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = -K y{\left(x \right)}$
In [46]:
dsolve(eq1, y(x))
Out[46]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(\sqrt{K} x \right)} + C_{2} \cos{\left(\sqrt{K} x \right)}$

$K < 0$ の場合

$K < 0$ は SymPy では var('K', negative = True)

In [47]:
var('K', negative = True)

eq2 = Eq(Derivative(y(x), x, 2), -K*y(x))
eq2
Out[47]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = -K y{\left(x \right)}$
In [48]:
dsolve(eq2, y(x))
Out[48]:
$\displaystyle y{\left(x \right)} = C_{1} e^{-x \sqrt{-K}} + C_{2} e^{x \sqrt{-K}}$

$K = 0$ の場合

In [49]:
eq3 = Eq(Derivative(y(x), x, 2), 0)
eq3
Out[49]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 0$
In [50]:
dsolve(eq3, y(x))
Out[50]:
$\displaystyle y{\left(x \right)} = C_{1} + C_{2} x$

人類の至宝：オイラーの公式

SymPy がオイラーの公式を知っているか確認する。

SymPy では

  • 虚数単位 $i$ は I
  • 円周率 $\pi$ は pi
  • ネイピア数 $e$ は E

で定義されています。

指数関数 $e^{i x}$ を .rewrite(cos) をつけて三角関数で書き直してみます。

In [51]:
Eq(exp(I*x), 
   exp(I*x).rewrite(cos))
Out[51]:
$\displaystyle e^{i x} = i \sin{\left(x \right)} + \cos{\left(x \right)}$

オイラーの項式 $e^{i x} = \cos x + i \sin x$ が確認できました。

オイラーの等式

$x = \pi$ のとき，$e^{i \pi} + 1 = 0$。指数関数 exp() ではなく，ネイピア数 E を使ってみます。

In [52]:
E**(I*pi) + 1
Out[52]:
$\displaystyle 0$

左辺と右辺が等しいかどうかを確認するときには == を使います。

In [53]:
E**(I*pi) + 1 == 0
Out[53]:
True

オイラーの公式からみた三角関数と双曲線関数の関係

三角関数 $\cos x, \sin x$ を .rewrite(exp) で（複素）指数関数表示します。

In [54]:
cos(x).rewrite(exp)
Out[54]:
$\displaystyle \frac{e^{i x}}{2} + \frac{e^{-i x}}{2}$
In [55]:
sin(x).rewrite(exp)
Out[55]:
$\displaystyle -\frac{i \left(e^{i x} -e^{-i x}\right)}{2}$

双曲線関数 $\cosh x, \sinh x$ も指数関数表示します。

In [56]:
cosh(x).rewrite(exp)
Out[56]:
$\displaystyle \frac{e^{x}}{2} + \frac{e^{-x}}{2}$
In [57]:
sinh(x).rewrite(exp)
Out[57]:
$\displaystyle \frac{e^{x}}{2} -\frac{e^{-x}}{2}$

…というわけで，以下のことが確認できました。

\begin{eqnarray}
\cos x &=& \frac{e^{ix} + e^{-ix}}{2} \\
\sin x &=& \frac{e^{ix} -e^{-ix}}{2 i} \\
\cosh x &=& \frac{e^x + e^{-x}}{2} \\
\sinh x &=& \frac{e^x -e^{-x}}{2}
\end{eqnarray}

さて，三角関数や双曲線関数の変数（引数）が虚数でもいいのだと拡張すると…

In [58]:
cosh(I*x)
Out[58]:
$\displaystyle \cos{\left(x \right)}$
In [59]:
sinh(I*x)
Out[59]:
$\displaystyle i \sin{\left(x \right)}$
In [60]:
cos(I*x)
Out[60]:
$\displaystyle \cosh{\left(x \right)}$
In [61]:
sin(I*x)
Out[61]:
$\displaystyle i \sinh{\left(x \right)}$

… ということで以下のことを確認できました。

\begin{eqnarray}
\cosh (i x) &=& \cos x \\
\sinh (i x) &=& i \sin x \\
\cos (i x) &=& \cosh x \\
\sin (i x) &=& i \sinh x
\end{eqnarray}

例題

微分方程式 $y” + K y = 0$ の $K > 0$ の場合の解は

\begin{eqnarray}
y &=& A \cos\left( \sqrt{K} x\right) + \frac{B}{\sqrt{K}} \sin\left( \sqrt{K} x\right)
\end{eqnarray}

のように書ける。このことを使って，$K<0$ の場合と $K = 0$ の場合の解を上記の解から直接求める。

In [62]:
# 念のために定数の宣言
var('A B x')
var('K', positve =True)

y1 = A*cos(sqrt(K)*x) + B/sqrt(K)*sin(sqrt(K)*x)
y1
Out[62]:
$\displaystyle A \cos{\left(\sqrt{K} x \right)} + \frac{B \sin{\left(\sqrt{K} x \right)}}{\sqrt{K}}$
$K < 0$ の場合

.subs(K, -abs(K)) で $K$ に $-|K|$ を代入すると…

In [63]:
y2 = y1.subs(K, -abs(K))
y2
Out[63]:
$\displaystyle A \cosh{\left(x \sqrt{\left|{K}\right|} \right)} + \frac{B \sinh{\left(x \sqrt{\left|{K}\right|} \right)}}{\sqrt{\left|{K}\right|}}$
$K = 0$ の場合

$K \rightarrow 0$ の極限をとります。

In [64]:
Eq(Limit(y1, K, 0), 
   Limit(y1, K, 0).doit())
Out[64]:
$\displaystyle \lim_{K \to 0^+}\left(A \cos{\left(\sqrt{K} x \right)} + \frac{B \sin{\left(\sqrt{K} x \right)}}{\sqrt{K}}\right) = A + B x$

… ということで，$K>0$ の場合の解を1つおぼえれば，$K < 0$ の場合も，$ K =0$ の場合もそれからわかる，ということが確認できました。

定数係数2階線形同次方程式

定数係数2階線形微分方程式（同次方程式）は以下のように書ける。

$$ \frac{d^2 y}{dx^2} + 2 b \frac{dy}{dx} + cy = 0$$

$b, c$ は定数。

In [65]:
y = Function('y')
var('x')
var('b c')

eq = Derivative(y(x), x, 2) + 2*b*Derivative(y(x), x) + c*y(x)
eq
Out[65]:
$\displaystyle 2 b \frac{d}{d x} y{\left(x \right)} + c y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)}$

特に仮定なしで eq = 0 を dsolve() したときの解。ルートの中身 $b^2 -c >0$ を仮定したときの解になっている。

In [66]:
ans1 = dsolve(eq, y(x))
ans1.simplify()
Out[66]:
$\displaystyle y{\left(x \right)} = \left(C_{1} e^{x \sqrt{b^{2} -c}} + C_{2} e^{-x \sqrt{b^{2} -c}}\right) e^{-b x}$

ひとつだけの解ですますには…

以下の表現ひとつで3つの場合の全てに対応した解になることを示しておきます。

In [67]:
def y(mu, x):
    return exp(-b*x)*(A*cos(mu*x) + B/mu * sin(mu*x))
y(mu, x)
Out[67]:
$\displaystyle \left(A \cos{\left(\mu x \right)} + \frac{B \sin{\left(\mu x \right)}}{\mu}\right) e^{-b x}$

$c-b^2 > 0$ の場合は $\mu \equiv \sqrt{c-b^2}$ として…

In [68]:
y(sqrt(c - b**2), x)
Out[68]:
$\displaystyle \left(A \cos{\left(x \sqrt{-b^{2} + c} \right)} + \frac{B \sin{\left(x \sqrt{-b^{2} + c} \right)}}{\sqrt{-b^{2} + c}}\right) e^{-b x}$

$c-b^2 < 0$ の場合は $\mu \equiv \sqrt{c-b^2} = i \sqrt{b^2-c}$ として…

In [69]:
y(I*sqrt(b**2-c), x)
Out[69]:
$\displaystyle \left(A \cosh{\left(x \sqrt{b^{2} -c} \right)} + \frac{B \sinh{\left(x \sqrt{b^{2} -c} \right)}}{\sqrt{b^{2} -c}}\right) e^{-b x}$

最後に $c-b^2 = 0$ の場合は $\mu \rightarrow 0$ の極限をとって…

In [70]:
limit(y(mu, x), mu, 0)
Out[70]:
$\displaystyle \left(A + B x\right) e^{-b x}$

定数係数2階線形非同次方程式

人力で解く際には，同次方程式と非同次方程式とでは，解く手間がずいぶん違ったが，SymPy ではどちらも同じ。dsolve() を使う。

例題

非同次方程式 $y” + a^2 y = \sin b x$ を解く。

In [71]:
y = Function('y')
var('x')
var('a b', positive = True)

eq = Eq(Derivative(y(x), x, 2) + a**2 * y(x), sin(b*x))
eq
Out[71]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \sin{\left(b x \right)}$
In [72]:
dsolve(eq, y(x))
Out[72]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(a x \right)} + C_{2} \cos{\left(a x \right)} + \frac{\sin{\left(b x \right)}}{a^{2} -b^{2}}$

$a = b$ の場合は上記の解の分母がゼロになってしまうので，別途計算する必要がある。

この状況は，力学では固有振動数 $a$
に等しい振動数の外力が加えられた時に起こる「共鳴（共振）」と呼ばれる現象である。

In [73]:
b = a
eq2 = Eq(Derivative(y(x), x, 2) + a**2 * y(x), sin(b*x))
eq2
Out[73]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \sin{\left(a x \right)}$
In [74]:
dsolve(eq2, y(x))
Out[74]:
$\displaystyle y{\left(x \right)} = C_{2} \sin{\left(a x \right)} + \left(C_{1} -\frac{x}{2 a}\right) \cos{\left(a x \right)}$

積分定数 $C_1, C_2 $がつくのは同次方程式の一般解。

非同次方程式の特殊解
$\displaystyle -\frac{x}{2a} \cos (a x)$ は，振幅が $x$ に比例して単調増加していく。

あえて非同次方程式のロンスキアンを使った公式で特殊解を求める

非同次方程式も dsolve() で解がもとまったわけだが，そこをあえてロンスキアンを使った公式で特殊解を求めてみる。

まず，同次方程式は eq2 の右辺をゼロとした式であるから…

In [75]:
eq0 = Eq(eq2.lhs, 0)
eq0
Out[75]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 0$

微分方程式 eq0 を y(x) について解くと…

In [76]:
dsolve(eq0, y(x))
Out[76]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(a x \right)} + C_{2} \cos{\left(a x \right)}$

上記のように同次方程式の1次独立な基本解はそれぞれ…

In [77]:
y1 = cos(a*x)
y2 = sin(a*x)

ロンスキアンは $W(x) \equiv y_1 y_2′ -y_1′ y_2 = \cdots$

In [78]:
W = (y1*diff(y2, x) - diff(y1, x)*y2).simplify()
W
Out[78]:
$\displaystyle a$

非同次方程式の特殊解 $y_s(x)$ は，非同次項を $R(x) = \sin b x$ として

$$y_s(x) = y_2(x) \int^x \frac{R(x’) y_1 (x’)}{W(x’)} dx’ –y_1(x) \int^x \frac{R(x’) y_2(x’)}{W(x’)} dx’$$

In [79]:
var('a b', positive = True)

R = sin(b*x)

ys = ( y2 * integrate(R*y1/W, x) 
      -y1 * integrate(R*y2/W, x))

ys.simplify().expand()
Out[79]:
$\displaystyle \begin{cases} \frac{\sin{\left(b x \right)}}{a^{2} -b^{2}} & \text{for}\: a \neq b \\-\frac{x \cos{\left(a x \right)}}{2 a} + \frac{\sin{\left(b x \right)} \cos{\left(a x -b x \right)}}{2 a b} & \text{otherwise} \end{cases}$

上記の otherwise の場合の答えが今ひとつです。直接 $a = b$ とすれば…

In [80]:
R = sin(a*x)

ys = ( y2 * integrate(R*y1/W, x) 
      -y1 * integrate(R*y2/W, x))

ys.simplify().expand()
Out[80]:
$\displaystyle -\frac{x \cos{\left(a x \right)}}{2 a} + \frac{\sin{\left(a x \right)}}{2 a^{2}}$

$a = b$ の場合の非同次方程式の特殊解を求めたはずなのに，$\sin (a x)$ に比例する項も出てきました。これは，同次方程式の一般解に吸収できる項でしたね。

ということで，ロンスキアンを使った公式で求めた特殊解は

$\displaystyle -\frac{x \cos (ax)}{2a}$ であることを確認できました。

非同次方程式の一般解

最終的に，この非同次方程式の一般解は，同次方程式の一般解（つまり，同次方程式の基本解 $y_1, \, y_2$ の線形和 $C_1 y_1 + C_2 y_2$）に，非同次方程式の特殊解 $y_s$ を足したものだから，

In [81]:
# 非同次方程式の一般解は...
dsolve(eq0, y(x)).rhs - x*cos(a*x)/(2*a)
Out[81]:
$\displaystyle C_{1} \sin{\left(a x \right)} + C_{2} \cos{\left(a x \right)} -\frac{x \cos{\left(a x \right)}}{2 a}$

注意：dsolve(eq0, y(x)) のままでは，以下のように等式を出力するので，欲しいのは右辺だから .rhs で右辺（right-hand-side）を取り出しています。

In [82]:
dsolve(eq0, y(x))
Out[82]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(a x \right)} + C_{2} \cos{\left(a x \right)}$