6

我想知道如何将 MATLAB 函数 ode45 导出到 python。根据文档应该如下：

 MATLAB:  [t,y]=ode45(@vdp1,[0 20],[2 0]);

 Python:  import numpy as np
          def  vdp1(t,y):
              dydt= np.array([y[1], (1-y[0]**2)*y[1]-y[0]])
              return dydt
          import scipy integrate 
          l=scipy.integrate.ode(vdp1([0,20],[2,0])).set_integrator("dopri5")

结果完全不同，Matlab 返回的维度与 Python 不同。

4

3 回答 3

9

正如@LutzL 提到的，您可以使用较新的 API solve_ivp，.

results = solve_ivp(obj_func, t_span, y0, t_eval = time_series)

如果t_eval未指定，则每个时间戳不会有一条记录，这主要是我假设的情况。

另一个注意事项是，对于odeint其他积分器来说，输出数组ndarray的形状通常是[len(time), len(states)]，但是对于solve_ivp，输出是一list(length of state vector)维 ndarray （长度等于t_eval）。

因此，如果您想要相同的顺序，则必须合并它。你可以这样做：

Y =results
merged = np.hstack([i.reshape(-1,1) for i in Y.y])

首先，您需要对其进行整形以使其成为一个[n,1]数组，然后将其水平合并。希望这可以帮助！

于 2018-04-19T15:15:33.883 回答
4

Integrated.ode的界面不像更简单的方法odeint那样直观，但是它不支持选择 ODE 积分器。主要区别在于它ode不会为您运行循环；如果你需要一个解决方案在一堆点上，你必须说在什么点上，并一次计算一个点。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def vdp1(t, y):
    return np.array([y[1], (1 - y[0]**2)*y[1] - y[0]])
t0, t1 = 0, 20                # start and end
t = np.linspace(t0, t1, 100)  # the points of evaluation of solution
y0 = [2, 0]                   # initial value
y = np.zeros((len(t), len(y0)))   # array for solution
y[0, :] = y0
r = integrate.ode(vdp1).set_integrator("dopri5")  # choice of method
r.set_initial_value(y0, t0)   # initial values
for i in range(1, t.size):
   y[i, :] = r.integrate(t[i]) # get one more value, add it to the array
   if not r.successful():
       raise RuntimeError("Could not integrate")
plt.plot(t, y)
plt.show()

解决方案

于 2018-01-25T18:01:19.650 回答
2

函数scipy.integrate.solve_ivp默认使用 RK45 方法，类似于 Matlab 的函数 ODE45 使用的方法，因为两者都使用具有四阶方法精度的 Dormand-Pierce 公式。

vdp1 = @(T,Y) [Y(2); (1 - Y(1)^2) * Y(2) - Y(1)];
[T,Y] = ode45 (vdp1, [0, 20], [2, 0]);

from scipy.integrate import solve_ivp

vdp1 = lambda T,Y: [Y[1], (1 - Y[0]**2) * Y[1] - Y[0]]
sol = solve_ivp (vdp1, [0, 20], [2, 0])

T = sol.t
Y = sol.y

常微分方程 (solve_ivp)

于 2021-01-27T23:56:47.777 回答