用 Python 反向求解 ODE

Solving ODE with Python reversely

在Python中我们求解一个初始点y0 = od0在特定点z的微分方程OD_H类似于下面的代码

def OD_H(od, z, c, b):
   ....
   return ....

od = solve_ivp(lambda od, z: OD_H(od, z, c, b), t_span = [z1, z], y0 = [od0])['y'][-1][-1]

od = odeint(OD_H, od0, [0, z], args=(c, b))[-1]

所以我们有

answer of ODE OD_H(y0 = 0.69, z=0.153) = 0.59

我的问题是,

现在如果我有OD_H = 0.59y0 = 0.69的答案,我怎么能得到z呢?它应该是 0.153 但考虑到我们不知道它的值,我们无法通过反复试验来找到它。

感谢您的帮助。

在这种情况下,您提出了一个根问题,其中计算的求解器函数减去所需答案是函数 f(x),其中 f(x)=0。
因为 ODE 求解器 returns 点数组而不是可调用函数,所以您需要先对解点进行插值。然后,这用于求根问题。

from scipy.integrate import solve_ivp # Recommended initival value problem solver
from scipy.interpolate import interp1d # 1D interpolation
from scipy.optimize import brentq # Root finding method in an interval
exponential_decay = lambda t, y: -0.5 * y # dy/dt = f(t, y)
t_span = [0, 10] # Interval of integration
y0 = [2] # Initial state: y(t=t_span[0])=2
desired_answer = 0.59
sol_ode = solve_ivp(exponential_decay, t_span, y0) # IVP solution
f_sol_ode = interp1d(sol_ode.t, sol_ode.y) # Build interpolated function
brentq(lambda x: f_sol_ode(x) - desired_answer, t_span[0], t_span[1])