python - 使用“dopri5”scipy.integrate.ode 求解器求解 ODE 系统时出现错误消息

Question

我正在尝试解决一个 ODE 系统，该系统将在势能下在 2D 空间中移动的粒子的位置和速度描述为时间的函数，使用dopri5in scipy.integrate.ode：

我查找了解决这种 ODE 的示例，然后我基本上复制了以下链接中的所有内容，只需切换变量和维度即可解决我的特定问题：

如何在 Python 中使用 dorpi5 或 dop853

这是代码：

from scipy.integrate import ode
from scipy.interpolate import RegularGridInterpolator
import numpy as np
import matplotlib.pyplot as plt

# Define V(x,y) and obtain the gradient function through interpolation:

def V(x, y):
    return np.exp((x**2+y**2))-0.5

x = np.linspace(-10, 10, 1000)
y = np.linspace(-10, 10, 1000)
X, Y = np.meshgrid(x, y)
v = V(X, Y)
gradv_array = np.gradient(v, x, y)

gradv_x = RegularGridInterpolator((x, y), gradv_array[0])
gradv_y = RegularGridInterpolator((x, y), gradv_array[1])


def fun(t, z, m):
    """
    Right hand side of the differential equations
      dx/dt = v_x
      dy/dt = v_y
      dv_x/dt = -gradv_x(x, y)/m 
      dv_y/dt = -gradv_y(x, y)/m
    """
    x, y, v_x, v_y = z
    f = [v_x, v_y, -gradv_x(x, y)/m, -gradv_y(x, y)/m]
    return f

# Create an `ode` instance to solve the system of differential
# equations defined by `fun`, and set the solver method to 'dopri5'.

solver = ode(fun)
solver.set_integrator('dopri5')

# Give the value of m to the solver. This is passed to
# `fun` when the solver calls it.
m = 1
solver.set_f_params(m)

# Set the initial value z(0) = z0.
t0 = 0.0
z0 = [0.5, 0.3, 0, 0]
solver.set_initial_value(z0, t0)

# Create the array `t` of time values at which to compute
# the solution, and create an array to hold the solution.
# Put the initial value in the solution array.
t1 = 2.5
N = 75
t = np.linspace(t0, t1, N)
sol = np.empty((N, 4))
sol[0] = z0

# Repeatedly call the `integrate` method to advance the
# solution to time t[k], and save the solution in sol[k].
k = 1
while solver.successful() and solver.t < t1:
    solver.integrate(t[k])
    sol[k] = solver.y
    k += 1

# Plot the solution...
plt.plot(t, sol[:,0], label='v_x')
plt.plot(t, sol[:,1], label='v_y')
plt.xlabel('t')
plt.grid(True)
plt.legend()
plt.show()

plt.plot(t, sol[:,2], label='x')
plt.plot(t, sol[:,3], label='y')
plt.xlabel('t')
plt.grid(True)
plt.legend()
plt.show()

我尝试执行原始代码和修改后的代码，前者有效，后者无效，并显示以下错误消息：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-39-2db0a4c50fa3> in <module>()
     42 k = 1
     43 while solver.successful() and solver.t < t1:
---> 44     solver.integrate(t[k])
     45     sol[k] = solver.y
     46     k += 1

~\Anaconda3\lib\site-packages\scipy\integrate\_ode.py in integrate(self, t, step, relax)
    430             self._y, self.t = mth(self.f, self.jac or (lambda: None),
    431                                   self._y, self.t, t,
--> 432                                   self.f_params, self.jac_params)
    433         except SystemError:
    434             # f2py issue with tuple returns, see ticket 1187.

~\Anaconda3\lib\site-packages\scipy\integrate\_ode.py in run(self, f, jac, y0, t0, t1, f_params, jac_params)
   1170     def run(self, f, jac, y0, t0, t1, f_params, jac_params):
   1171         x, y, iwork, istate = self.runner(*((f, t0, y0, t1) +
-> 1172                                           tuple(self.call_args) + (f_params,)))
   1173         self.istate = istate
   1174         if istate < 0:

<ipython-input-39-2db0a4c50fa3> in fun(t, z, m)
     10     """
     11     x, y, v_x, v_y = z
---> 12     f = [v_x, v_y, -gradv_x(x, y)/m, -gradv_y(x, y)/m]
     13     return f
     14 

~\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py in __call__(self, xi, method)
   2460         method = self.method if method is None else method
   2461         if method not in ["linear", "nearest"]:
-> 2462             raise ValueError("Method '%s' is not defined" % method)
   2463 
   2464         ndim = len(self.grid)

ValueError: Method '0.3' is not defined

显然yinsolver.y被解释为该列表中的第二个元素：

z0 = [0.5, 0.3, 0, 0]

也许应该做出改变，因为有四个 ODE 而不是两个？但是怎么做？

编辑：发现错误。生成的函数scipy.interpolate.RegularGridInterpolator接受参数列表，但不接受多个参数，即gradv_x([x, y])代替grad_x(x, y).

标签： pythonnumpyscipyodesolver