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) = x^y, 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.

Example 1: Drag

A ball of mass $m$ is thrown and experiences a frictional drag force as it moves through a gas or a fluid. The forces acting on this ball are gravity $-mg\hat{z}$ the drag force. The drag force can often has the form,

$\large \vec{F}_{fric} = -a\vec{v} - b\vec{v}|\vec{v}|,$

where $a$ and $b$ are constants. For low Reynolds number, the linear term $-a\vec{v}$ usually dominates whereas for high Renolds number, the quadratic term $-b\vec{v}|\vec{v}|$ dominates.

If the ball falls for a long time the velocities in the $x-$ and $y-$directions will approach zero while the velocity in the $z-$direction will approach a constant called the terminal velocity when the drag force balances the gravitational force, $-mg-av_z - bv_z|v_z|=0$. Solving this equation for $v_z$ yields the terminal velocity.

If $a=0$ and $b=0$, the ball will follow a parabolic trajectory and analytic solutions to the differential equation are possible. If $a>0$ and $b=0$, there is linear damping and it is also possible to find an analytic solution.

Example 2: A charged particle in constant electric and magnetic fields

A charged particle moves in a constant electric field $\vec{E}$ and a constant magnetic field $\vec{B}$. The force on the electron is,

$\large \vec{F} = q(\vec{E} + \vec{v}\times \vec{B}),$

where $q$ is the charge of the particle and $m$ is the mass.

When the electric field is perpendicular to the magnetic field, the average velocity of the charged particle will be in a direction perpendicular to both the electric field and the magnetic field.

Example 3: Satellite

A satellite orbits the earth. The graviational force on the satellite is,

$\large \vec{F} = -\frac{Gm_e m_s}{r^2} \hat{r}$,

where $G= 6.6726 \times 10^{-11}$ N m²/kg² is the graviational constant, $m_e = 5.97219 \times 10^{24}$ kg is the mass of the earth, $m_s$ is the mass of the satellite, and $\vec{r}$ is the position of the satellite measured from the center of the earth.

If the orbit falls below about 6400000 m, the satellite will crash into the earth. There are various kinds of orbits such as geosynchronous orbits, geostationary orbits, low earth orbits, elliptical orbits, and graveyard orbits. The difference just depends on the initial conditions of the satellite.