SciPy で数値解析

# NumPy も使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']

from scipy.integrate import quad

# quad を使うときは変数名 x 決め打ち
def sinsin(x):
    return np.sin(np.sin(x))

# quad で数値積分
quadans = quad(sinsin, 0, np.pi)
print(quadans)

# 積分値のみを出力するには...
print("積分値は", quadans[0])

# 別のやり方
ans, err = quadans
print("積分値は", ans)
print("推定誤差は", err)

from scipy.special import jv

def J0(x):
    return jv(0,x)

# pyplot.plot でグラフを描く
x = np.linspace(0, 10, 101)
plt.plot(x, J0(x))

# pyplot.plot で x 軸を描く
plt.plot(x, 0*x, color='black')

plt.xlabel('x')
plt.ylabel('J0(x)')
plt.grid()
plt.xlim(0, 10)
plt.title('ベッセル関数 J0(x)');

from scipy.optimize import root_scalar

# x が 2 ~ 3 の間で j0(x) = 0 の根を探す
ans = root_scalar(J0, bracket = [2, 3])
print(ans)

# 根のみを出力するには...
print(ans.root)

# ほんとに根になってるか確認。ゼロに近ければよい。
J0(ans.root)

# どうしてももっと精度が欲しいというのであれば...
# xtol に，小さい値をいれてみる
ans = root_scalar(J0, bracket = [2, 3], xtol=1e-14)
ans.root

# ほんとに根になってるか確認。ゼロに近ければよい。
J0(_)

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-14)
print(answer)

# グラフにするために，t の値と y の値を取り出す
# t は answer.t, y は answer.y[0] として取り出す
# scatter では点のサイズは s=10 などと設定
plt.scatter(answer.t, answer.y[0], s=10)
plt.grid()
plt.xlabel("$t$")
plt.ylabel("$x(t)$")
plt.title("ロジスティック方程式の数値解");

# 解析解
def x(t):
    return np.exp(t)/(9+np.exp(t))

# 解析解のラインも描く
t = np.linspace(tini, tend, 201)
plt.plot(t, x(t), color='red', linewidth=1, label='解析解')

plt.scatter(answer.t, answer.y[0], s=10, label='数値解')
plt.grid()
plt.xlabel("$t$")
plt.ylabel("$x(t)$")
plt.title("ロジスティック方程式の数値解")
plt.legend();

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, 201)
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-14)

print(answer1)

# 数値解 x は .y[0] のようにして取り出す
# 数式を label に設定する際は $ で囲んで LaTeX 式に
# a, b の値を凡例に
key = '$a=%3.1f, b=%3.1f$' % (a, b)
plt.plot(answer1.t, answer1.y[0], label=key)
plt.grid()
plt.xlabel("$t$")
plt.ylabel("$x(t)$")
plt.legend()
plt.title("減衰＋強制振動")
plt.xlim(0, 20);