4

Is it possible to integrate any Ordinary Differential Equation backward in time using scipy.integrate.odeint ? If it is possible, could someone tell me what should be the arguement 'time' in 'odeint.

ADK
  • 259
  • 2
  • 17
  • 1
    I suspect you need to consult a mathematician. Mathematically speaking, backwards integration of many classes of ODEs is *ill-posed*, ie there is no guarantee that a solution exists – talonmies Nov 05 '12 at 07:49

3 Answers3

4

odeint handles negative values of the t argument. No special treatment is needed.

Here's an example:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def mysys(z, t):
    """A slightly damped oscillator."""
    return [z[1] - 0.02*z[0], -z[0]]


if __name__ == "__main__":
    # Note that t starts at 0 and goes "backwards"
    t = np.linspace(0, -50, 501)

    z0 = [1, 1]
    sol = odeint(mysys, z0, t)

    plt.plot(t, sol)
    plt.xlabel('t')
    plt.show()

The plot: Solve an ODE backward in time

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • I was considering writing a code of a forward/reverse shooting method to find a stable manifold when an equilibrium is saddle unstable. But, this answer has solved all problems I had. Just adding a minus sign in code enables me to find a stable manifold! You made my day! – T_T Nov 04 '18 at 09:50
0

You can make a change of variables s = t_0 - t, and integrate the differential equation with respect to s. odeint doesn't do this for you.

pv.
  • 33,875
  • 8
  • 55
  • 49
  • Thanks a lot buddy. However, how would you modify your above code to introduce time delay (say of time_delay=1.5) in slightly damped oscillator ? This is what I exactly want to do with my system. – ADK Nov 06 '12 at 09:01
0

It is not necessary to make a change of variables. Here an example:

import math
import numpy
import scipy
import pylab as p

from math import *
from numpy import *
from scipy.integrate import odeint
from scipy.interpolate import splrep
from scipy.interpolate import splev

g1=0.01
g2=0.01
w1=1
w2=1 
b1=1.0/20.0
b2=1.0/20.0
b=1.0/20.0
c0=0
c1=0.2
wf=1

def wtime(t):
  f=1+c0+c1*cos(2*wf*t)
  return f

def dv(y,t):
  return array([y[1], -(wtime(t)+g1*w1+g2*w2)*y[0]+w1*y[2]+w2*y[3], w1*y[2]-g1*w1*y[0], w2*y[3]-g2*w2*y[0]])

tv=linspace(100,0,1000)
v1zero=array([1,0,0,0])
v2zero=array([0,1,0,0])
v1s=odeint(dv,v1zero,tv)
v2s=odeint(dv,v2zero,tv)

p.plot(tv,v1s[:,0])
p.show()

I check the result with Wolfram Mathematica (that program can solve backward odes).

Andrés Felipe Estrada Guerra
  • 1
  • 1