Numerical integration, with examples

Table of contents:

Introduction

 Back to top

If is often the case that we want to write computer programs to simulate various physics processes. For example, we might want to write a program that shows the path of a projectile through the air, including air resistance. The projectile motion defines a plane and so will have 2 coordinates, $x$ and $y$ for horizontal and vertical. What the program will do is to simulate the dynamics, and perhaps even to draw $x(t)$, $y(t)$ (or equivalently, $y(x)$) with some initial conditions. In our ballistics example, a ball thrown up into the air will have some initial velocity vector $\vec {v_0}$ (initial speed $v_0$ at angle $\theta_0$), and an initial position $(x_0,y_0)$. If we take our initial position to be at the coordinate $(x_0,y_0)=(0,0)$, then our initial conditions will be: $$\begin{align} x(t\!=\!0) & =0 \nonumber \\ y(t\!=\!0) & =0 \nonumber \\ \dot x(t\!=\!0) & =v_{0x} \nonumber \\ & =v_0\cos\theta_0 \nonumber \\ \dot y(t\!=\!0) & = v_{0y} \nonumber \\ & =v_0\sin\theta_0 \nonumber \end{align}\nonumber$$

What we want to do is to write a program that would start at the initial conditions, and "swim" the projectile to the new position, and plot it. And repeat as appropriate. This is easy to do conceptually, we just use the differential equation that governs the motion to do the propagation, or "swimming".

To be general, imagine we have an equation that governs the motion such as this: $\dot y(t) = f(y,t)$, and we want to "swim" from some value $t_s$ to some other value $t_f$. Using the calculus definition of $dy$, we can write $dy \equiv y(t_f) - y(t_s)$, with $t_f = t_s + \delta t$, write the equation of motion as $dy = f(y,t)\cdot dt$, and integrate it over the interval between $t_s$ and $t_f$. This gives us the equation

$$y(t_f)=y(t_s) + \int_{t_s}^{t_f}f(y,t)dt\label{eintegrate}$$

Given that sometimes we have to deal with functions $f(y,t)$ that are not analytically integrable, we have to do the integration numerically on a computer. Voila, now you see why the subject of simulating physical systems where all you have are the initial conditions and rate of change is called "numerical integration": it all depends on integrating the differential equation over some interval that you choose.

The trick is to turn $$\int_{t_s}^{t_f}f(y,t)dt\nonumber$$

into something you can use a computer to solve.

Note however that sometimes we start with a differential equation that is second order in the time derivative: $$\ddot y(t) = f(y,t)\nonumber$$ How do we then come up with a way to integrate to get the change in position? The answer is wonderful: we can take any 2nd order differential equation and break it down into 2 1st order differential equations by implementing the following "trick":

Define $x\equiv y$ and $z\equiv\dot y$, and get the following 2 equations: $$\begin{align} \dot z & = f(y,t) = f(x,t) \nonumber \\ \dot x & = \dot y = z \nonumber \end{align}\nonumber$$ Now we have 2 coupled differential equations: $$\begin{align} \dot z & = f(x,t) \nonumber \\ \dot x & = z \nonumber\end{align}$$ We can then use the above perscription and integrete to "swim" $x$ and $z$: $$\begin{align} z(t+\dt) & =z(t)+\dt\cdot \dot z(t) \nonumber \\ & =z(t)+\dt\cdot f(x,t) \nonumber \\ x(t+\dt) & =x(t)+\dt\cdot \dot x(t)\nonumber \\ & =x(t)+\dt\cdot z(t)\nonumber \\ \end{align}\nonumber $$ which is equivalent to the following 2 equations, going back to the original variable $y(t)$: $$\begin{align} \dot y(t+\dt) & = \dot y(t) + \dt\cdot f(y,t)\label{eq1} \\ y(t+\dt) & = y(t) + \dt\cdot \dot y(t)\label{eq2} \end{align}$$


Euler method

 Back to top

This method is based on the simple idea that if you make the interval $[t_s,t_f]$ small enough ($\delta t$ is small), then if $f(y,t)$ is approximately constant over the interval, solving equation $\ref{eintegrate}$ is relatively easy: $$\begin{align} y(t_f) & = y(t_s)+\int_{t_s}^{t_f}f(y,t)dt \nonumber \\ & = y(t_s)+\dt\cdot f(y_s,t_s) \nonumber \end{align}\nonumber $$ where $\dt\equiv t_f - t_s$, and $\delta t$ is small enough so that $f(y_s,t_s)\sim f(y_f,t_f)$ over that interval. The condition for this approximation to be "accurate" is that $f(y,t)$ doesn't change "fast" over the interval $\dt$. More on this later.

This approximation allows us to write code that will start at the initial conditions, step through a time $\dt$, and calculate $y(t)$ and $\dot y(t)$ each time so that we can plot the motion.

