# Implementing Euler's method for initial value ODEs

In my physics class, I had to calculate the trajectory of a projectile that was fired (very fast) with $v_0$ in an angle off a planet (radius $R$, mass $M$) from the surface. The projectile would reach a maximum distance and start to circle around the planet. The question was whether it would hit the planet. We solved that analytically, but I wonder how to simulate such a problem numerically.

So basically I start off with

$$\vec{r}_0 = \begin{pmatrix} 0 \\ 0 \\ R \end{pmatrix}$$

and

$$\vec{v}_0 = v_0 \begin{pmatrix} \cos(30°) \\ 0 \\ \sin(30°) \end{pmatrix}$$

The acceleration on the projectile is just

$$\vec{a}(\vec{r}) = \frac{\vec{F}(\vec{r})}{m} = G \frac{M}{r^3} \vec{r}$$

I could implement a simple C++/Python program which basically does:

$$\vec{r}_{t+1} = \vec{r}_{t} + \vec{v}_{t} \cdot \mathrm dt$$ $$\vec{v}_{t+1} = \vec{v}_{t} + \vec{a}_{t} \cdot \mathrm dt$$

But this seems a little unsuited as in either language I do not have "natural" vector classes.

Would Octave or Mathematica be better suited for this task? How would I implement this problem to obtain a list of x, y, z coordinates in order to plot the trajectory in gnuplot?

• I asked this on Physics, and I was told to ask over here. If it is a basic problem, a hint to the right direction should be easy I assume. – Martin Ueding Jan 2 '12 at 13:34
• Sure, take a look at any book on numerical analysis or scientific computing and read the chapter on initial value ODEs. – David Ketcheson Jan 2 '12 at 15:13
• Hint: Instead of "natural vector classes", for n-dimensional vectors use n-dimensional arrays of scalars. Also, what David said. – milancurcic Jan 2 '12 at 15:34
• The problem in Python is that if r = [1, 0], r+r will be [1, 0, 1, 0], not [2, 0]. I tried it in NumPy now and it works kinda usable. – Martin Ueding Jan 2 '12 at 17:22
• @DavidKetcheson: Is there a concensus that hw questions are off-topic here? My impression was that level of difficulty was not considered a criterion for exclusion. – Faheem Mitha Jan 2 '12 at 17:48

Just use Python + NumPy. Here is the code

from math import pi, sin, cos
from numpy import array, vstack
from numpy.linalg import norm

from pylab import plot, savefig

R = 10.
M = 1.
G = 1.
phi = 30 * 180. / pi
v0_magnitude = 0.1
r = array([0, 0, R])
v = v0_magnitude * array([cos(pi), 0, sin(pi)])
trajectory = []

dt = 0.1
t = 0
while t < 100:
trajectory.append(r.copy())
a = -G * M * r / norm(r)**3
v += a * dt
r += v * dt
t += dt

trajectory = vstack(trajectory)
plot(trajectory[:, 0], trajectory[:, 2])
savefig("trajectory.png")


And here is the graph:

I have also uploaded the above code and the graph at this gist.

• Cool that one can plot directly in Python. I already built something similar, but mine just prints the trajectory. – Martin Ueding Jan 2 '12 at 19:53
• It should be noted (for completeness's sake) that such a simple solution method will only work because the equations of motion here are so simple. For a more complicated system, or for a system with constraints, this wouldn't be possible, and a more complicated method (such as DASSL or ode15s in MATLAB, for instance) would need to be used. Also, to plot in Python, one needs to install the matplotlib package in Python, which can be a pain. I'd recommend the Enthought Python Distribution for an integrated solution with precompiled libraries. – aeismail Jan 2 '12 at 20:18
• Yes, the Euler method is not very precise -- the graph should be an exact ellipse, but as you can see, the second loop doesn't run on top of the first one due to numerical errors. You can implement rk4, or use some solver from SciPy if more precise method is needed. Or just decrease the 'dt' a lot. – Ondřej Čertík Jan 2 '12 at 23:04
• If you want it to be an ellipse, use a symplectic method. You'll still have phase error, but you can't see that in this kind of plot. – David Ketcheson Jan 3 '12 at 6:18
• If I use a smaller $\mathrm dt$, it gets an ellipse for all practical purposes. The initial question was just whether the projectile would hit the planet again, and this did the job. – Martin Ueding Jan 3 '12 at 11:45

One option to look at for systems that involve evolution over time would be to consider using either MATLAB or Modelica. MATLAB is of course proprietary, but Modelica exists in many implementations, some of which, such as OpenModelica, are definitely open-source.

As someone who is relatively new to Modelica (I only learned of it about a year ago), I find that it is very good at implementing problems based on dynamical systems: sets of ordinary differential equations or differential-algebraic equations. In particular, the setup of problems is highly modularized, and you can write equations in "natural" form, rather than having to rewrite them in arbitrary forms to massage them into whatever format is required by a language or implementation such as what MATLAB requires. OpenModelica also offers some integrated graphing capabilities.

• I have octave here, which says to be matlab compatible. When I search the Debian index for modelica, I get scilab, which I thought was a matlab clone as well. Could one simulate this in octave? – Martin Ueding Jan 2 '12 at 16:16
• Yes, it should be possible. Octave supports solving DAE's, so this kind of system should be solvable as well. Just look up solving ODE's and DAE's in your Octave user documentation. As for Modelica, as I mentioned, it's implemented in a number of different packages, including one that's wrapped into Mathematica. So many solution routes are possible. – aeismail Jan 2 '12 at 18:04

Both Maple and Mathematica have easy interfaces to numerically solve differential equations and plot their results. Here http://www.math.tamu.edu/~bangerth/teaching/2010-fall-442/project-tools.pdf http://www.math.tamu.edu/~bangerth/teaching/2010-fall-442/project-tools.mw is an example with Maple, complete with cool graphs.