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…

Accepted Answer

```
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()
```

Numerical ODE solving in Python

Solution

Problem

Related problems