But what does "a small enough interval" mean? What's "good enough"? To some degree it might involve one's sensibilities, but we can also use one of the most powerful things in physics to help guide us: energy conservation! More on that later.


Euler method and projectile motion

 Back to top

Let's take the simple example of projectile motion in a constant gravity field ("ballistic" motion, where you start with an initial energy and let it fly). Here the vertical acceleration is constant, so that $\ddot y = -g$ ($y$ measured positive up, acceleration is down), with initial conditions $y(0) = 0$ (start from the ground) and $\dot y(0) = v_0$. (We are ignoring the horizontal motion for now, but it's pretty easy since there is no horizontal acceleration.)

This motion is pretty easy to integrate, and we know the real (analytic) answer: $$\dot y(t) = v_0 - gt\label{eq3}$$ $$y(t) = v_0 t - \half g t^2\label{eq4}$$ The energy will be given by $E = \half m\dot y^2 + mgy$, and since energy is conserved and the initial energy is $E_0 = \half mv_0^2$, we should be able to verify this directly:

$$\begin{align} E & = \half m\dot y^2 + mgy \nonumber \\ & = \half m(v_0-gt)^2 + mg(v_0t-\half gt^2) \nonumber \\ & = \half mv_0^2-mgv_0t+\half mg^2t^2 + mgv_0t-\half mg^2t^2 \nonumber \\ & = \half mv_0^2$\nonumber \\ \end{align}\nonumber $$ As expected!

But let's pretend we don't know how to integrate $\ddot y=-g$, and instead solve it numerically.

So if we start at $t=0$, and use equations (\ref{eq1}) and (\ref{eq2}), we would have after the 1st iteration at time $\dt$: $$\begin{align} \dot y(\dt) & = \dot y(0) + \ddot y(0)\cdot\dt \nonumber \\ & = v_0 - g\dt\label{eq5} \\ y(\dt) & = y(0) + \dot y(0)\cdot\dt \nonumber \\ &= v_0\dt\label{eq6}\end{align}$$ We can calculate how far off this approximation is from the correct calculation we get from directly integrating $\ddot y=-g$, but we can also see it by calculating the energy after a time $\dt$: $$\begin{align} E(\dt) & = \half m\dot y^2 + mgy \nonumber \\ & = \half m(v_0-g\dt)^2 + mgv_0\dt \nonumber \\ & = \half mv_0^2-mgv_0\dt+\half mg^2\dt^2 + mgv_0\dt \nonumber \\ & = \half mv_0^2 + \half mg^2\dt^2\nonumber\end{align}$$ We can see that this is off by an amount $\half mg^2\dt^2$ (2nd term in the last equation), which is also the amount that equation (\ref{eq6}) differs from equation (\ref{eq4}).

Taking another step in $\dt$ we get: $$\begin{align} \dot y(2\dt) & = \dot y(\dt) + \ddot y(\dt)\cdot\dt \nonumber \\ & = v_0 - g(2\dt) \nonumber \\ y(2\dt) & = y(\dt) + \dot y(\dt)\cdot\dt \nonumber \\ & = v_0(2\dt)-\frac{1}{4}g(2\dt)^2\nonumber \end{align}$$ Again we see that the equation for $y(2\dt)$ is off,this time by an amount $\frac{1}{4}g(2\dt)^2$. When we calculate the energy we get $$\begin{align} E(2\dt) & =\half mv_0^2 + mg^2\dt^2 \nonumber \\ & = \half mv_0^2 + \frac{1}{2^2}mg^2(2\dt)^2\nonumber\end{align}$$ For $3\dt$, we get $$\begin{align} \dot y(3\dt) & = \dot y(2\dt) + \ddot y(2\dt)\cdot\dt \nonumber \\ & = v_0 - g(3\dt)\nonumber \\ y(3\dt) & = y(2\dt) + \dot y(2\dt)\cdot\dt \nonumber \\ & = v_0(3\dt)-\frac{1}{3^2}g(3\dt)^2\nonumber\end{align}$$ and the energy is given by $$\begin{align} E(3\dt) & = \half mv_0^2 + \frac{3}{2}mg^2\dt^2 \nonumber \\ & = \half mv_0^2 + \frac{1}{6}mg^2(3\dt)^2\nonumber \end{align}$$ You can see that the energy per step is diverging, getting further away from energy conservation by adding another $\frac {1}{2}mg^2\dt^2$ at every step. So the bottom line is that the Euler technique is not so good - it violates energy conservation! Except it violates is to second order in $\dt$, so if you make $\dt$ small, you are still wrong, but only a little bit.


Modified Euler method (aka Euler-Cromer)

 Back to top

