Return to Python でコンピュータ演習

SymPy で常微分方程式

「理工系の数学C」の「SymPy で常微分方程式」の抜粋バージョン。詳細は以下のページをご覧ください。

常微分方程式の解法 dsolve()

Python では，SymPy の dsolve() を使って1階および2階の常微分方程式を解くことができます。

まず必要なパッケージを import します。

In [1]:
# a Python library for symbolic mathematics
from sympy import *
# 1文字変数の定義が省略できる
from sympy.abc import *
# π，ネイピア数，虚数単位
from sympy import pi, E, I

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

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

1階線形常微分方程式の例

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

In [2]:
# y を関数として宣言
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)}$
In [3]:
dsolve(eq, y(x))
Out[3]:
$\displaystyle y{\left(x \right)} = C_{1} e^{- x^{2}}$

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

マルサスの人口モデル

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

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

In [4]:
# ギリシャ文字も from sympy.abc import * で定義済み
# var('gamma')
N = Function('N')

eqm = Eq(Derivative(N(t), t), gamma*N(t))
eqm
Out[4]:
$\displaystyle \frac{d}{d t} N{\left(t \right)} = \gamma N{\left(t \right)}$
In [5]:
# 初期条件変数は2文字変数だから定義しよう
var('t0, N0')

dsolve(eqm, N(t), ics={N(t).subs(t,t0):N0})
Out[5]:
$\displaystyle N{\left(t \right)} = N_{0} e^{\gamma t} e^{- \gamma t_{0}}$
参考：米国の人口とマルサスモデル

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

