Solving ODE with Python reversely

Question

In Python we solve a differential equation OD_H with an initial point y0 = od0 in a specific point z similar to the following code

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]

or

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

So we have

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

My question is,

Now If I have the answer of OD_H = 0.59 and y0 = 0.69, how could I obtain z? It should be 0.153 but consider we don't know its value and we cannot do trial and error to find it.

I appreciate your help.

Answer 1

In this case, you are proposing a root problem where the solver function evaluated minus the desired answer is the function f(x) where f(x)=0.
Because ODE solver returns point arrays and not callable functions, you need to interpolate first the solution points. Then, this is used in root finding problem.

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