One obvious variation in the Euler technique is to use the slope at the end of the interval instead of at the beginning. The equations below show this explicitly: $$\dot y(t+\dt) = \dot y(t) + \dt\cdot f(y,t)\label{eq7}$$ $$y(t+\dt) = y(t) + \dt\cdot \dot y(t+\dt)\label{eq8}$$ This approximation is also imperfect, since as shown in the figure above, the slope at the end of the interval overestimates the average slope over the interval, and as you would expect, the value of $y(t+\dt)$ overestimates the true value. And as above, energy is not conserved. But for some examples, Euler-Cromer is better than straight Euler.


Modified Euler method (aka Runge-Kutta)

 Back to top

So both Euler and Euler-Cromer are problematic. What to do? In our example above ($y=ax^3$), we found that using the slope at the beginning of the interval to swim to the endpoint (Euler approximation) will undershoot the correct answer, and using the slope at the endpoint (Euler-Cromer) will overshoot. The logical thing to do is to consider using the value of the slope at the midpoint of the interval, $\dot y(t+\half \dt)$, as a better approximation to the average slope over the interval.

Some definitions that will help with the notation and formulas below:

So if we know $y_s$ and want to swim to $y_f$, the RK2 approximation is to use $\dot y(t_h)$, the slope at the halfway point $t_h$:

$$y_f = y_s + \dt\cdot \dot y(t_h)$$

The problem is that we know $\dot y(t_s)$, but we don't know $\dot y(t_h)$. But we do have the differential equation that governs how $\dot y(t)$ behaves, and we have the Euler approximation. So we can make a 2 step process:

  1. Use the Euler approximation to swim to the halfway point:
    $$y_h = y_s + \half\dt\cdot\dot y(y_s,t_s)$$ $$\dot y(t_h) = \dot y_s + \half \dt\cdot f(y_s,t_s)\label{eq9}$$
  2. Then use the RK2 approximation involving the slope at the halfway to swim to the final point $y_f$:
    $$y_f = y_s + \dt\cdot \dot y_h\label{eq10}$$ $$\dot y_f = \dot y_s + \dt\cdot f(y_h,...)\label{eq11}$$
Equations (\ref{eq10}) and (\ref{eq11}) describe what is called "Runge-Kutta, second order", aka RK2. In words, we swim the position and slope to the halfway point $t_{1/2}$, and use those to calculate how to swim (integrate) from $y(t)$ to $y(t+\delta t)$.

How well will the RK2 method work? Let's combine equations (\ref{eq9}) and (\ref{eq10}) to get: $$\begin{align} y(t+\dt) & = y(t) + \dt\cdot \dot y(y_h,t_h) \nonumber \\ & = y(t) + \dt\cdot [\dot y(t)+\half \dt\cdot f(y_h,t_h)] \nonumber \\ & = y(t) + \dt\cdot \dot y(t) + \half \dt^2\cdot f(h_h,t_h) \nonumber \end{align}$$ There's your answer: RK2 will be more accurate than Euler because it's really a 2nd order correction to Euler. In fact, it's really nothing more than an expansion of $y(t+\delta t)$ around $\delta t$, keeping the terms up to order $\delta t^2$. As you will see below, sometimes RK2 is all you need to do to get reasonable results.

To illustrate the limitations of the 3 method graphically, the figure below shows an arbitrary function $y=ax^3$ in blue. The starting and ending points are labeled $y(t)$ and $y(t+\delta t)$ are drawn at times $t$ and $t+\delta t$ respectively, drawn with solid black circles. The yellow line shows the slope calculated at some time, and drawn starting at the point $y(t)$. The radio buttons allow you to specify which of the 3 times you want to use for the slope: $\dot y(t)$ for Euler,$\dot y(t+\delta t)$ for Euler-Cromer, and $\dot y(d+\delta t)$for Runge-Kutta (2). The value of $y(t+\delta t)$ using one of the 3 approximations is drawn with an open blue circle. By changing the selection and hitting the "Redraw" button, you can see how close each of the approximations are to the true value (the one on the curve). You can select The Euler method starts at the beginning of the interval, at $y(t)$, and uses the slope at the beginning of the interval, $\dot y(t)$, to integrate the function (to swim to the end of the interval at $y(t+\dt)$). You can also grab point $y(t+\delta t)$ with the mouse (or touch, on a touch pad or smart phone) and change the value for $\delta t$, and see how much the approximation changes as a function of $\delta t$.

As you can see in the figure, when the slope is changing continually over the interval, using the slope at the beginning is not such a great approximation: the slope at the beginning, $\dot y(t)$, underestimates what the average slope would be over the interval, and so the value of $y(t+\dt)$ underestimates the true value at the endpoint, which is what we are after. And of course, this is what is contributing to violation of energy conservation. Similarly, using $\dot y(t+\delta t)$ overestimates the approximation. Using Runge-Kutta (2) allows much better accuracy, as long as $\delta t$ doesn't get too large.

