An answer to this question on Stack Overflow.

Question

I have been given two second order ODEs and I've been asked to solve them with odeint in python.

These are the equations:

d^x(t)/dt^2 = 10dy(t)/dt + x(t) - (k + 1)(x(t))/z^3
d^2y(t)/dt^2 = - 10dy(t)/dt + y(t) - ((k+1)(y(t) + k))/z^3

where z = np.sqrt((y+k)^2+x^2))

I've been given initial variables (x, y, dxdt, dydt) I know their values but I'm not stuck on typing them so I won't put them here.

def function(init, time, k):
    xt = init[0]
    yt = init[1]
    z = np.sqrt((init[1]+k)^2+init[0]^2))
    dxdt = init[2]
    dydt = init[3]
    ddxddt = 10*dydt + xt - ((k+1)(xt))/z^3
    ddyddt = -10*dxdt + xt - ((k+1)(yt + k))/z^3
    return(dxdt, ddxddt, dydt, ddyddt)
init = [0.921, 0, 0, 3.0]
values = odeint(function, initial, time, args(k,))

After this I define initial, and define time, k, and put them in the odeint.

But I can see that I am doing something really wrong with my actual set up function. I don't understand how to split up 2nd order odes.

Answer

You have a few mistakes here.

First: z^3 is not a power, it's the exclusive-or operation. In Python powers are done using the ** operator, so you'll want to be writing z**3.

Second: You've misnamed the arguments of your function. Instead of:

def function(init, time, k):

you should have

def function(state, time, k):

since state evolves according to the derivatives the function returns. It will only have the initial values in the first timestep.

Third: your state interpretation and state deltas are inconsistent. You write:

xt   = init[0]
yt   = init[1]
dxdt = init[2]
dydt = init[3]

But later

return dxdt, ddxddt, dydt, ddyddt

This implies, among other things, that dydt=ddxddt. You should instead write:

xt, yt, dxdt, dydt = state
[....]
return dxdt, dydt, ddxddt, ddyddt

Note that you then have to ensure your initial conditions agree with the way you've ordered your state.

A minimum working example of a correct implementation might look like this:

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
def function(state, time, k):
  xt,yt,dxdt,dydt = state
  z               = np.sqrt((yt+k)**2+xt**2)
  ddxddt          = 10*dxdt + xt - ((k+1)*(xt    ))/z**3
  ddyddt          = -10*dydt + yt - ((k+1)*(yt + k))/z**3
  return dxdt, dydt, ddxddt, ddyddt
init = [
  0.921, #x[0]
  0,     #y[0]
  0,     #x'[0]
  3.0    #y'[0]
]
k = 1
times  = np.linspace(0,1,1000)
values = scipy.integrate.odeint(function, init, times, args=(k,), tfirst=False)
plt.plot(values)
plt.show()

and gives this output: