如何在指定时间点获得 ODE 解?

计算科学 算法 Python 数值分析
2021-12-18 09:56:05

下面的代码基本上说明了我的问题。这是一个钟摆的测试代码。我使用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()
3个回答

另请参阅:Python:带有步进控制 ODE 求解器的网格

GertVdE 的回答很好地涵盖了面向用户的实现。

本质上,您希望发生两种不同的事情:

  • 积分器应采取尽可能大的时间步长
  • 无论如何,您都希望在固定间隔的时间点输出

以下是它在实践中的工作方式:

  • 积分器有自己的内部时间步长信息。如果您要在“一步”模式下使用 VODE(在 SUNDIALS 中调用 Fortran 库 DVODE,或其后继者 CVODE),则需要尽可能大的一步(受稳定性和错误的限制)。然后它将输出该步骤。如果您要重复执行此操作,您的输出将是积分器的时间步长。但是,该输出不是您想要的——您想要定期间隔的输出。
  • 积分器还可以使用由积分器构造的多项式近似在内部时间步长之间进行插值。在 Hairer 和 Wanner 中,此功能称为“密集输出”。因此,当您在内部时间步之间查询积分器的输出(因为它及时生成)时,它真正做的是在这些时间步之间进行插值以获得输出。
  • GertVdE 描述的接口似乎使用“正常”模式,它自动进行密集输出并在内部采取时间步长,因为大多数用户关心他们得到的输出,而不关心实现中的数字细节,只要他们正常工作。

您可以在SUNDIALS 使用说明和各种SUNDIALS 用户指南中查看所有这些详细信息,其中涵盖了积分器背后的理论及其使用。尽管这些文档来源涵盖了 VODE 的后续版本(CVODE 和 CVODES),但理论应该相同,接口也应该相似,只是 CVODE 和 CVODES 具有附加功能。

scipy.integrate.ode 文档中给出的示例提示在您的兴趣点之间使用连续集成:

>>> r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
>>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
>>> t1 = 10
>>> dt = 1
>>> while r.successful() and r.t < t1:
>>>     r.integrate(r.t+dt)
>>>     print("%g %g" % (r.t, r.y))

scipy.integrate.odeint 函数采用输入参数“t”:求解 y 的时间点序列。初始值点应该是这个序列的第一个元素。

自 SciPy 1.0.0(2017 年末)以来的新功能:

另一种解决方案是解决solve_ivp并使用dense_output允许在解决方案步骤之间进行插值的选项:

import numpy as np
from scipy.integrate import ode, solve_ivp
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 = 100

t = np.linspace(t0,t1,2000)
sol = solve_ivp(f, (t0, t1), y0, dense_output=True)
y_real = sol.y
y_interp = sol.sol(t)

plt.clf()
plt.plot(y_real[0,:],y_real[1,:],'or')
plt.plot(y_interp[0,:],y_interp[1,:],'.-b')
plt.show()