Click and drag points $y(t)$ and $y(t+\dt)$
Euler Euler-Cromer Runge-Kutta (2)

The next figure shows what is happening visually, where the gray line is the slope at the midpoint, and the yellow line is what you get when you "swim" from $y(t)$ to $y(t+\dt)$ using the slope at the midpoint. As you can see from the fact that the two circles at $t+\dt$ (filled in is exact, open is from the technique) are very close compared to the figure above.


RK2 Formally

 Back to top

As discussed above, often in physics we will start with a 2nd order differential equation, not a first order one, for instance $\ddot x=f(x,t)$. It is easy to use RK2 here, with the following substitution:

$$\alpha \equiv \dot x\nonumber$$ $$\beta = x\nonumber$$ Then our equation of motion is $$\dot\alpha = f(\beta,t)\nonumber$$ $$\dot\beta = g(\alpha)=\alpha\nonumber$$

The RK2 method for $\alpha$ and $\beta$ would be: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f(\beta_h,t_h)\label{eqnalpha}$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot g(\alpha_h,t_h)\label{eqnbeta}$$ with $$\alpha_h\equiv\alpha(t+\half \dt)=\alpha(t)+\half \dt\cdot f(\beta)\label{eqnalphah}$$ $$\beta_h \equiv \beta(t+\half \dt)=\beta(t)+\half \dt\cdot g(\alpha)\label{eqnbetah}$$ $$t_h\equiv t+\half \dt\label{eqnth}$$ We are now ready to tackle some real processes.


RK2 and projectile motion

 Back to top

In the case of projectile motion, we have the two equations: $$\ddot y = -g\nonumber$$ $$\ddot x = 0\nonumber$$ with initial conditions $x(0)=0$, $y(0)=0$, $\dot x(0)=v_{0x}$ and $\dot y(0)=v_{0y}$. Motion along the horizontal $x$ direction is trivial, since zero acceleration means constant velocity. The equation for $x(t+\dt)$ is:

$$\begin{align} \dot x(t+\dt) & = \dot x(t) = v_{0x} \nonumber \\ x(t+\dt) & = x(t) + v_{0x}\dt\nonumber\end{align}$$ For vertical motion we use the RK2 equations as above. Just to show how to follow the perscription, we will reduce the equation 2nd order equation $\ddot y=-g$ to two first order equations using the substitutions $$\alpha \equiv \dot y\nonumber$$ $$\beta \equiv y\nonumber$$ The equations of motion are then: $$\dot\alpha = f(\beta)\equiv -g\nonumber$$ $$\dot\beta = g(\alpha)\equiv \alpha\nonumber$$ Next we calculate $\alpha_h$ and $\beta_h$: $$\begin{align} \alpha_h\equiv\alpha(t+\half\dt) & =\alpha(t)-\half g\dt \nonumber\\ \beta_h\equiv\beta(t+\half\dt) & =\beta(t)+\half\dt\cdot\alpha(t)\nonumber\end{align}$$ and plug into equations (\ref{eqnalpha}) and (\ref{eqnbeta}): $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f(\beta_h)=\alpha(t)-g\dt\nonumber$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot g(\alpha_h)=\beta(t)+\dt\cdot\alpha(t)-\half g\dt^2\nonumber$$ You can see explicitly here how the 2nd order correction for $\beta$ (which is also $y$) now comes into play. Switching back to the varialbles $y$ and $\dot y$ gives: $$\begin{align} \dot y(t+\dt) & =\dot y(t)-g\dt \nonumber \\ y(t+\dt) & =y(t)+\dt\cdot \dot y(t) - \half g\dt^2\nonumber\end{align}$$ We now have something we can code: the new variables $y$ and $\dot y$ at each interval are functions of the old variables and the interval $\dt$.

It is worth doing this explicitly once, starting at $t=0$ and using the initial conditions as above, we can swim to $y(\dt)$ using equations (\ref{eq9}) and (\ref{eq10}) above: $$\begin{align} \dot y(\dt) & =\dot y(0) -g\dt \nonumber\\ & = v_0-g\dt\nonumber\\ y(\dt) & =y(0)+ \dt\cdot \dot y(\half \dt) \nonumber\\ & = \dt(v_0\!-\!\half g\dt) \nonumber\\ & = v_0\dt - \half g\dt^2\nonumber\end{align}$$ Voila, this looks like the exact analytical equations!

Going to the next interval, $t=\dt$ (so $t+\dt=2\dt$) gives: $$\begin{align} \dot y(2\dt) & =\dot y(\dt) -g(2\dt) \nonumber\\ & = v_0-g(2\dt)\nonumber\\ y(2\dt) & =y(\dt)+ \dt\cdot\dot y(\frac{3}{2}\dt) \nonumber\\ & = v_0\dt - \half g\dt^2 + \dt(v_0 - \frac{3}{2}g\dt) \nonumber\\ & = v_0(2\dt) -2g(\dt)^2 = v_0(2\dt) -\half g(2\dt)^2\nonumber\end{align}$$ which is also exactly correct. One can then calculate all other intervals, and we find that they are also correct. This particular technique works much better.

