5

I have a second order differential equation that I want to solve it in python. The problem is that for one of the variables I don't have the initial condition in 0 but only the value at infinity. Can one tell me what parameters I should provide for scipy.integrate.odeint ? Can it be solved?

Equation: enter image description here

Theta needs to be found in terms of time. Its first derivative is equal to zero at t=0. theta is not known at t=0 but it goes to zero at sufficiently large time. all the rest is known. As an approximate I can be set to zero, thus removing the second order derivative which should make the problem easier.

rowman
  • 1,516
  • 1
  • 16
  • 26
  • You can find the answer from this question http://stackoverflow.com/questions/1824751/infinity-generated-in-python-code – Mihai8 Jan 20 '13 at 16:28
  • You cannot give `scipy.integrate.odeint` conditions at two different values of `t`. If instead of a second condition at infinity, you had it at `t = t1`, you could nest your solution with `scipy.integrate.odeint` inside a call to `scipy.optimize.root` to find the value of tetha at `t = 0` that gave you the desired behavior at `t = t1`. Maybe choosing a large enough `t1` would allow you to use that idea. You may also want to try http://scicomp.stackexchange.com/ for help figuring the right strategy to tackle your problem. – Jaime Jan 21 '13 at 08:09
  • @Jaime, Could you please provide a rough answer by using `scipy.optimize.root` and predicting `t1` value? – rowman Jan 21 '13 at 20:35

1 Answers1

4

This is far from being a full answer, but is posted here on the OP's request.

The method I described in the comment is what is known as a shooting method, that allows converting a boundary value problem into an initial value problem. For convenience, I am going to rename your function theta as y. To solve your equation numerically, you would first turn it into a first order system, using two auxiliary function, z1 = y and z2 = y', and so your current equation

I y'' + g y' + k y = f(y, t)

would be rewitten as the system

z1' = z2
z2' = f(z1, t) - g z2 - k z1

and your boundary conditions are

z1(inf) = 0
z2(0) = 0

So first we set up the function to compute the derivative of your new vectorial function:

def deriv(z, t) :
    return np.array([z[1],
                     f(z[0], t) - g * z[1] - k * z[0]])

If we had a condition z1[0] = a we could solve this numerically between t = 0 and t = 1000, and get the value of y at the last time as something like

def y_at_inf(a) :
    return scipy.integrate.odeint(deriv, np.array([a, 0]),
                                  np.linspace(0, 1000, 10000))[0][-1, 0]

So now all we need to know is what value of a makes y = 0 at t = 1000, our poor man's infinity, with

a = scipy.optimize.root(y_at_inf, [1])
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • I haven't read in detail, but what I typically do with equations like that is solve them backward from the stationary state back to the past, running time backwards, i.e. ``t --> - \infty`` – PatrickT Jun 18 '16 at 11:55