Louis Moresi

Professor of Geophysics University of Melbourne Victoria 3010 Australia +61 4 0333 1413

Computational Geodynamics: Decay Equation


Sections

We start with the numerical solution of a very simple differential equation. In fact we choose something simple enough that we already know the answer.

$$\frac{d\theta}{dt} = - k \theta \label{eq:ode_decay}$$

This is the equation which governs radioactive decay, in which case $$\theta$$ is the amount of the radioactive isotope remaining and $$d\theta / dt$$ is the activity that we can measure. $$k$$ is closely related to the half life.

The solution to this equation is $\theta(t) = \theta_0 e^{-kt}$ where $$\theta_0$$ is the amount of the radioactive material remaining. The same equation also describes the cooling of, say, a cup of coffee. In this case we interpret $$\theta$$ as the excess temperature (above room temperature).

We want to march forward in time from our starting point where $$\theta = \theta_0$$ to obtain the value of $$\theta$$ at later times. To do this, we need to approximate the original differential equation, and, in particular, the value of the time derivative at each time. There are a number of ways to do this.

A first order numerical approximation

Assume that the variation in $$\theta(t)$$ is linear, i.e. $\theta(t’) = \theta _ n + \beta t’$ where we use a local time coordinate $$t’ = t - n\Delta t$$, so that when we differentiate $\frac{d \theta}{dt} = \beta$ To determine the approximation for the derivative therefore becomes the solution to the following equation: $\begin{split} & \theta _ {n+1} = \theta _ n + \beta \Delta t & \Rightarrow \beta = \frac{d \theta}{dt} = \frac{\theta _ {n+1} - \theta _ n}{\Delta t} \end{split}$

This is a first order difference expression for the derivative which we substitute into the original differential equation ($$\ref{eq:ode_decay}$$ ) at the current timestep $\frac{\theta _ {n+1} - \theta _ n}{\Delta t} = - k \theta _ n$

This rearranges to give us a time-marching algorithm:

$\theta _ {n+1} = \theta _ n (1-k \Delta t)$

It is an indication of the fact that this problem is really not all that difficult that this difference equation can be written recursively to give: $\theta_{n+1} = \theta_0 (1-k \Delta t)^n$ In a moment we will compute some values for this expression to see how accurate it is. First we consider whether we can improve the accuracy of the approximation by doing a bit more work.

Higher order expansion

First we try fitting the local expansion for $$\theta$$ through an additional point. This time we assume that the variation in $$\theta(t)$$ is quadratic, i.e. $\theta(t’) = \theta _ {n-1} + \beta t’ + \gamma {t’}^2$ The local time coordinate is $$t’ = t - (n-1)\Delta t$$, and when we differentiate $$\frac{d \theta}{dt} = \beta + 2 \gamma t’ \label{eq:dthdt2}$$

To solve for $$\beta$$ and $$\gamma$$ we fit the curve through the sample points: $\begin{split} \theta _ n &= \theta _ {n-1} + \beta \Delta t + \gamma (\Delta t)^2 \theta _ {n+1} &= \theta _ {n-1} + 2 \beta \Delta t + 4 \gamma (\Delta t)^2 \end{split}$ Which we can solve to give this: $\begin{split} \beta &= \left( 4 \theta _ n - \theta_{n+1} - 3\theta _ {n-1} \right) \frac{1}{2\Delta t} \gamma &= \left( \theta _ {n+1} + \theta _ {n-1} -2 \theta _ n \right) \frac{1}{2\Delta t^2} \end{split}$ By substituting into $$(\ref{eq:dthdt2})$$, and then into the original differential equation $$(\ref{eq:ode_decay})$$ we obtain the following $\left. \frac{d\theta}{dt} \right| _ {t=n\Delta t} = \beta + 2\gamma \Delta t = \frac{1}{2\Delta t} \left( \theta _ {n+1} - \theta _ {n-1} \right) = -k \theta _ n$ The difference approximation to the derivative turns out to be the average of the expressions for the previous derivative and the new derivative. We have now included information about the current timestep and the previous timestep in our expression for the value of $$\theta$$ at the forthcoming timestep: $\theta _ {n+1} = \theta _ {n-1} -2k \theta _ n \Delta t$

