Numerical ODE solving in Python

Question

How do I numerically solve an ODE in Python?

Consider

\ddot{u}(\phi) = -u + \sqrt{u}

with the following conditions

u(0) = 1.49907

and

\dot{u}(0) = 0

with the constraint

0 <= \phi <= 7\pi.

Then finally, I want to produce a parametric plot where the x and y coordinates are generated as a function of u.

The problem is, I need to run odeint twice since this is a second order differential equation. I tried having it run again after the first time but it comes back with a Jacobian error. There must be a way to run it twice all at once.

Here is the error:

odepack.error: The function and its Jacobian must be callable functions

which the code below generates. The line in question is the sol = odeint.

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


def f(u, t):
    return -u + np.sqrt(u)


times = linspace(0.0001, 7 * np.pi, 1000)
y0 = 1.49907
yprime0 = 0
yvals = odeint(f, yprime0, times)

sol = odeint(yvals, y0, times)

x = 1 / sol * np.cos(times)
y = 1 / sol * np.sin(times)

plot(x,y)

plt.show()

Edit

I am trying to construct the plot on page 9

Classical Mechanics Taylor

Here is the plot with Mathematica

In[27]:= sol = 
 NDSolve[{y''[t] == -y[t] + Sqrt[y[t]], y[0] == 1/.66707928, 
   y'[0] == 0}, y, {t, 0, 10*\[Pi]}];

In[28]:= ysol = y[t] /. sol[[1]];

In[30]:= ParametricPlot[{1/ysol*Cos[t], 1/ysol*Sin[t]}, {t, 0, 
  7 \[Pi]}, PlotRange -> {{-2, 2}, {-2.5, 2.5}}]

Answer 1

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

pi = np.pi
sqrt = np.sqrt
cos = np.cos
sin = np.sin

def deriv_z(z, phi):
    u, udot = z
    return [udot, -u + sqrt(u)]

phi = np.linspace(0, 7.0*pi, 2000)
zinit = [1.49907, 0]
z = integrate.odeint(deriv_z, zinit, phi)
u, udot = z.T
# plt.plot(phi, u)
fig, ax = plt.subplots()
ax.plot(1/u*cos(phi), 1/u*sin(phi))
ax.set_aspect('equal')
plt.grid(True)
plt.show()

Answer 2

The code from your other question is really close to what you want. Two changes are needed:

• You were solving a different ODE (because you changed two signs inside function deriv)
• The y component of your desired plot comes from the solution values, not from the values of the first derivative of the solution, so you need to replace u[:,0] (function values) for u[:, 1] (derivatives).

This is the end result:

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

def deriv(u, t):
    return np.array([u[1], -u[0] + np.sqrt(u[0])])

time = np.arange(0.01, 7 * np.pi, 0.0001)
uinit = np.array([1.49907, 0])
u = odeint(deriv, uinit, time)

x = 1 / u[:, 0] * np.cos(time)
y = 1 / u[:, 0] * np.sin(time)

plt.plot(x, y)
plt.show()

However, I suggest that you use the code from unutbu's answer because it's self documenting (u, udot = z) and uses np.linspace instead of np.arange. Then, run this to get your desired figure:

x = 1 / u * np.cos(phi)
y = 1 / u * np.sin(phi)
plt.plot(x, y)
plt.show()