Python で数値解析するためのライブラリである SciPy の使用例を示します。グラフは Matplotlib の pyplot.plot() を使って描きます。三角関数などの数学関数については，SymPy を使わず，全て NumPy の関数を使います。
ライブラリの import と設定
# NumPy を使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt
# グラフの目盛設定用
from matplotlib.ticker import MultipleLocator
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
# mathtext font の設定
plt.rcParams['mathtext.fontset'] = 'cm'
SciPy による数値微分
SciPy には数値微分をしてくれる関数 scipy.misc.derivative がありましたが，非推奨となり，SciPy 1.12.0 からは使えなくなったようです。
解析的な微分の定義は（前方差分の極限として）
$$\frac{df}{dx} \equiv \lim_{h \rightarrow 0} \frac{f(x+h) -f(x)}{h}$$
でした。滑らかな関数であれば，（後方差分の極限として）
$$\frac{df}{dx} = \lim_{h \rightarrow 0} \frac{f(x) -f(x-h)}{h}$$
としても良いです。
現実世界では $ h\rightarrow 0$ の極限はとれませんから十分小さい値として h を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。
\begin{eqnarray}
\frac{df}{dx} &\simeq& \frac{1}{2} \left\{\frac{f(x+h) – f(x)}{h} + \frac{f(x) – f(x-h)}{h} \right\} \\
&=& \frac{f(x+h) – f(x-h)}{2h}
\end{eqnarray}
関数 $f(x)$ が初等関数（および初等関数の合成関数）である場合は解析的に微分できますが，そうでない場合，例えば $f(x)$ が数値的にしか与えられていないとか，何か解析的に微分できない場合には，数値微分によって近似的な値を求める必要があります。
また，解析的に微分できそうだが，人力ではけっこう大変そうな場合にも，数値微分によって近似的な値を使うほうが便利な場合があるかもしれません。
以下のように，関数 f(x) を数値微分してくれる関数 ndiff() を定義します。
def ndiff(f, x, h=1e-5):
# h を省略すると h=1e-5 として計算する
return (f(x + h) - f(x - h))/(2*h)
例として，$f(x) = \sin x$ の $x = 3$ における数値微分を求めてみます。（答えは $\cos 3$ です。）
# h を省略した場合
print(' ndiff(f, 3) = %.15f' % ndiff(np.sin, 3))
for h in [1e-4, 1e-5, 1e-6, 1e-7, 1e-8]:
print('h = %3.1e, ndiff = %.15f' % (h, ndiff(np.sin, 3, h)))
print(' '*10 + 'cos(%3.1f) = %.15f' % (3, np.cos(3)))
以上の結果からわかるように，h をやみくもに小さくすればするほど，数値微分の精度があがる，というわけではないです。
○練習：関数の数値微分
$\displaystyle f(x) = \frac{x^3}{e^x -1}$ の極大値を求めるための準備として，まず，$f(x)$ とその数値微分 ndiff(f, x) を $0 \leq x \leq 5$ の範囲でグラフにする。
def f(x):
return x**3/(np.exp(x) - 1)
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# x の範囲
xrange = np.linspace(0.001, 5, 200)
# y = f(x) のグラフ
ax.plot(xrange, f(xrange), label = '$f(x)$')
# y = f'(x) のグラフ
ax.plot(xrange, ndiff(f, xrange), label = '$f^{\prime}(x)$')
# x 軸のラベル
ax.set_xlabel('$x$')
# グリッド
ax.grid()
# 凡例の表示
ax.legend();
上のグラフから，数値微分した微分係数 $f^{\prime}(x)$ の値がゼロになる $x$ は，$2 < x < 3$ のあたりであること，そのあたりで確かに $f(x)$ が極大（最大）となっていることが見てとれます。
SciPy による数値積分
scipy.integrate.quad()
SciPy には数値積分してくれる関数 scipy.integrate.quad() が用意されています。それを使ってみましょう。結果には積分値と推定誤差が出力されます。
$\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ の計算をします。
数値計算の際は，三角関数等は NumPy の関数を使い，円周率 $\pi$ も NumPy の np.pi を使います。
from scipy.integrate import quad
# quad を使うときは変数名 x 決め打ち
def sinsin(x):
return np.sin(np.sin(x))
# quad で数値積分
# quad(f, a, b): 被積分関数 f(x) を x = a から b まで数値積分
quad(sinsin, 0, np.pi)
scipy.integrate.quad() は上記のように，数値積分の値と絶対誤差の推定値をタプル (tuple) で返します。
数値積分値（と誤差）を個別に取り出すには，例えば以下のように。
y, abserr = quad(sinsin, 0, np.pi)
print("積分値は", y)
print("誤差は", abserr)
参考：シンプソン法による数値積分
シンプソン法で $\displaystyle \int_a^b f(x)\, dx$ を求める。
積分区間 $[a, b]$ を $N (=2n)$ 等分（偶数等分）し，
$$h = \frac{b-a}{N}, \quad f_i = f(a + i\,h)$$
とする。シンプソン法の公式は
\begin{eqnarray}
\int_a^b f(x)\, dx &\simeq&
\frac{h}{3}( f_0 + 4 f_1 + f_2) + \frac{h}{3}( f_2 + 4 f_3 + f_4) + \\
&& \frac{h}{3}( f_4 + 4 f_5 + f_6) + \frac{h}{3}( f_6 + 4 f_7 + f_8)+ \cdots\\
&&+ \frac{h}{3}( f_{N-2} + 4 f_{N-1} + f_N)\\
&=&\frac{h}{3}
\left(f_0 + 2 \sum_{i=1}^{n-1} f_{2i} + 4 \sum_{i=1}^{n} f_{2i-1}+ f_{2n}\right)
\end{eqnarray}
SciPy では scipy.integrate.quad() を使って数値積分することができるのでわざわざ自作の関数をつくる必要はないのですが，教育的見地からあえて，シンプソン法の公式を使って数値積分の値を返す関数 simpson() を以下のように定義してみます。
def simpson(f, a, b, Ndiv):
h = (b - a)/Ndiv
simp = 0
for i in range(0, Ndiv, 2):
simp += f(a + i*h) + 4*f(a + (i+1)*h) + f(a + (i+2)*h)
return simp*h/3
N の値を変えて $\displaystyle \int_0^{\pi} \sin(\sin x)\, dx$ を simpson() で数値積分して，精度を確認してみます。
def f(x):
return np.sin(np.sin(x))
for N in [10*2**i for i in range(9)]:
print('N = %4d, %.15f' %(N, simpson(f, 0, np.pi, N)))
以上の結果から，どんなことがわかるかな？
SciPy による方程式の数値解
例としてベッセル関数 $J_0(x)$ のゼロ点を求めてみます。SciPy におけるベッセル関数 $J_0(x)$ は scipy.special.jv(0,x) です。
まずは $J_0(x) = 0$ となる $x$ のあたりをつけるためにグラフを描いてみます。
from scipy.special import jv
def J0(x):
return jv(0,x)
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
xrange = np.linspace(0, 10, 200)
ax.plot(xrange, J0(xrange))
# x 軸
ax.plot(xrange, 0*xrange)
ax.set_xlabel('x')
ax.set_ylabel('J0(x)')
ax.grid()
ax.set_xlim(0, 10)
ax.set_title('ベッセル関数 J0(x)');
scipy.optimize.root_scalar
方程式の数値解を求めてくれる関数は scipy.optimize.root_scalar です。$2 < x < 3$ のあたりで $x$ 軸と交わっていますから，このあたりの解を探します。
from scipy.optimize import root_scalar
# x が 2 ~ 3 の間で j0(x) = 0 の根を探す
root_scalar(J0, bracket = [2, 3])
scipy.optimize.root_scalar() は結果を上記のように表示します。解 root のみを表示したい場合は以下のようにします。
# 根のみを出力するには...
ans = root_scalar(J0, bracket = [2, 3])
print(ans.root)
# ほんとに根になってるか確認。ゼロに近ければよい。
J0(ans.root)
# どうしてももっと精度が欲しいというのであれば...
# xtol に，小さい値をいれてみる
ans = root_scalar(J0, bracket = [2, 3], xtol=1e-15)
ans.root
# ほんとに根になってるか確認。ゼロに近ければよい。
J0(_)
○練習：関数の極大値
- $\displaystyle f(x) = \frac{x^3}{e^x -1}$ が極大値をとるときの $x$ の値，およびそのときの極大値を求める。
- $\displaystyle g(x) = \frac{x^5}{e^x -1}$ が極大値をとるときの $x$ の値，およびそのときの極大値も同様にして求めてみよ。
まず $f(x)$ の数値微分 ndiff(f, x) をグラフにしてみる。
def f(x):
return x**3/(np.exp(x) - 1)
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
xrange = np.linspace(0, 5, 200)
ax.plot(xrange, ndiff(f, xrange), label = '$f^{\prime}(x)$')
# x 軸
ax.plot(xrange, 0*xrange)
ax.set_xlabel('$x$')
ax.set_ylabel(' ')
ax.grid()
ax.legend();
ndiff(f, x) は $2 < x <3$ のあたりでゼロになっているので，そのあたりを root_scalar() で探してみる。
def df(x):
return ndiff(f, x)
ans = root_scalar(df, bracket = [2, 3])
ans.root
# そのときの f(x) の極大値は...
f(_)
SciPy による1階常微分方程式の数値解法
scipy.integrate.solve_ivp()
SciPy には常微分方程式を数値的に解いてくれる関数 scipy.integrate.solve_ivp() が用意されています。それを使ってみましょう。
例として，簡単なロジスティック方程式
$$ \frac{dx}{dt} = f(x) = (1-x)\,x$$
を，初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 0.1$，$t = t_0 = 0$ から $t_1 = 10$ まで解きます。
solve_ivp() を使う場合は変数名が決め打ちで t, y です。連立常微分方程式の場合は y がリストになります。
from scipy.integrate import solve_ivp
# scipy.integrate.solve_ivp() は
# dy/dt = f(t, y) の形を解く
# solve_ivp を使う際は変数名 t, y は決め打ち
# solve.ivp() に渡す関数の引数の順番に注意。t が先
# 連立でないときは，1次元リストとして返す
def f(t, y):
x = y[0]
return [(1-x)*x]
# t の範囲
tini = 0
tend = 10
# 計算する点は t の範囲を N 等分して N+1 点
N = 20
t_list = np.linspace(tini, tend, N+1)
# y の初期値
y_ini = [0.1]
answer = solve_ivp(f, [tini, tend], y_ini, t_eval = t_list,
# デフォルトの許容誤差は大きめなので
# 小さい値に設定
rtol = 1.e-12,
atol = 1.e-15)
print(answer)
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
ax.scatter(answer.t, answer.y[0], s = 10)
ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel("$x(t)$")
ax.set_title("ロジスティック方程式の数値解");
参考までに，解析解は
$$x(t) = \frac{\exp(t)}{9 + \exp(t)}$$
です。（解析解がわかっているときには，わざわざ数値計算する必要はないのですが，練習ということで。）
# 解析解
def xsol(t):
return np.exp(t)/(9+np.exp(t))
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# 解析解
ts = np.linspace(tini, tend, 200)
ax.plot(ts, xsol(ts), c = 'tab:blue', label='解析解')
# 数値解
ax.scatter(answer.t, answer.y[0], zorder = 2,
c = 'tab:orange', s = 10, label = '数値解')
ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel("$x(t)$")
ax.set_title("ロジスティック方程式の解")
ax.legend();
SciPy による2階常微分方程式の数値解法
scipy.integrate.solve_ivp()
SciPy には常微分方程式を数値的に解いてくれる関数 scipy.integrate.solve_ivp() が用意されています。この函数は，1階常微分方程式だけでなく，連立1階常微分方程式も，従って高階の常微分方程式も解いてくれるので，それを使ってみましょう。
例として，簡単な減衰＋強制振動の方程式
$$ \frac{d^2 x}{dt^2} = -x – a \frac{dx}{dt} + b \cos(t) $$
を，
- $a , b $ にいくつか値を入れて
- 初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 3, v_0 = v(t_0) = 0$，
- $t = t_0 = 0$ から $t_1 = 20$ までを解きます。
2階常微分方程式は，
\begin{eqnarray}
\frac{dx}{dt} &=& F_1(t,x,v) = v \\
\frac{dv}{dt} &=& F_2(t,x,v) = -x – a v + b \cos t
\end{eqnarray}
として連立1階微分方程式の形にして解きます。
from scipy.integrate import solve_ivp
# scipy.integrate.solve_ivp() は
# dy/dt = F(t, y) の形を解く
def F(t, y, a, b):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
x = y[0]
v = y[1]
F1 = v
F2 = -x - a*v + b*np.cos(t)
return [F1, F2]
tini = 0
tend = 20
t_list = np.linspace(tini, tend, 200)
y_ini = [3, 0]
a = 0.5; b = 0.2
answer1 = solve_ivp(F, [tini, tend], y_ini, t_eval = t_list,
args=(a, b),
# デフォルトの許容誤差は大きめなので
# 小さい値に設定
rtol = 1.e-12,
atol = 1.e-15)
print(answer1)
横軸に $t = $ answer1.t，縦軸に $x(t) = $ answer1.y[0] をとってグラフにしてみます。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# 数値解 x は .y[0] のようにして取り出す
# 数式を label に設定する際は $ で囲んで LaTeX 式に
# a, b の値を凡例に
key = '$a=%3.1f, b=%3.1f$' % (a, b)
ax.plot(answer1.t, answer1.y[0], label=key)
ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel("$x(t)$")
ax.legend()
ax.set_title("減衰＋強制振動")
ax.set_xlim(0, 20);
○練習：減衰＋強制振動
減衰＋強制振動の方程式を以下の場合に解き，グラフにする。
まず $a$ の影響を調べるため，
- $a = 0, b = 0$
- $a = 0.2, b = 0$
- $a = 0.5, b = 0$
- $a = 0.8, b = 0$
次に，$b$ の影響を調べるため，
- $a = 0, b = 0$
- $a = 0, b = 0.2$
- $a = 0, b = 0.5$
- $a = 0, b = 0.8$
○練習：単振り子
2階常微分方程式
$$ \frac{d^2\theta}{dt^2} = -4\pi^2 \sin\theta$$
を初期条件
$$ \theta = 45^{\circ}, \ \ \frac{d\theta}{dt} = 0 \quad \mbox{at}\ \ t = 0$$
のもとで $t=0$ から $t = 5$ まで数値的に解け。
$\theta \Rightarrow $ y[0]，$\dot{\theta} \Rightarrow $ y[1] として
def F(t, y):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
theta = y[0]
dottheta = y[1]
F1 = dottheta
F2 = -4*np.pi**2*np.sin(theta)
return [F1, F2]
tini = 0
tend = 5
t_list = np.linspace(tini, tend, 200)
# 初期条件はラジアンになおして
y_ini = [np.radians(45), 0]
answer2 = solve_ivp(F, [tini, tend], y_ini, t_eval = t_list,
# デフォルトの許容誤差は大きめなので
# 小さい値に設定
rtol = 1.e-12,
atol = 1.e-15)
一方，$\theta$ があまり大きくないとすれば，テイラー展開の1次の近似で $\sin\theta \simeq \theta$ であるから
$$ \frac{d^2\theta}{dt^2} \simeq -4\pi^2 \theta$$
$\simeq$ を $=$ に変えて方程式にすると，これは単振動の方程式であるから，同じ初期条件で以下のように解析的に解ける。
$$\theta_s = \theta_0 \cos (2 \pi t) = \frac{\pi}{4} \cos (2 \pi t)$$
また，このとき
$$\dot{\theta}_s = – \frac{\ \pi^2}{2} \sin (2 \pi t)$$
この単振動の解も重ねてグラフにしてみよ。周期にはどのような違いがあるか。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# 数値解 x は .y[0] のようにして取り出す
# 縦軸の値をラジアンから度に変換
theta = np.degrees(1)*(answer2.y[0])
# 数値解
ax.plot(answer2.t, theta, label='数値解')
# 同じ初期条件の単振動解
ax.plot(t_list, 45*np.cos(2*np.pi*t_list), label = '単振動解')
ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel(r"$\theta$ (度)")
ax.legend();
参考：ルンゲ・クッタ法による2階常微分方程式の数値解法
さて，SymPy で自作のルンゲ・クッタ法 rk2() を使った数値解法の際は，かなり時間がかかりました。
NumPy の関数を使った場合はどうでしょうか。
def rk2(Func, x0, v0, t0, t1, N):
""" d^2 x/dt^2 = Func(t, x, v) を
初期条件 t = t0 で x = x0, v = v0 とし，
[t0, t1] を N 分割して
4次のルンゲ・クッタ法で解く
"""
h = (t1 - t0)/N
T = [t0]
X = [x0]
V = [v0]
t = t0
x = x0
v = v0
for i in range(N):
k1 = h*v
m1 = h*Func(t, x, v)
k2 = h*(v+m1/2)
m2 = h*Func(t+h/2, x+k1/2, v+m1/2)
k3 = h*(v+m2/2)
m3 = h*Func(t+h/2, x+k2/2, v+m2/2)
k4 = h*(v+m3)
m4 = h*Func(t+h, x+k3, v+m3)
x = x + (k1 + 2*k2 + 2*k3 + k4)/6
v = v + (m1 + 2*m2 + 2*m3 + m4)/6
t = t + h
T.append(t)
X.append(x)
V.append(v)
return [T, X, V]
def F(t, x, v):
# NumPy の関数・定数を使う
return -4*np.pi**2*np.sin(x)
x0 = np.radians(45)
v0 = 0
t0 = 0
t1 = 3
for N in [100, 200, 400, 800, 1600, 3200]:
T2, X2, V2 = rk2(F, x0, v0, t0, t1, N)
# X2 の最後の要素 X2[-1] を表示
print("N =%5d, x(t1) =%15.12f" % (N, X2[-1]))
%%time
T2, X2, V2 = rk2(F, x0, v0, t0, t1, 1600)
rk2() は同じコードですので，F(t, x, v) に NumPy の関数を使っただけで，約 600 倍 (!) 程度高速になっていることがわかります。
ですので，解析的に計算をしたいなら SymPy ですが，数値計算を高速にしたいなら NumPy，という具合に使い分けが肝心です。
○練習：単振り子の周期
数値計算で求めた単振り子の周期を調べるため，$\dot{\theta} =$ y[1] のグラフを描いてみます。比較のために，単振動の場合の $\dot{\theta}_s$ も一緒にグラフにします。
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
# 数値解
ax.plot(answer2.t, answer2.y[1], label='数値解')
# 同じ初期条件の単振動解
ax.plot(t_list, -np.pi**2/2*np.sin(2*np.pi*t_list), label = '単振動解')
ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.set_xlim(0, 5)
ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel(r"$\dot{\theta}$")
ax.legend();
上のグラフから，単振動の場合は $t = 0$ で $\dot{\theta} = 0$ から始まり，$t = 0.5$ で折り返し，$t = 1$ で再び $\dot{\theta} = 0$ となることから，周期が $1$ であることがわかります。
単振り子の数値解の場合は $t = 1$ ではなく，$1 < t_1 < 1.5$ にあることがわかります。この値を数値的に求めてください。
ヒント：以下のように，x を入れると $t = $ 0 から $t = $ x まで数値的に微分方程式を解き，最後の値 $\dot{\theta}(x)$ を返す関数を dottheta(x) として定義してみます。
def F(t, y):
theta = y[0]
dottheta = y[1]
F1 = dottheta
F2 = -4*np.pi**2*np.sin(theta)
return [F1, F2]
def dottheta(x):
tini = 0
tend = x
t_list = np.linspace(tini, tend)
# 初期条件 theta0 はラジアンになおして
y_ini = [np.radians(theta0), 0]
answer2 = solve_ivp(F, [tini, tend], y_ini, t_eval = t_list,
# デフォルトの許容誤差は大きめなので
# 小さい値に設定
rtol = 1.e-12,
atol = 1.e-15)
return (answer2.y[1])[-1]
theta0 = 45
ans1 = root_scalar(dottheta, bracket = [1, 1.5])
t1 = ans1.root
print(t1)
念のため，次に「ひとまわり」して戻ってくる時刻は $2 < t_2 < 2.5$ にあります。$t_2$ を求め，時間間隔 $t_2 – t_1$ と $t_1$ が等しいかどうか，確認してください。
ans2 = root_scalar(dottheta, bracket = [2, 2.5])
t2 = ans2.root
print('t1 =', t1)
print('t2 - t1 =', t2 - t1)
念には念をいれて，$0 < t < 5$ の範囲で $\dot{\theta}(t)$ が正の側からゼロになる時刻 $t_1, t_2, t_3, \dots$ を求め，時間間隔が一定であるか（つまり，単振り子の運動が確かに周期的であるか）確認してください。
t0 = [0]
for i in range(1, 5):
ans = root_scalar(dottheta, bracket = [i, i+0.5], xtol=1e-15)
t0.append(ans.root)
print('t%1d - t%1d = %15.13f' %(i, i-1, t0[i]-t0[i-1]))
時間間隔が（ほぼ）一定であることが確認できたら，これを $t=0$ における初期条件 $\theta_0 = 45^{\circ}, \ \dot{\theta}_0 = 0$ の場合の周期とします。
同様にして，$t=0$ における初期条件
- $\theta_0 = 1^{\circ}, \ \dot{\theta}_0 = 0$ （単振動解との比較も兼ねて）
- $\theta_0 = 60^{\circ}, \ \dot{\theta}_0 = 0$
- $\theta_0 = 70^{\circ}, \ \dot{\theta}_0 = 0$
- $\theta_0 = 80^{\circ}, \ \dot{\theta}_0 = 0$
の場合の単振り子の周期を求めなさい。この結果から，初期条件としての振幅 $\theta_0$ と周期の間にはどのような関係があるか，考察してください。