A comparison of the two schemes

How well do these two methods perform, and what difference does the choice of $$\Delta t$$ make on the accuracy ?

t $$^1\theta _ H$$ $$^1\theta _ h$$ $$^2\theta _ H$$ $$^2\theta _ h$$ Exact
0.0 1.0000 1.0000 1.0000 1.0000 1.0000
0.1 - 0.9000 - 0.9000 0.9048
0.2 0.8000 0.8100 0.8000 0.8200 0.8187
0.3 - 0.7290 - 0.7360 0.7408
0.4 0.6400 0.6561 0.6800 0.6728 0.6703
0.5 - 0.5905 - 0.6014 0.6063
0.6 0.5120 0.5314 0.5280 0.5525 0.5488
0.7 - 0.4783 - 0.4909 0.4966
0.8 0.4096 0.4305 0.4688 0.4543 0.4483
0.9 - 0.3874 - 0.4001 0.4066
1.0 0.3277 0.3487 0.3405 0.3743 0.3679

Results of the numerical simulation of the decay equation. $$^1\theta_H$$ is from the first order expansion of $$\theta(t)$$, at a timestep of $$\Delta t=0.2$$, and $$^2\theta_H$$ is the second order expansion at $$\Delta t=0.2$$; other calculations are at $$\Delta t = 0.1$$. $$^1\theta_h$$ is the linear expansion, $$^2\theta_h$$ results from the quadratic expansion of $$\theta(t)$$. The actual solution is also shown.

The results are more accurate when a smaller timestep is used although it requires more computation to achieve the greater accuracy. Higher order expansion also increases the accuracy and may be more efficient in terms of the number of computations required for a given level of accuracy. Note, however, that the supposedly better quadratic expansion produces an error which oscillates as time increases. Does this error grow ? Does this make second order expansions useless ?

Second Order Runge-Kutta

The Runge-Kutta approach to higher order integration methods is illustrated above. The idea is to estimate the gradient $$d \theta / d t$$ at the half way point between two timestep values. This is done in two stages. Initially a first order estimate, $$\hat{\theta}$$ is made for the value of the function $$\theta$$ at $$t=t+\Delta t /2$$ in the future. This value is then subsituted into the differential equation to obtain the estimate for the gradient at this time. The revised gradient is then used to update the original $$\theta(t)$$ by an entire timestep.

The first order step is $\begin{split} \hat{\theta}(t+\Delta t /2) & = \theta(t) + \left. \frac{d \theta}{d t} \right| _ t \frac{\Delta t}{2} &= \theta(t) \left[ 1-\frac{k\Delta t}{2} \right] \end{split}$ Substitute to estimate the gradient at the mid-point $\left. \frac{d \theta}{d t} \right| _ {t+\Delta t /2} \approx -k \theta(t) \left[ 1-\frac{k\Delta t}{2} \right]$ Use this value as the average gradient over the interval $$t\rightarrow t+\Delta t$$ to update $$\theta$$ $\begin{split} \theta(t+\Delta t) & \approx \theta(t) + \delta t \left( -k \theta(t) \left[ 1-\frac{k\Delta t}{2} \right] \right) & \approx \theta(t) \left( 1 - k \Delta t + k^2 \frac{\Delta t^2}{2} \right) \end{split}$ It’s worth noting that the Taylor expansion of the solution should look like $e^{-kt} = 1 - kt + \frac{k^2 t^2}{2!} - \frac{k^3 t^3}{3!} + \ldots$

The Runge Kutta method can be extended by repeating the estimates on smaller regions of the interval. The usual choice is fourth order RK. This is largely because, obviously, it’s accurate to fourth order, but also because the number of operations to go higher than fourth order is disproportionately large. See Numerical Recipes for a discussion on this and better methods for ODE’s.