Why does this work so well for motion in a constant gravitational field? Think back about what the equation $\ddot y=-g$ is telling us: the 2nd derivative is constant, so that means that the first derivative is changing linearly with time. Which means that making the approximation of using the derivative at the midpoint of the interval is actually exactly the right thing to do in this situation.

In the graph below, the solid line shows the exact solution of the parabolic path for a particle in a constant downward gravitational field. For comparision with the numerical integration, you can select using the radio buttons either: 1) the Euler technique (also known as the Explicit Euler); 2) the Cromer-Euler (also known as the Implicit Euler); or 3) the Runge-Kutta RK2. When you change the selection, remember to hit the Start button again. The Euler technique draws the path in blue circles, the Euler-Cromer in yellow, and the RK2 in red. Below that is a graph of the total energy, with arbitrary units that bring out the difference. You can easily see how the 3 techniques differ.

Euler Euler-Cromer Runge-Kutta 2

$|\vec {v_0}|=$ 90
$\theta_0=$ 65

Energy (arbitrary units):


RK2 and projectile motion with air resistance

 Back to top

In the above example, we started with a simple (constant) 2nd derivative, $\ddot y = -g$. The 1st derivative would then be a function of time $t$, and so on for $y(t)$. However nature often throws us more difficult equations to integrate. The next one to study would be equations where the differentials are functions of external variables like $t$ and also of themselves, e.g. $\dot z = f(z(t),t)$.

For example, ballistic motion through the air should take into account air resistance. If the velocity is not too high, then the air resistance will be linear (against) the velocity, and so the equation for the motion in the vertical direction would be something like

$$\ddot y = -g - \omega_0 \dot y\nonumber$$

where the constant $\omega_0 \equiv k/m$, and $k$ is the proportionality constant for the air resistance such that the force is given by $\vec{F_r}=-k\vec{\dot y}$.

Following the same line of reasoning that led to the RK2 method, we first define $z\equiv \dot y$ and extend equations (\ref{eq7}) and (\ref{eq8}) by defining $z(t+\dt) = z(t) + \dt\cdot \dot z(z_h,t_h)$. That is, we swim from $z(t)$ to $z(t+\dt)$ by using the slope $\dot z(z_h,t_h)$ evaluated at the midpoint $z_h, t_h$ (instead of evaluating $\dot z$ at the beginning of the interval, $\dot z(z,t))$, as with the Euler method). So, from the previous RK2 chapter above, we get $$z_h \equiv z(t+\half \dt) = z(t) + \half \dt\cdot \dot z(t)\nonumber$$ That is, we use the slope at the beginning of the interval to get us to the midpoint ($z_h,t_h$), and then use the slope at the midpoint to get us to the endpoint ($z(t+\dt),t+\dt$): $$\begin{align} z(t+\dt) & = z(t) + \dt\cdot \dot z(z_h,t_h)\nonumber\\ & = z(t) + \dt\cdot \dot z(z(t)+\half \dt\cdot \dot z(t))\nonumber\end{align}$$ Or in the usual notation for Runge-Kutta: $$z(t+\dt) = z(t) + \dt\cdot \dot z(z(t)+\half k_1)\label{eqn11}$$ where $$k_1 \equiv \dt \dot z(z(t))\nonumber$$ And of course, we do the same thing for $y(t+\dt)$: $$y(t+\dt) = y(t) + \dt\cdot \dot y(y(t) + \half k_1)\label{eqn12}$$ where $$k_1 = \dt\dot y(t)\nonumber$$ Note that there is no explicit $t$ dependence in the equation for $\ddot y$.

We are now ready to simulate projectile motion with air resistance, where the force from air resistance is linear in the velocity (this is true as long as the velocity is not "large" compared to the limiting viscosity of the air): $\vec {F_r}=-k\vec v$. So there will be a resistance in both the $x$ and $y$ direction, which means the equations for the acceleration will be: $$\begin{align} \ddot y & = -g-\frac{k}{m}\dot y\nonumber\\ \ddot x & = -\frac{k}{m}\dot x\nonumber\end{align}$$ To make things simpler, define $\tau \equiv \frac{m}{k}$, and note that $\tau$ has the units of time. Then we can write the acceleration equations as: $$\begin{align} \ddot y & = -g-\dot y/\tau\nonumber\\ \ddot x & = -\dot x/\tau\nonumber\end{align}$$ For the initial conditions, assume $y(0)=x(0)=0$, the initial angle of the velocity is $\theta_0$, and then $\dot y(0)=v_{0y}=v_0\sin\theta_0$ and $\dot x(0)=v_{0x}=v_0\cos\theta_0$.

