Python の SciPy で数値解析

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

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

# 横軸 x 座標の値と... 
x = [0, 1, 2, 3, 4, 5]
# 対応する縦軸 y 座標の値を準備する
y = [0**2, 1**2, 2**2, 3**2, 4**2, 5**2]

plt.plot(x, y);

# 線でつながなず，点々だけを描く例
plt.plot(x, y, linestyle='', marker='o');

# 線でつながなず，点々だけを描くもう一つの例
plt.scatter(x, y);

# 関数 y = sin(x) を -2π から 2π まで描く例
# π は numpy で定義されている定数 np.pi を使う
# -2π から 2π の範囲を 100 等分した点を生成する
x = np.linspace(-2*np.pi, 2*np.pi, 100)

# sin も numpy で定義されている関数 np.sin() を使う
plt.plot(x, np.sin(x));

from scipy.misc import derivative

def f(x):
    return np.sin(x)

# 数値微分の値
print("数値微分", derivative(f, 3, dx=1E-6))

# 解析的微分の値
print("解析解　", np.cos(3))

def f(x):
    return np.sin(x)

def df(x, h):
    return derivative(f, x, dx=h)

x = np.linspace(-2*np.pi, 2*np.pi, 101)
plt.plot(x, np.cos(x) - df(x, 1E-6));

plt.plot(x, np.cos(x) - df(x, 1E-7));

from scipy.integrate import quad

def sinsin(x):
    return np.sin(np.sin(x))

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

# 積分値のみを取り出す
print(ans[0])

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])
ans.root

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

# どうしてももっと精度が欲しいというのであれば...
# 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 = fun(t, y) の形を解く

def fun(t, y):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
# 連立でないときは，1次元リストとして返す
    return [(1-y)*y]

t0 = 0
t1 = 10
t_span = [t0, t1]
y_ini = [0.1]

t_list = np.linspace(t0, t1, 21)
answer = solve_ivp(fun, t_span, y_ini, 
                   t_eval = t_list, 
                   rtol = 1.e-12, 
                   atol = 1.e-14)

# scatter では点のサイズは s=10 などと設定
# 計算結果は .y として取り出します。
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(t0, t1, 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 = Fun(t, y) の形を解く

def Fun(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]

t0 = 0
t1 = 20
t_span = [t0, t1]
y_ini = [3, 0]

t_list = np.linspace(t0, t1, 201)

a = 0.5; b = 0.2
answer1 = solve_ivp(Fun, t_span, y_ini, 
                    t_eval = t_list, 
                    args=(a, b), 
                    rtol = 1.e-12, 
                    atol = 1.e-14)

# 数値解は .y[0] のようにして取り出す
plt.plot(answer1.t, answer1.y[0], label="a=0.5, b=0.2");
plt.grid();
plt.xlabel("t");
plt.ylabel("x(t)");
plt.legend();
plt.title("減衰＋強制振動");
plt.xlim(0, 20);