代码之家 › 专栏 › 技术社区 › Migui Mag

在Python中模拟MATLAB中的ode45函数

ode45 differential-equations scipy matlab python

6

Migui Mag · 技术社区 · 7 年前

我想知道如何将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不同。

3 回复 | 直到 7 年前

1

9

Paul Wintz 5 年前

正如@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\u ivp解决方案 ,输出为 list(length of state vector) 一维Ndaray(长度等于 t\u评估 ).

因此,如果您想要相同的顺序,就必须将其合并。您可以通过以下方式执行此操作:

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

首先,您需要重塑,使其成为 [n,1] 数组,并将其水平合并。希望这有帮助!

2

4

user6655984 user6655984 7 年前

的接口 integrate.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()

3

2

Lost Fool David Contreras 3 年前

功能 西皮。整合solve\u ivp解决方案 使用方法 RK45 默认情况下,类似于Matlab函数使用的方法 ODE45 因为两者都使用四阶精度的休眠皮尔斯公式。

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

Ordinary Differential Equations (solve_ivp)