CODEWITHSUNDEEP
pythonnumpyscipyodeodeint
Hersh Bhargava
Hersh Bhargava

Reputation: 615

Override values in ODE integration with scipy.integrate.odeint

I am conducting numerical integration of a system of ODEs using scipy.integrate.odeint. In the abstract, I want to be able to give the solver an arbitrary condition to apply to its output values. For instance: if the value of a term (e.g. X[0]) is less than 1, override the value to 0, rather than using the value computed by numerical integration. Then carry on integration.

Here is an example of how I approach the integration.

# Toy ODE system to be integrated
def dXdt(X, t, args):
    return [args[0]*X[0], t*args[1]*X[1]]

# Setup initial conditions and integrate the system.
from scipy.integrate import odeint
import numpy as np
t = np.linspace(0, 100, 100) # timepoints to integrate
x0 = [100, 100] # initial conditions
args = [-.01, -.02]
X = odeint(dXdt, x0, t, args=(args,))

A workaround to achieve approximately the behavior I want would be to modify dXdt as follows:

def dXdt(X, t, args):
    if X[0] > 1: # make dX[0]/dt extremely negative to force X[0] toward 0
        return [-1e200, t*args[1]*X[1]] 
    return [args[0]*X[0], t*args[1]*X[1]]

This is obviously not generalizable to other conditions that might be desirable, and also does not solve the problem precisely (a precise solution would require knowledge of the solver's step size in dXdt.

Upvotes: 1

Views: 212

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

This is not possible. You have to consider odeint as black box where you have no access to the internal state and the steps taken during the integration. In other contexts (languages, libraries) you have an event-action mechanism, where you can define the event of crossing some value downwards and take the action to modify the state.

In scipy.integrate.solve_ivp you have one-half of that implemented, the event mechanism. You would have to define the event as terminal, extract the last state at the event, modify it and restart the integration with the modified state. In the end you can then concatenate the solution segments.

Upvotes: 2

Related Questions