3.4. Implicit Runge-Kutta methods exercises#

Exercise 3.1

Determine the order of the DIRK method shown below.

\[\begin{split} \begin{array}{c|cc} \frac{1}{4} & \frac{1}{4} \\ \frac{3}{4} & \frac{1}{2} & \frac{1}{4} \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{split}\]
Solution

This DIRK method is a second-order method.

Exercise 3.2

Derive a third-order Radau IIA method. Present your method in a Butcher tableau.

Solution
\[\begin{split} \begin{array}{c|cc} \frac{1}{3} & \frac{5}{12} & -\frac{1}{12} \\ 1 & \frac{3}{4} & \frac{1}{4} \\ \hline & \frac{3}{4} & \frac{1}{4} \end{array} \end{split}\]

Exercise 3.3

The Robertson [1966] model describes the following chemical reaction between three chemicals

\[\begin{split} \begin{align*} A &\xrightarrow{k_1} B, \\ B &\xrightarrow{k_2} C, \\ B + B &\xrightarrow{k_3} C, \end{align*} \end{split}\]

where \(A\), \(B\) and \(C\) are the concentration of the three chemicals and \(k_1 = 0.04\), \(k_2 = 3 \times 10^{7}\) and \(k_3 = 10^{4}\) are parameters determining the reaction rates. Writing these as a system of ODEs gives

\[\begin{split} \begin{align*} \frac{\mathrm{d}A}{\mathrm{d}t} &= -k_1A, \\ \frac{\mathrm{d}B}{\mathrm{d}t} &= k_1A - k_2B - k_3B^2, \\ \frac{\mathrm{d}C}{\mathrm{d}t} &= k_2B + k_3B^2. \end{align*} \end{split}\]

(a)   Write a Python or MATLAB program to solve the Robertson model over \(t \in [0, 100]\), with initial values \(A(0) = 1\) and \(B(0) = C(0) = 0\) using the Radau IIA model from Exercise with a step length of \(h = 1\). Produce a plot of the concentrations of the three chemicals over \(t\) on the same axes, with \(B\) multiplied by a scaling factor of \(10^5\) so that it is visible on the plot.

Solution
../_images/282a0ae1efe964e96553da738e5fe66f760c7468099e85c7f51c75533782a4d1.png

(b)   Solve the same IVP using the Fehlberg 5(4) ERK method. Determine the smallest step length used throughout the solution.

Solution

The smallest step length used by the Fehlberg 5(4) method was \(h = 1.48 \times 10^{-4}\)

(c)   Comment on the applicability of the two methods used here to solve this IVP. Why do your think one method took longer than the other?

Solution

The Fehlberg 5(4) ERK method took much longer than the Radau IIA IRK method to compute the solution. This was because the chemical reaction is violent as indicated by the plot of the solution which shows the concentration of \(B\) increasing rapidly shortly after the start of the reaction. This means the Fehlberg 5(4) ERK method is having to use a very small step length to ensure a stable solution. The Robertson model is an example of a stiff problem which the Radau class of IRK methods are more suited to solving than an ERK method (stiffness is covered in the next chapter).