It turns out that both of these equations can be solved analytically in a very straight forward manner by direct integration and application of the initial conditions. Starting with the horizontal motion, we get: $$\begin{align} \dot x(t)&=v_{0x}e^{-t/\tau}\nonumber\\ x(t) & = v_{0x}\tau (1-e^{-t/\tau})\nonumber\end{align}$$ Note that $v_{0x}\tau$ is the horizontal position at $t=\infty$ so we can define $x_\infty\equiv v_{0x}\tau$ for later and write $$x(t) = x_\infty (1-e^{-t/\tau})\nonumber$$ So at $t=0$, the projectile starts off with some initial velocity, and as it increases, the air resistance starts to dominate ultimately limiting the horizontal distance traveled.

Solving for time $t$ (we will need this later) gives $$t = -\tau\ln (1-x/x_\infty)\nonumber$$ For the vertical motion, we can integrate the vertical acceleration $\ddot y(t)$ to get: $$\dot y(t) = v_{0y}e^{-t/\tau}+g\tau(1-e^{-t/\tau})\nonumber$$ Note that as $t\to\infty$, $y(t)\to g\tau$ which means we can interpret $g\tau$ as the terminal velocity in the vertical direction: $v_{ter}\equiv g\tau$. So we can write the equation for $\dot y(t)$ as: $$\dot y(t) = v_{0y}e^{-t/\tau}+v_{ter}(1-e^{-t/\tau})\nonumber$$ To get the horizontal motion we integrate the above equation to get: $$y(t) = v_{ter} t + (v_{0y}+v_{ter})\tau(1-e^{-t/\tau})\nonumber$$ Putting this altogether, here are the equations of motion in a convenient form. $$\begin{align} \dot x(t) & =v_{0x}e^{-t/\tau}\nonumber\\ x(t) & = v_{0x}\tau (1-e^{-t/\tau})\nonumber\\ \dot y(t) & = v_{0y}e^{-t/\tau}+v_{ter}(1-e^{-t/\tau})\nonumber\\ y(t) & = v_{ter} t + (v_{0y}+v_{ter})\tau(1-e^{-t/\tau})\nonumber \end{align}$$ Solving for $y(x)$ gives: $$y(x) = x\frac{v_{0y}+v_{ter}}{v_{0x}}+v_{ter}\tau \ln(1\!-\!x/x_\infty)\nonumber$$ Notice that at small $x$, $y\sim x$, but as $x$ increases, the logarithm term begins to dominate, and becomes negative (log of a number less than 0 is negative). In the figure below, you can change the ratio of $k/m$, and rescale the plot so that you can see the full path. As $k/m$ gets large, you can see the sharp fall off of the position: the air is just not letting it get very far!

To set up the numerical approximation, for the vertical motion the Euler method would then give: $$\begin{align} \dot y(t+\dt) & = \dot y(t) + \dt\cdot \ddot y(\dot y(t)) \nonumber\\ & = \dot y(t) - \dt\cdot (g+\dot y(t)/\tau)\nonumber\\ y(t+\dt) & =y(t)+\dt\cdot \dot y(t)\nonumber\\ \end{align}$$ and for the horizontal motion: $$\begin{align} \dot x(t+\dt) & = \dot x(t) + \dt\cdot \ddot x(\dot x(t)) \nonumber\\ & = \dot x(t) - \dt\cdot \dot x(t)/\tau \nonumber\\ x(t+\dt) & =x(t)+\dt\cdot \dot x(t)\nonumber\\ \end{align}$$ The computer code would look something like what's below. The parameters are gravity=9.8 and tau=m/k, the interval $\dt$ is dt, and of course (y,ydot)=$(y,\dot y)$. You start with some initial value for y and ydot and iterate:

		yaccel = -gravity - ydot/tau;
		ydot_new = ydot + dt*yaccel;
		y_new = y + dt*ydot;
		ydot = ydot_new;
		y = y_new;
		//
		xaccel = -xdot/tau;
		xdot_new = xdot + dt*xaccel;
		x_new = x + dt*xdot;
		xdot = xdot_new;
		x = x_new;
This solution will have the same problem (inaccurate and violates energy conservation) as the above simple example ($\ddot y=-g$).

The modified (aka implicit) Euler formula (aka Euler-Cromer) differs only by which version of the slope you use, so that the only difference is in the equation for the position $y(t)$:

Euler $y(t+\dt)=y(t)+\dt\cdot \dot y(t)$
Euler-Cromer $y(t+\dt)=y(t)+\dt\cdot \dot y(t+\dt)$

