下面的代码基本上说明了我的问题。这是一个钟摆的测试代码。我使用https://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode上建议的方法解决它。现在我想在一个均匀间隔的时间范围内绘制点,以显示钟摆摆动的动态特性。例如我想得到解决方案并绘制相应的解决方案状态空间点
[0,0.2,0.4,0.6,0.8,1.0]
Python 是否有任何自动方式可以做到这一点(使用在 ODE 集成期间获得的内部解决方案点的插值)?
这是主要代码:
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
import warnings
def f(t,y):
l = 1
m = 1
d = 1
g = 9.8
return [y[1], -np.sin(y[0])*g/l-y[1]*d/m]
y0, t0 = [np.pi/2, 0], 0
t1 = 500
backend = 'vode'
solver = ode(f).set_integrator(backend, nsteps=1)
solver.set_initial_value(y0, t0)
# suppress Fortran-printed warning
solver._integrator.iwork[2] = -1
y, t = [], []
warnings.filterwarnings("ignore", category=UserWarning)
while solver.t < t1:
solver.integrate(t1, step=True)
y.append(solver.y)
t.append(solver.t)
warnings.resetwarnings()
y = np.array(y)
t = np.array(t)
plt.plot(y[:,0], y[:,1], 'b-')
plt.plot(y[0::5,0], y[0::5,1], 'b.')
plt.show()