Numerical solutions of sixth order differential equations

An object moving in three dimensions is described by six variables: $x$, $y$, $z$, $v_x$, $v_y$, and $v_z$. If the force on the object is known, then the motion can be described by six first order differential equations,

$\large \frac{dx}{dt}=v_x$   $\large \frac{dv_x}{dt}=F_x(x,y,z,v_x,v_y,v_z,t)/m$

$\large \frac{dy}{dt}=v_y$   $\large \frac{dv_y}{dt}=F_y(x,y,z,v_x,v_y,v_z,t)/m$

$\large \frac{dz}{dt}=v_z$   $\large \frac{dv_z}{dt}=F_z(x,y,z,v_x,v_y,v_z,t)/m$

Here $F_x$, $F_y$, and $F_z$ are the three components of the force, $m$ is the mass, and $t$ is the time. The form below can be used to numerically integrate these equations for a total number of $N_{steps}$ steps using a step size of $\Delta t$. The force can be specified as a function of $x$, $y$, $z$, $v_x$, $v_y$, $v_z$ and $t$.

 Numerical 6th order differential equation solver 

$ \large \frac{dx}{dt}=$

$v_x$

$ \large \frac{dv_x}{dt}=$

$ \large \frac{dy}{dt}=$

$v_y$

$ \large \frac{dv_y}{dt}=$

$ \large \frac{dz}{dt}=$

$v_z$

$ \large \frac{dv_z}{dt}=$

Initial conditions: 

$t_0=$

 

$\Delta t=$

$x(t_0)=$

 

$N_{steps}$

$v_x(t_0)=$

 

Plot:

vs.

$y(t_0)=$

 

$v_y(t_0)=$

 

$z(t_0)=$

 

$v_z(t_0)=$

 

 

 $t$   $x$   $v_x$   $y$   $v_y$   $z$   $v_z$

Numerical integration is based on the idea that if the initial conditions for $\vec{x}(t_0)$ and $\vec{v}(t_0)$ are known, a good estimate for $\vec{x}$ and $\vec{v}$ a short time $\Delta t$ later is,

$\large \vec{x}(t_0+\Delta t) \approx \frac{d\vec{v}}{dt}|_{t_0}\Delta t$ and $\large \vec{v}(t_0+\Delta t) \approx \vec{F}(x(t_0),v_x(t_0),t_0)\Delta t/m.$

Once an estimate for the position and the velocity of the object at time $t_0 + \Delta t$ is calculated, they can be used to estimate the position and the velocity at time $t_0 + 2\Delta t$.

In the form above, the force can be specified using the standard mathematical functions abs(x), acos(x), asin(x), atan(x), cos(x), exp(x), pi = 3.141592653589793, pow(x,y) = xy, round(x), sin(x), sqrt(x), and tan(x). Note that multiplication must be specified with a '*' symbol, 3*cos(x) not 3cos(x). Powers are specified with the 'pow' function: x² is pow(x,2) not x^2.

The numerical integration will become unstable if the time step is too long. The time step should be a factor of 10 to 100 times smaller than any characteristic time of the motion that is being calculated. If the orbit of the earth around the sun is being calculated, then a time step of 3 days (about 4E6 seconds) is reasonable. If the motion of an electron moving past a ion is being calculated, a time step of 1 ps is reasonable. If the routine becomes unstabe, nothing will be plotted. In this case, try a shorter time step.

There are many routines available to perform numerical integration. The procedure outlined above is called the Euler method. However this page really uses the more accurate fourth order Runge-Kutta method with a fixed step size. Some numerical integration routines calculate the optimum time step automatically.