and the computer code would differ from the block above only in the calculation of y_new and x_new

		.
		.
		y_new = y + dt*ydot_new;
		.
		.
		//
		.
		.
		x_new = x + dt*xdot_new;
		.
		.
The RK2 would be the following:

$$\begin{align} \dot x_h & \equiv \dot x(t) + \half \dt\cdot\ddot x(\dot x(t))\nonumber \\ \dot x(t+\dt) & = \dot x(t) + \dt \cdot\ddot x(\dot x_h)\nonumber \\ x(t+\dt) & = x(t) + \dt \cdot\dot x_h$\nonumber \\ \end{align}$$ In words, we "swim" the slope to the halfway mark between $t$ and $t+\dt$, forming $\dot x_h$, and then use that in place of $\dot x$ in the Euler equations above. This works rather well, as you can see from the interactive figure below. You can use the sliders to change the initial conditions, the "Rescale" will allow the curve to fit in the window, and you can toggle the air resistance on and off (on will draw a black curve, off will draw a white one).


Euler Euler-Cromer Runge-Kutta 2
$k/m$ 0.1
$|\vec {v_0}|$ 130
$\theta_0$ 60
Rescale 1
Air Resistance: None $\vec{F_r}=k\vec v$

Motion of a pendulum

 Back to top

The simple pendulum consists of a mass $m$, a length $L$, and angle $\theta$ measured with respect to the vertical downward direction. The diagram below details the coordinates and parameters.

It's easy to use Newton's law to calculate the force components, but it's also easy to use Lagrangians, and this will warm you up for the inevitable consideration of the double pendulum.

We take the coordinate origin $(0,0)$ to be the point where the pivot is connected to the support. Then from the diagram, the cartesian coordinates of the mass $m$ are given by: $$(x,y) = [L\sin\theta,-L\cos\theta]\nonumber$$ The derivatives, which will be needed later, are given by $$(\dot x,\dot y) = (L\dot\theta\cos\theta,L\dot\theta\sin\theta)\nonumber$$ The total kinetic energy $KE$ and gravitational potential energy $PE$ is: $$\begin{align} KE & = \half m (\dot x^2 + \dot y^2) \nonumber\\ & = \half mL^2\dot{\theta^2}\nonumber\\ PE & = mgy\nonumber\\ & = -mgL\cos\theta\nonumber\\ \end{align}$$ This gives us a Langrangian: $$\mathcal{L} = \half mL^2\dot{\theta^2} + mgL\cos\theta\nonumber$$

The motion is totally described by the angle $\theta$, so to calculate the motion, we then solve the differential equation $$\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot\theta} - \frac{\partial \mathcal{L}}{\partial\theta} = 0\nonumber$$ to get: $$\ddot\theta = -\omega_0^2\sin\theta\label{eq13}$$ where $\omega_0 \equiv g/L$ is the natural resonant frequency of oscillation.

The constants of the motion are: $\theta(0) = \theta_0$ and $\dot\theta(0) = 0$.

This equation cannot be solved analytically, so we will have to employ numerical techniques to be able to plot the motion.

