Solve an implicit ODE (differential algebraic equation DAE)

Question

I'm trying to solve a second order ODE using odeint from scipy. The issue I'm having is the function is implicitly coupled to the second order term, as seen in the simplified snippet (please ignore the pretend physics of the example):

import numpy as np
from scipy.integrate import odeint

def integral(y,t,F_l,mass):
    dydt = np.zeros_like(y)
    x, v = y
    F_r =  (((1-a)/3)**2 + (2*(1+a)/3)**2) * v # 'a' implicit 
    a  = (F_l - F_r)/mass

    dydt = [v, a]
return dydt


y0 = [0,5]
time = np.linspace(0.,10.,21)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon,mass))

in this case I realise it is possible to algebraically solve for the implicit variable, however in my actual scenario there is a lot of logic between F_r and the evaluation of a and algebraic manipulation fails.

I believe the DAE could be solved using MATLAB's ode15i function, but I'm trying to avoid that scenario if at all possible.

My question is - is there a way to solve implicit ODE functions (DAE) in python( scipy preferably)? And is there a better way to pose the problem above to do so?

As a last resort, it may be acceptable to pass a from the previous time-step. How could I pass dydt[1] back into the function after each time-step?

Answer 1

if algebraic manipulation fails, you can go for a numerical solution of your constraint, running for example fsolve at each timestep:

import sys
from numpy import linspace
from scipy.integrate import odeint
from scipy.optimize import fsolve

y0 = [0, 5]
time = linspace(0., 10., 1000)
F_lon = 10.
mass = 1000.

def F_r(a, v):
    return (((1 - a) / 3) ** 2 + (2 * (1 + a) / 3) ** 2) * v

def constraint(a, v):
    return (F_lon - F_r(a, v)) / mass - a

def integral(y, _):
    v = y[1]
    a, _, ier, mesg = fsolve(constraint, 0, args=[v, ], full_output=True)
    if ier != 1:
        print "I coudn't solve the algebraic constraint, error:\n\n", mesg
        sys.stdout.flush()
    return [v, a]

dydt = odeint(integral, y0, time)

Clearly this will slow down your time integration. Always check that fsolve finds a good solution, and flush the output so that you can realize it as it happens and stop the simulation.

About how to "cache" the value of a variable at a previous timestep, you can exploit the fact that default arguments are calculated only at the function definition,

from numpy import linspace
from scipy.integrate import odeint

#you can choose a better guess using fsolve instead of 0
def integral(y, _, F_l, M, cache=[0]):
    v, preva = y[1], cache[0]
    #use value for 'a' from the previous timestep
    F_r = (((1 - preva) / 3) ** 2 + (2 * (1 + preva) / 3) ** 2) * v 
    #calculate the new value
    a = (F_l - F_r) / M
    cache[0] = a
    return [v, a]

y0 = [0, 5]
time = linspace(0., 10., 1000)
F_lon = 100.
mass = 1000.

dydt = odeint(integral, y0, time, args=(F_lon, mass))

Notice that in order for the trick to work the cache parameter must be mutable, and that's why I use a list. See this link if you are not familiar with how default arguments work.

Notice that the two codes DO NOT produce the same result, and you should be very careful using the value at the previous timestep, both for numerical stability and precision. The second is clearly much faster though.

Answer 2

Quite Old , but worth updating so it may be useful for anyone, who stumbles upon this question. There are quite few packages currently available in python that can solve implicit ODE. GEKKO (https://github.com/BYU-PRISM/GEKKO) is one of the packages, that specializes on dynamic optimization for mixed integer , non linear optimization problems, but can also be used as a general purpose DAE solver.

The above "pretend physics" problem can be solved in GEKKO as follows.

m= GEKKO()
m.time = np.linspace(0,100,101)
F_l = m.Param(value=1000)
mass = m.Param(value =1000)
m.options.IMODE=4
m.options.NODES=3
F_r = m.Var(value=0)
x = m.Var(value=0)
v = m.Var(value=0,lb=0)
a = m.Var(value=5,lb=0)
m.Equation(x.dt() == v)
m.Equation(v.dt() == a)
m.Equation (F_r ==  (((1-a)/3)**2 + (2*(1+a)/3)**2 * v)) 
m.Equation (a == (1000 - F_l)/mass)
m.solve(disp=False)
plt.plot(x)