python - solve_ivp differential equation solver, way to not integrate all return values?

Question

Hey i have a model which gives diferential return values. But also a value which isnt a differential (z)

def model(t, f):
x = f[0]
y = f[1]
dx_dt = 2x*y
dy_dt = x**3
z = 2*x

return dx_dt, dy_dt, z

My solver can solve this equations and give me x and y at the respective time values.

t = np.linspace(0, 10, 100)
f0 = [2, 1, 0]
result = solve_ivp(model, [np.min(t), np.max(t)], f0, t_eval=t)

But now i want also my solution for z which should NOT be integrated by the solver. Is there any possibility to just integrate the first 2 values of my solver? But not the third?

question from:https://stackoverflow.com/questions/65902080/solve-ivp-differential-equation-solver-way-to-not-integrate-all-return-values

Answer 1

You have the situation where the ODE function could be written as

def model(t,u):
    v = f(t,u)
    du = g(t,u,v) # v may contain components of du that get simply copied
    return du

and you are additionally interested in a quantity z=h(t,u,v). All variables can also be tuples or vectors. Presently the functions stand for code blocks, but in many cases they can also be easily separated out as functions. Thus the first variant would be just to do that, so that the ODE function has minimal extra functionality, and the values z are computed at the end from the solution values.

Another variant, if transforming the code blocks to separate functions appears not optimal, you can also construct the model function as in the question,

def model(t,u):
    v = f(t,u)
    du = g(t,u,v) # v may contain components of du that get simply copied
    z = h(t,u,v)
    return du,z

Then in the call of the solver you use a lambda expression to separate out the derivatives vector

result = solve_ivp(lambda t,u: model(t,u)[0], t[[0,-1]], f0, t_eval=t)

and from the result by calling the same model function again you get

z = model(t,result.y.T)[1]

or use an iterator formulation if the model function can not be automatically vectorized

z = [ model(t_,u_)[1] for t_, u_ in zip(t,result.y.T) ]
z = np.array(z).T # to get the same index ordering as in result.y

Use appropriate array slicing operations if the model function returns just one list of values and not a pair of tuples.