Speed issues, python not responding: Solving ODEs

Question

f = odefun(lambda x,y: (x)/((x**2) + (y**2))**(3.0/2.0),0,1)
for x in range(500):
    listx.append(f(x))

For a quick background I'm trying to do a gravitation/orbital simulation. I create two lists listx and listy then I put them into matplot and plot the results (hopefully and ellipse arises).

The issue is that to solve this equation 500 times takes a noticeably long time plus, it doesn't give me an adequate amount of plotting points. So, if I step it up to 1000 it takes waaay too long, and in some cases I have to kill python due to it being unresponsive.

I'm just wondering if there's a more efficient way (I'm sure there is) of populating these lists.

Answer 1

Use the odeint function included in scipy.integrate. The following code snippets are equivalent. This is the version using mpmath (included as part of sympy)

from sympy.mpmath import odefun
from matplotlib import pyplot as plt
f = odefun(lambda x, y: x / (x**2 + y**2)**1.5, 0, 1)
X = range(1000)
listx = [f(x) for x in X]
plt.plot(X, listx)
plt.show()

This is the version using NumPy and SciPy

import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
T = np.arange(1000)
xarr = odeint(lambda t, x: x / (x**2 + t**2)**1.5, 1, T)
plt.plot(T, xarr)
plt.show()

Depending on which points you want to plot, it may be easier to use NumPy's linspace function that lets you take equispaced points between two values. As it is, you are just plotting Notice that I did have to swap the order of the arguments. Be aware that the second version uses NumPy arrays and not Python lists. If you want to learn more about NumPy and SciPy, good places to start are the scipy lecture notes and the SciPy wiki's numpy tutorial.

Depending on what you need, you could also take a look at the ode class that is included in scipy.integrate.

I would only recommend using mpmath when you need the arbitrary precision arithmetic it provides. If ordinary floating point arithmetic is good enough, mpmath will be much slower than many of your other options. If you still need to use mpmath, I would recommend installing gmpy to speed things up. That will help, but it will still be much faster to use SciPy.