scipyのsolve_ivpでローレンツ方程式を解く

以下のページを参考にした． https://org-technology.com/posts/ordinary-differential-equations.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def lorenz(t, state, p, r, b):
    x, y, z = state
    return [-p * x + p * y, -x * z + r * x - y, x * y - b * z]
    pass


p = 10
r = 28
b = 8 / 3

v0 = [0.1, 0.1, 0.1]
t_end = 100
t = np.arange(0, t_end, 0.01)

sol = solve_ivp(lorenz, [0, t_end], v0, args=(p, r, b), dense_output=True)

z = sol.sol(t)

plt.plot(t, z.T)
plt.show()