Our First Code

The simplest way to answer our earlier question is to code the methods and see. The following is a script which computes the first and second order expressions and the Runge-Kutta method, and compares with the “exact” solution (which is numerical as well, of course, but computed through series expansions and the like).

Clearly, the quadratic expansion does not do a good job in this particular case, since the error is a growing instability which makes the solution useless almost immediately. However, the second order Runge-Kutta method is very accurate. This gives a clear warning that numerical solution of equations is partly an art (though the stability or otherwise of a given method can actually be proven formally in many cases).

Explicit versus Implicit Methods}

The use of implicit solution methods is crucial when we come to the mantle flow problem (and to elasticity). It is a deceptively simple concept which we use to eliminate oscillatory errors and other problems associated with stiff systems, and in other cases where we want to take large timesteps. Again, Numerical Recipes provides a good example of the way in which implicit methods avoid the shuffling timesteps which result when an otherwise unimportant component of the system blows up. To obtain an implicit solution, we express our updated timestep in terms of values of derivatives etc which are evaluated at that timestep. This may result in dependencies which mean large systems of simultaneous equations need to be solved.

But our simple example can also be solved implicitly as follows $\theta(t+\Delta t) = \left. \frac{d \theta}{d t} \right| _ {t+\Delta t} + \theta(t)$ Easy to write down, but $$\left. d \theta/ dt \right| _ {t+\Delta t}$$ is not known until $$\theta(t+\Delta t)$$ is known. In this case, however, we simply substitute $$\left. d \theta/ dt \right| _ {t+\Delta t}= -k \theta(t+\Delta t)$$ from the differential equation and write $\begin{split} \theta(t+\Delta t) & = -k \theta(t+\Delta t) \Delta t + \theta(t) & = \frac{\theta(t) }{ 1 + k\Delta t} \end{split}$

A higher order implicit method can be found by considering how the second order Runge-Kutta method constructs its estimates. In place of the first order approximation above, we try $\theta(t+\Delta t) = \left. \frac{d \theta}{d t} \right| _ {t+\Delta t/2} + \theta(t)$ where we assume $\left. \frac{d \theta}{d t} \right| _ {t+\Delta t/2} = \frac{1}{2} \left. \frac{d \theta}{d t} \right| _ {t+\Delta t} + \frac{1}{2} \left. \frac{d \theta}{d t} \right| _ {t}$ Substituting from our original differential equation now gives $\begin{split} \theta(t+\Delta t) &= -k\frac{\theta(t+\Delta t) + \theta(t)}{2} \Delta t + \theta(t) &= \theta(t) \frac{1-k\Delta t / 2}{1+k\Delta t / 2} \end{split}$

$$t$$ $$^1\theta_i$$ $$^2\theta_i$$ Exact
0.0 1.0000 1.0000 1.0000
0.1 0.9091 0.9048 0.9048
0.2 0.8264 0.8186 0.8187
0.3 0.7513 0.7406 0.7408
0.4 0.6830 0.6701 0.6703
0.5 0.6209 0.6063 0.6063

Results of the two implicit schemes derived above for the numerical integration of the decay equation

The first order scheme is of approximately the same accuracy as the explicit scheme, not surprisingly, but the mid-point implicit method has very impressive results considering that there is very little extra work involved in deriving the formulation. Can this strategy stabilize the quadratic expansion which proved to be so unstable in explicit form ? Evaluating the derivatives at time $$t=(n+1)\Delta t$$ gives

$\left. \frac{d \theta}{d t} \right| _ {(n+1)\Delta t} = 3 \theta _ {n+1} - 4\theta _ n + \theta _ {n-1}$

which we substitute into the differential equation at $$t=(n+1)\Delta t$$ to give (eventually)

$\theta _ {n+1} = \frac{4\theta _ n - \theta _ {n-1}}{3+2k\Delta t}$

This method is stable although not as accurate as the second order Runge-Kutta Scheme. It is trivial to modify the Python script to demonstrate this.