In [6]:
usa = [
# 西暦, 人口（百万人）
  [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 [7]:
t0 = usa[0][0]
N0 = usa[0][1]

残りのパラメーター $\gamma_1$ は，別の時刻 $t_1$ における $N_1 = N_m(t_1)$ から求める。

\begin{eqnarray}
N_1 &=& N_0 e^{\gamma_1 (t_1 – t_0)} \\
\log \frac{N_1}{N_0} &=& \gamma_1 (t_1 – t_0) \\
\therefore\ \ \gamma_1 &=& \frac{1}{t_1 – t_0} \log \frac{N_1}{N_0}
\end{eqnarray}

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

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

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

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

In [9]:
def Nm(t):
    return N0 * exp(gamma1*(t-t0))

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

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

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

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

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

一方，リスト usa には

usa = [
 [x1, y1],
 [x2, y2],
 ...,
 [xn, yn]
]

のように数値データが格納されているため，$x$ 座標のみ，$y$ 座標のみのリストを作成するために少しだけ工夫が必要です。やり方を2つほど…

  1. for 文を使ってリストを作成する例。
  2. numpy.array() で配列に変換して「列」成分を指定する例。
In [11]:
# SPB で2次元リスト usa を plot_list() する例 1
# for 文で x 成分リスト usaX と y 成分リスト usaY を作成
usaX = []
usaY = []
for dat in usa:
    usaX.append(dat[0])
    usaY.append(dat[1])

p2 = plot_list(usaX, usaY, "アメリカの人口", 
               # 線で繋がす点で, 凡例を表示
               is_point = True, legend = True, 
               xlabel = "年", ylabel = "人口（単位：百万人）");

In [12]:
# SPB で2次元リスト usa を plot_list() する例 2
# NumPy 配列に変換して1列目 [:,0] 2列目 [:,1] を指定
import numpy as np
npusa = np.array(usa)

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

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

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

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

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

In [14]:
N = Function('N')
var('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}}}$
参考：米国の人口とヴェアフルストモデル

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

$t_1$（年）と $t_2$（年）の値を使って $\gamma$ と $N_{\rm max}$ を求めます。指数関数を含む連立方程式はなかなか解いてくれないかもしれないので，簡単な代数方程式の形にして解きます。

$N(t)$ の分子分母を $e^{\gamma t}$ で割り，さらに

\begin{eqnarray}
n_0 &\equiv& \frac{N_0}{N_{\rm max}}
\end{eqnarray}

とします。

In [16]:
t0 = usa[0][0]
N0 = usa[0][1]
var('n0 gamma2')

def Nv1(t):
    return N0/(n0 + (1-n0)*exp(gamma2*(t0-t)))
Nv1(t)
Out[16]:
$\displaystyle \frac{3.9}{n_{0} + \left(1 – n_{0}\right) e^{\gamma_{2} \cdot \left(1790 – t\right)}}$

$t_1 = 1850$（年）と $t_2 = 1910$（年）の値を入れて連立方程式の形にします。簡単な連立方程式にするため，さらに $T \equiv e^{-60\, \gamma_2}$ とし，$n_0$ と $T$ に関するシンプルな連立方程式の形にして，solve() で解きます。

In [17]:
t1 = usa[6][0];
eq1 = Eq(usa[6][1], Nv1(t1)) 

# var('T') は省略可
eq1 = eq1.subs(exp(-60*gamma2), T)
eq1
Out[17]:
$\displaystyle 23.2 = \frac{3.9}{T \left(1 – n_{0}\right) + n_{0}}$
In [18]:
t2 = usa[12][0];
eq2 = Eq(usa[12][1], Nv1(t2))
eq2 = eq2.subs(exp(-60*gamma2), T)
eq2
Out[18]:
$\displaystyle 92.0 = \frac{3.9}{T^{2} \cdot \left(1 – n_{0}\right) + n_{0}}$
In [19]:
ans = solve([eq1, eq2], [n0, T])
ans
Out[19]:
[(0.0200125277046207, 0.151115116017121)]

$n_0, \ T$ からもとの $N_{\rm max}, \ \gamma_2$ の値になおすと…

In [20]:
n0, T = ans[0]
Nmax = N0/n0
gamma2 = -1/60*log(T)

def NV(t):
    return N0*Nmax*exp(gamma2*t)/(
        N0*exp(gamma2*t)-(N0-Nmax)*exp(gamma2*t0))
In [21]:
p5 = plot(NV(t), (t, 1790, 1940), "ヴェアフルストモデル", 
          legend = True, 
          xlabel = "年", ylabel = "人口（単位：百万人）");

In [22]:
p6 = p1+p2+p5
p6.show();

最も簡単な2階線形常微分方程式の例

$y^{\prime\prime} + K y = 0$ あるいは移項して $y^{\prime\prime} = – K y$ を解く。

$K > 0$ の場合
In [23]:
y = Function('y')
# var('x')
var('K', positive = True)

eq1 = Eq(Derivative(y(x), x, 2), -K*y(x))
eq1
Out[23]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = – K y{\left(x \right)}$
In [24]:
dsolve(eq1, y(x))
Out[24]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(\sqrt{K} x \right)} + C_{2} \cos{\left(\sqrt{K} x \right)}$
$K < 0$ の場合
In [25]:
var('K', negative = True)

eq2 = Eq(Derivative(y(x), x, 2), -K*y(x))
eq2
Out[25]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = – K y{\left(x \right)}$
In [26]:
dsolve(eq2, y(x))
Out[26]:
$\displaystyle y{\left(x \right)} = C_{1} e^{- x \sqrt{- K}} + C_{2} e^{x \sqrt{- K}}$
$K = 0$ の場合
In [27]:
eq3 = Eq(Derivative(y(x), x, 2), 0)
eq3
Out[27]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 0$
In [28]:
dsolve(eq3, y(x))
Out[28]:
$\displaystyle y{\left(x \right)} = C_{1} + C_{2} x$

定数係数2階線形常微分方程式

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

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

$b, c$ は定数。一般解 $y$ は以下のように場合分けして…

In [29]:
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[29]:
$\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 [30]:
ans1 = dsolve(eq, y(x))
ans1.simplify()
Out[30]:
$\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}$

非同次方程式の例

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

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

In [31]:
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[31]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \sin{\left(b x \right)}$
In [32]:
dsolve(eq, y(x))
Out[32]:
$\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 [33]:
b = a
eq2 = Eq(Derivative(y(x), x, 2) + a**2 * y(x), sin(b*x))
eq2
Out[33]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \sin{\left(a x \right)}$
In [34]:
dsolve(eq2, y(x))
Out[34]:
$\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$ に比例して単調増加していく。

練習 2-1

重力場中の投げ上げ運動

  1. 高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y$ を求めなさい。運動方程式は，

$$ \frac{d^2 y}{dt^2} = -g$$

  1. 速度に比例する空気抵抗がある場合，運動方程式は以下のようになる。
    このとき，高さ $0$ から初速度 $v_0$ で鉛直上方に投げ上げた物体の時刻 $t$ における高さ $y$ を
    求めなさい。

$$ \frac{d^2 y}{dt^2} = -g – \beta \frac{dy}{dt}$$

  1. 前問 2. で求めた解が $\beta \rightarrow 0$ のとき，前問 1. の答えに一致することを示しなさい。
In [35]:
# 1. 
y = Function('y')
# var('t')
var('g', positive = True)

eq = Eq(Derivative(y(t), t, 2), -g)
eq
Out[35]:
$\displaystyle \frac{d^{2}}{d t^{2}} y{\left(t \right)} = – g$
In [36]:
dsolve(eq, y(t))
Out[36]:
$\displaystyle y{\left(t \right)} = C_{1} + C_{2} t – \frac{g t^{2}}{2}$
In [37]:
# 初期条件を与えて解く例
# 1文字変数以外を使うときは var() で宣言しておく
var('v0')

dsolve(eq, y(t), 
       ics = {y(t).subs(t, 0):0,          # y(0) = 0
              diff(y(t), t).subs(t, 0):v0 # dy/dt(0) = v0
             }
      )
Out[37]:
$\displaystyle y{\left(t \right)} = – \frac{g t^{2}}{2} + t v_{0}$
In [38]:
# 2. も同様に
In [39]:
# 3. 
# 等式 eq の右辺は eq.rhs
# 極限の例

limit(sin(x)/x, x, 0)
Out[39]:
$\displaystyle 1$
In [ ]: