0

我正在使用 scipy.integrate 来解决 ODE 系统。我的代码摘要如下所示：

# step 1: import all necessary modules and functions imports here
import numpy as np
from scipy.integrate import ode

# step 2: define ode sys  
def ode_rhs(t, y, p1, p2, ..., pn):   
  ''' do some stuff to define rhs''' 
  return ydot

# step 3: set parameters p1, p2, ..., pn 
"""
p1 = ...
p2 = ...
...
pn = ...
"""

# step 4: set time parameters
ts = 0.0
tf = 340.

# step 5: set initial conditions
y0 = np.zeros(some_shape) 

# step 6: set ode solver
#backend = 'dopri5'
backend = 'vode'
r = ode(ode_rhs).set_integrator(backend, method='BDF', atol=1e-12, rtol=1e-6)
r.set_initial_value(y0,ts)

# step 7: store solution
t   = [ts]
y  =  [y0]

# step 8: call ode rhs and integrate
while r.successful() and r.t < tf:
  r.set_f_params(p1, p2, ..., pn)
  r.integrate(tf, step=True)

  t.append(r.t)  
  y.append(r.y)  

# step 9: process output for plotting
t   = array(t)
y = array(y).T

求解器在用于推进求解的任何时间步长输出 y 值。我需要的是特定时间（0、1、2、... 339）的输出，而不是求解器内部使用的时间步长。我本可以使用 r.integrate(t+dt)。但是，我必须使用可变步进来减少仿真时间。

4

0 回答 0