Solve_ivp Differential Equation Solver, Way To Not Integrate All Return Values?

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 retur

Solution 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 copiedreturndu

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)
    returndu,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.

Solution 2:

Hey thank you for your solution. i think this solution may also work:

defmodel(t,y):

global inner_func
definner_func(t,y):
    
    #differential Eq Systemreturn dx_dt, dy_dt, z

dx_dt, dy_dt, z = inner_func(t,y)
return dx_dt, dy_dt

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

z = np.array([inner_func(t[i], result[i,:]) for i inrange(len(t))])[:,2]