PHY.K02UF Molecular and Solid State Physics

Numerical solutions of second order differential equations

An object moving in one-dimension can be described in terms of its position $x$ and and its velocity $v_x$. If the force on the object is known, then the motion can be described by two first order differential equations,

$\large \frac{dx}{dt}=v_x$ and $\large \frac{dv_x}{dt}=a_x=F_x(x,v_x,t)/m.$

Here $F$ is 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 acceleration $a_x$ can be given as a function of $x$, $v_x$, and $t$.

 Numerical 2nd order differential equation solver 

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

$v_x$

$ \large a_x=\frac{F_x}{m}=\frac{dv_x}{dt}=$

Intitial conditions: 

$x(t_0)=$

 

$\Delta t=$

$v_x(t_0)=$

 

$N_{steps}$

$t_0=$

 

Plot:

vs.

 

 $t$       $x$      $v_x$

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

$\large x(t_0+\Delta t) \approx \frac{dv_x}{dt}|_{t_0}\Delta t$ and $\large v_x(t_0+\Delta t) \approx 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.