scipyのodeintを理解したかった。

概要

scipyのodeintを理解したかった。

写真

サンプルコード

import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def F(d, t):
	rho, theta, z = d
	drhodt = 0
	dthetadt = 1
	dzdt = -1
	return [drhodt, dthetadt, dzdt]

rho0 = 1
theta0 = 0
z0 = 100
tspan = np.linspace(0, 20, 100)
sol = odeint(F, [rho0, theta0, z0], tspan)
rho = sol[ : , 0]
theta = sol[ : , 1]
z = sol[ : , 2]
X = rho * np.cos(theta)
Y = rho * np.sin(theta)
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.plot(X, Y, z)
plt.savefig('./odeint0.png')
plt.show()

写真

サンプルコード

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def F(v, t, p, r, b):
	return [-p * v[0] + p * v[1], -v[0] * v[2] + r * v[0] - v[1], v[0] * v[1] - b * v[2]]

v0 = [0.3, 0.3, 0.3]
t = np.arange(0, 100, 0.01)
p = 10
r = 28
b = 3 / 3
v = odeint(F, v0, t, args = (p, r, b))
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.plot(v[ : , 0], v[ : , 1], v[ : , 2])
plt.savefig('./odeint1.png')
plt.show()

写真

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint

def F(X, t):
	x, y, z = X
	return [-y - z, x + a * y, b + z * (x - c)]

dt = 5e-3
x = 0
y = 0
z = 0
x0 = (1, 1, 1)
a = 0.1
b = 0.1
c = 10
t_werte = np.arange(0, 500, dt)
erg = odeint(F, x0, t_werte)
x_, y_, z_ = erg[len(erg) / 4 : ].transpose()
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.plot(x_, y_, z_, "k,", alpha = 0.2, rasterized = True)
plt.savefig('./odeint10.png')
plt.show()

以上。