In the section above on projectile motion, we considered 2nd order differential equations (aka Newton's laws) where the 2nd derivative of the position was either constant, or a function of the velocity. The 2nd order equation $\ddot y = f(\dot y,t)$ can be turned into a first order equation and integrated numerically via the substitution $z\equiv \dot y$, giving $$\dot z = f(z,t)\nonumber$$ This is essentially what we did above.

Now we consider 2nd order equations like the one for the pendulum (or harmonic oscillator) by making the following definitions: $$y\equiv \dot \theta\nonumber$$ $$x\equiv \theta\nonumber$$ This changes the 2nd order equation (\ref{eq13}) into two first order equations: $$\begin{align} \dot y & = f(x)\nonumber\\ & = -\omega_0^2\sin(x)\label{eq14}\\ \dot x & = g(y)\nonumber\\ & = y\label{eq15}\end{align}$$ We can then apply the Euler, Euler-Cromer, and RK2 technique directly.

Starting with Euler, we can integrate the equation for the new variables $x$ and $y$ to get: $$\begin{align} y(t+\dt)& =y(t)+\dt\cdot f(x(t))\nonumber\\ & = y(t)-\dt\cdot\omega_0^2\sin(x(t))\label{eq16}\\ x(t+\dt)& =x(t)+\dt\cdot g(y(t))\nonumber\\ & = x(t)+\dt\cdot y(t)\label{eq17} \end{align}$$ Note that we leave the time dependence explicitly to differentiate using (e.g.) $y(t)$ from $y(t+\dt)$.

Going back from $x,y$ to $\theta,\dot\theta$ gives the following 2 equations for the Euler approximation: $$\begin{align} \dot\theta(t+\dt)&=\dot\theta(t)-\dt\cdot\omega_0^2\sin(\theta)\nonumber\\ \theta(t+\dt)&=\theta(t)+\dt\cdot\dot\theta(t)\nonumber\end{align}$$ The Euler-Cromer approximation is then trivial to write down. We would want to use the value of the slope at the end of the interval, $\dot\theta(t+\dt)$, instead of the value at the beginning, in order to swim $\theta$ from $\theta(t)$ to $\theta(t+\dt)$: $$\dot\theta(t+\dt)=\dot\theta(t)-\dt\cdot\omega_0^2\sin(\theta)\nonumber\\$$ $$\theta(t+\dt)=\theta(t)+\dt\cdot\dot\theta(t+\dt)\nonumber$$ The 2nd order Runge-Kutta (RK2) method is also straightforward, but it's worth being explicit. We start with the definitions of equations (\ref{eq14}) and (\ref{eq15}), and integrate equations (\ref{eq16}) and (\ref{eq17}), and employ the Runge-Kutta suggestion to use the slope at the midpoint for the swimming: $$y(t+\dt)=y(t)+\dt\cdot f(x(t+\half \dt))\label{eq18}$$ $$x(t+\dt)=x(t)+\dt\cdot g(y(t+\half \dt))\label{eq19}$$ where we use the Euler method to calculate $y(t+\half \dt)$ and $x(t+\half \dt)$: $$y(t+\half \dt)=y(t)+\half \dt\cdot f(x(t))\nonumber\\$$ $$x(t+\half \dt)=x(t)+\half \dt\cdot g(y(t))\nonumber$$ Converting back from $x,y$ to $\theta,\dot\theta$, substituting for $f(x)$ and $y(x)$ in equations (\ref{eq14}) and (\ref{eq15}), and combining the above 4 equations give: $$\dot\theta(t+\dt)=\dot\theta(t)-\dt\cdot\omega_0^2 \sin(\theta(t)+\half \dt\dot\theta(t))\nonumber$$ $$\theta(t+\dt)=\theta(t)+\dt\cdot\dot\theta(t)-\half \dt^2 \omega_0^2\sin(\theta(t))\nonumber$$ Here you can see clearly the accuracy of the RK2 method comes from the fact that it contains terms of order $\dt^2$, as pointed out above. This is basically just saying that when you do a Taylor expansion of the integral you are trying to calculate numerically, with RK2 you are keeping higher order terms.

A comment about 4th order Runge-Kutta (aka RK4): this involves breaking up the interval into 4 segments instead of 2 with RK2. It is definitely more accurate, however whether it's worthwhile depends on the problem being tackled. One can easily get a similar accuracy with RK2 by just making the interval $\dt$ smaller, approaching the same accuracy for RK4. However if the process you are simulating is taking lots of CPU time, then RK4 might get you there faster.

In the simulation below, we use blue for the Euler method, yellow for the Euler-Cromer, and red for the RK2 method. The first box below the pendulum shows the energy (somewhat arbitrary units) for the 3 different methodds. The Euler method violates energy conservation maximally and diverges, the Euler-Cromer violates it but not so bad, and it oscillates. The RK2 method is quite good, you can't even see the energy violation without a large scale factor like 100 (try it!).


Euler Euler-Cromer Runge-Kutta (RK2)
$\dt$ 0.25
Time between calls 5

Scale by  1

Energy (arbitrary units)

To make the RK2 coding easy (or easier), it is probably best to structure things according to what's in equations (\ref{eq18}) and (\ref{eq19}), where you first find $\dot\theta(t+\half \dt)$ and use those as arguments $\theta(t+\half \dt)$ to the functions you define using equations (\ref{eq14}) and (\ref{eq15}). The code might look something like this, where dt =$\dt$, theta =$\theta$, thetaDot =$\dot\theta$, f1 and f1 are as defined in equations (\ref{eq14}) and (\ref{eq15}) respectively(and below, omega2=$\omega_0^2)$.

	theta_half = phi + 0.5*dt*g1(theta,thetaDot);
	thetaDot_half = thetaDot + 0.5*dt*f1(theta,thetaDot);
	thetaDot_new = thetaDot + dt*f1(theta_half,thetaDot_half);
	theta_new = theta + dt*g1(theta_half,thetaDot_half);
Then the code for the functions would look something like this:
function f1(x,y) {
	//
	// x is the variable, y is the 1st derivative
	//
	return -omega2*Math.sin(x);
}
function g1(x,y) {
	//
	// x is the variable, y is the 1st derivative
	//
	return y;
}
Back to top


Copyright © 2019 by Drew Baden

All rights reserved. No part of this publication may be reproduced, distributed, or transmitted in any form or by any means, including photocopying, recording, or other electronic or mechanical methods, without the prior written permission of the publisher, except in the case of brief quotations embodied in critical reviews and certain other noncommercial uses permitted by copyright law.