solve_ivp 返回不同的 odeint 结果？

Question

我正在尝试求解一个简单的 ODE，以便理解 Scipy 的新 API。

我为 4 阶的 Runge Kutta 编写了一个例程来编写它并用旧的 API odeint 确认它并且它匹配得很好。但现在我正试图绕过 solve_ivp，它似乎不起作用。我哪里错了？

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp, odeint
import time

freq = np.arange(1, 10000, 100)

def g(q, t):
    return -q ** 3 + np.sin(t)


a = 0
b = 10
npoints = 100
h = (b - a) / npoints
t = np.arange(a, b, h)

output1 = np.zeros(t.shape)
x = 0
for i in range(len(t)):
    output1[i] = x
    k1 = h * g(x, t[i])
    k2 = h * g(x + 0.5 * k1, t[i] + 0.5 * h)
    k3 = h * g(x + 0.5 * k2, t[i] + 0.5 * h)
    k4 = h * g(x + k3, t[i] + 0.5 * h)
    x = x + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)


# ---------------Solving using odeint (old API)---------------#

y1_odeint = odeint(g, 0, t)

#---------------Solving using new API-------------#

y2=solve_ivp(g,(a,b),[0],t_eval=t)


# --------------------Representação gráfica--------------------------#
fig = plt.figure()
ax = fig.add_subplot(121)
ax1=fig.add_subplot(122)

ax.plot(t, output1,label="my own")
ax.plot(t,y1_odeint,label="odeint")
ax.plot(y2.t,np.squeeze(y2.y),label="new API")
ax.legend()
ax.set_title("Output")

ax1.plot(t,output1-np.squeeze(y1_odeint),label="|odeint-my own|")
ax1.legend()


plt.tight_layout()
plt.show()

Answer 1

再看一下文档字符串，因为 solve_ivp. It expects the first argument of g to be t. By default, odeint 使用了相反的约定。如果你有最新版本的 scipy，你可以告诉 odeint 第一个参数是 t，方法是给它参数 tfirst=True.

solve_ivp 返回不同的 odeint 结果？

solve_ivp returning different outcome of odeint?

python

scipy

ode