用 Python 逆向求解 ODE答案

【问题标题】：Solving ODE with Python reversely用 Python 逆向求解 ODE
【发布时间】：2019-06-07 06:03:40
【问题描述】：

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

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.59和y0 = 0.69的答案，我怎么能得到z？它应该是0.153，但考虑到我们不知道它的价值，我们无法通过反复试验来找到它。

感谢您的帮助。

【问题讨论】：

您的目标是导数函数还是解的特定值？有了长长的介绍，我不确定它是否真的是第一个。这两种情况的答案都是“智能试错法”，也称为数值非线性求解器。可以使用 scipy.optimize.fsolve (?)。在第二种情况下也是边界值求解器或solve_ivp的事件机制。
@lutzl 我想获得z。我们有dh/dz，所以我们选择一个z 来查找该点的值。现在我有了答案，但我不知道在哪个z。
fsolve 需要初步猜测，我设置了z0=0 是的，你需要减去目标值z=scipy.optimze.fsolve(lambda vz: OD_H(y0,vz,c,b)-target, 0)。使用 odeint，您可以计算解决方案并检查它通过给定值的位置，并通过插值进行优化。一种粗略的方法是将 fsolve 应用于lambda vz: odeint(OD_H, y0, [0,vz])[-1,0]-target。这是昂贵的，因为它集成了每个函数调用，而使用solve_ivp的事件机制只需要一个集成。
是的，target 是目标值，0.59 是问题中唯一给出的值。仍然不完全清楚您到底想要计算什么。您将如何测试求解器返回的值实际上是一个解决方案？鉴于正向问题，更容易确定如何解决正确的逆向问题。
这是变量z的lambda-local值。你不应该需要它，因为它的作用域是 lambda 表达式，它也被声明了。

标签： python ode

【解决方案1】：

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

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])

【讨论】：

你不知道这个。只是我知道我们有初始点，我们有 ODE 的答案，但我们没有 z，即获得结果的点。
没错，x 就是那个点，作为您需要使用 ODE 求解器创建的可调用函数的变量。在抽象形式中，f(z) = your_ode_solver(args, z) - 0.59。然后，求根方法在 z 上迭代，直到 f(z) 接近于零。发生这种情况时，求根方法返回的 z 值是 ODE 求值的解函数为 0.59 的点。
如果您的答案正是我想要的，我可以请您编写代码吗？非常感谢。
exponential_decay = lambda t, y: -0.5 * y # dy/dt = f(t, y) 可以考虑这个吗？ exponential_decay = lambda t, y: (-0.5 * y)/(1+t)
问题是因为 odr[i] 小于解的最小值 screenshots.firefox.com/18zQqsfdGSvNmLkF/localhost 。 ODE 解中传播的 odr 错误或数值错误。