Processing math: 3%

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.

$v_x$
 
-1.0
-0.5
0.0
0.5
1.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5

$x$

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.