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 からは使えなくなったようです。
解析的な微分の定義は(前方差分の極限として)
でした。滑らかな関数であれば,(後方差分の極限として)
としても良いです。
現実世界では
関数
また,解析的に微分できそうだが,人力ではけっこう大変そうな場合にも,数値微分によって近似的な値を使うほうが便利な場合があるかもしれません。
以下のように,関数 f(x)
を数値微分してくれる関数 ndiff()
を定義します。
def ndiff(f, x, h=1e-5):
# h を省略すると h=1e-5 として計算する
return (f(x + h) - f(x - h))/(2*h)
例として,
# 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
をやみくもに小さくすればするほど,数値微分の精度があがる,というわけではないです。
○練習:関数の数値微分
ndiff(f, x)
を
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();
上のグラフから,数値微分した微分係数
SciPy による数値積分
scipy.integrate.quad()
SciPy には数値積分してくれる関数 scipy.integrate.quad()
が用意されています。それを使ってみましょう。結果には積分値と推定誤差が出力されます。
数値計算の際は,三角関数等は 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)
参考:シンプソン法による数値積分
シンプソン法で
積分区間
とする。シンプソン法の公式は
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
の値を変えて 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 による方程式の数値解
例としてベッセル関数 SciPy
におけるベッセル関数 scipy.special.jv(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
です。
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(_)
○練習:関数の極大値
が極大値をとるときの の値,およびそのときの極大値を求める。 が極大値をとるときの の値,およびそのときの極大値も同様にして求めてみよ。
まず 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)
は 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()
が用意されています。それを使ってみましょう。
例として,簡単なロジスティック方程式
を,初期条件
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("ロジスティック方程式の数値解");
参考までに,解析解は
です。(解析解がわかっているときには,わざわざ数値計算する必要はないのですが,練習ということで。)
# 解析解
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階常微分方程式も,従って高階の常微分方程式も解いてくれるので,それを使ってみましょう。
例として,簡単な減衰+強制振動の方程式
を,
にいくつか値を入れて- 初期条件
で , から までを解きます。
2階常微分方程式は,
として連立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)
横軸に answer1.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);
○練習:減衰+強制振動
減衰+強制振動の方程式を以下の場合に解き,グラフにする。
まず
次に,
○練習:単振り子
2階常微分方程式
を初期条件
のもとで
y[0]
,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)
一方,
また,このとき
この単振動の解も重ねてグラフにしてみよ。周期にはどのような違いがあるか。
# 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,という具合に使い分けが肝心です。
○練習:単振り子の周期
数値計算で求めた単振り子の周期を調べるため,y[1]
のグラフを描いてみます。比較のために,単振動の場合の
# 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();
上のグラフから,単振動の場合は
単振り子の数値解の場合は
ヒント:以下のように,x
を入れると 0
から 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)
念のため,次に「ひとまわり」して戻ってくる時刻は
ans2 = root_scalar(dottheta, bracket = [2, 2.5])
t2 = ans2.root
print('t1 =', t1)
print('t2 - t1 =', t2 - t1)
念には念をいれて,
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]))
時間間隔が(ほぼ)一定であることが確認できたら,これを
同様にして,
(単振動解との比較も兼ねて)
の場合の単振り子の周期を求めなさい。この結果から,初期条件としての振幅