Examination Questions#
The following are example examination questions that are similar in style and scope to the ones that you can expect to see in your examination. You are advised to work through these in preparation for the examination. The examination takes place in a PC lab where you will be permitted access to Python (Jupyter Notebook or Spyder), MATLAB and Excel. You are advised to familiarise yourself with one or more to help you answer the questions.
Explicit Runge-Kutta Methods#
(a) An explicit Runge-Kutta method is expressed using the Butcher tableau below.
\[\begin{split} \begin{align*} \begin{array}{c|cccc} 0 & 0 \\ \frac{1}{4} & \frac{1}{4} \\ \frac{1}{2} & 0 & \frac{1}{2} \\ 1 & 1 & -2 & 2 \\ \hline & \frac{1}{6} & 0 & \frac{2}{3} & \frac{1}{6} \end{array} \end{align*} \end{split}\]Use this Runge-Kutta method with a step length of \(h = 0.2\) to solve the following IVP
\[ \begin{align*} y' &= t - 2 y, & t &\in [0, 1], & y(0) &= 2. \end{align*} \][10 marks]Solution
\[\begin{split} \begin{align*} \mathbf{t} &= (0, 0.2, 0.4, 0.6, 0.8, 1), \\ \mathbf{y} &= (2, 1.3584, 0.9612, 0.7279, 0.6045, 0.5547). \end{align*} \end{split}\]import numpy as np def f(t, y): return t - 2 * y h = 0.2 nsteps = 5 t = np.arange(0, 1 + h, h) y = np.zeros(t.shape) for n in range(nsteps): k1 = f(t[n], y[n]) k2 = f(t[n] + 1/4 * h, y[n] + 1/4 * h * k1) k3 = f(t[n] + 1/2 * h, y[n] + 1/2 * h * k2) k4 = f(t[n] + h, y[n] + h * (k1 - 2 * k2 + 2 * k3)) y[n+1] = y[n] + h / 6 * (k1 + 4 * k3 + k4) print(t) print(y)
f = @(t, y) t - 2 * y; nsteps = 5; h = 0.2; t = 0 : h : 1; y = zeros(size(t)); y(1) = 2; for n = 1 : nsteps k1 = f(t(n), y(n)); k2 = f(t(n) + 1/4 * h, y(n) + 1/4 * h * k1); k3 = f(t(n) + 1/2 * h, y(n) + 1/2 * h * k2); k4 = f(t(n) + h, y(n) + h * (k1 - 2 * k2 + 2 * k3)); y(n+1) = y(n) + h / 6 * (k1 + 4 * k3 + k4); end t y
Note: solutions may also be calculated using a pen and calculator or using Excel.
(b) The exact solution to the IVP from part (a) is \(y(t) = \frac{1}{9}(3t + 19 e^{-3t} - 1)\). Calculate the global truncation error for the numerical solution to the IVP for \(y(1)\) using the Runge-Kutta method from part (a) with \(h = 0.2\) and \(h = 0.1\). Use your global truncation errors to estimate the order accuracy of this Runge-Kutta method.
[10 marks]Solution
The exact solution is \(y(1) = 0.554504\) and numerical solutions are \(y(1, h= 0.2) = 0.554686\) and \(y(1,h = 0.1) = 0.554514\) so the GTE are
\[\begin{split} \begin{align*} E(h = 0.2) = |0.554504 - 0.554686| = 1.82 \times 10^{-4}, \\ E(h = 0.1) = |0.554504 - 0.554515| = 9.60 \times 10^{-5}, \end{align*} \end{split}\]therefore
\[ \text{order} \approx \frac{\log \left( \dfrac{E(h1)}{E(h2)}\right)}{\log \left( \dfrac{h_1}{h_2}\right)} = \frac{\log \left( \dfrac{1.82 \times 10^{-4}}{9.60 \times 10^{-5}} \right)}{\log \left( \dfrac{0.2}{0.1} \right)} = 4.24. \](c) Explain in words how step size control is implemented within a Runge-Kutta method.
[5 marks]Solution
(a) The order conditions for a fourth-order explicit Runge-Kutta method are:
\[\begin{split} \begin{align*} b_1 + b_2 + b_3 + b_4 &= 1, \\ b_2c_2 + b_3c_3 + b_4c_4 &= \frac12, \\ b_2c_2^2 + b_3c_3^2 + b_4c_4^2 &= \frac13, \\ b_2c_2^3 + b_3c_3^3 + b_4c_4^3 &= \frac14, \\ b_3a_{32}c_2 + b_4a_{42}c_2 + b_4a_{43}c_3 &= \frac{1}{6}, \\ b_3a_{32}c_2c_3 + b_4a_{42}c_2c_4 + b_4a_{43}c_3c_4 &= \frac18, \\ b_3a_{32}c_2^2 + b_4a_{42}c_2^2 + b_4a_{43}c_3^2 &= \frac{1}{12}, \\ b_4a_{43}a_{32}c_2 &= \frac{1}{24}, \\ c_2 &= a_{21}, \\ c_3 &= a_{31} + a_{32}, \\ c_4 &= a_{41} + a_{42} + a_{43}. \end{align*} \end{split}\]Derive a fourth-order explicit Runge-Kutta method with \(c_2 = \frac13\), \(c_3 = \frac23\) and \(c_4=1\). Write down your method in the form of a Butcher tableau.
[10 marks]Solution
\[\begin{split} \begin{array}{c|cccc} 0 & 0 \\ \frac13 & \frac13 \\ \frac23 & -\frac13 & 1 \\ 1 & 1 & -1 & 1 \\ \hline & \frac18 & \frac38 & \frac38 & \frac18 \end{array} \end{split}\]import sympy as sp b1, b2, b3, b4, a21, a31, a32, a41, a42, a43 = sp.symbols("b1, b2, b3, b4, a21, a31, a32, a41, a42, a43") c2 = sp.Rational(1,3) c3 = sp.Rational(2,3) c4 = 1 eq1 = b1 + b2 + b3 + b4 - 1 eq2 = b2 * c2 + b3 * c3 + b4 * c4 - sp.Rational(1,2) eq3 = b2 * c2**2 + b3 * c3**2 + b4 * c4**2 - sp.Rational(1,3) eq4 = b2 * c2**3 + b3 * c3**3 + b4 * c4**3 - sp.Rational(1,4) eq5 = b3 * a32 * c2 + b4 * a42 * c2 + b4 * a43 * c3 - sp.Rational(1,6) eq6 = b3 * a32 * c2 * c3 + b4 * a42 * c2 * c4 + b4 * a43 * c3 * c4 - sp.Rational(1,8) eq7 = b3 * a32 * c2**2 + b4 * a42 * c2**2 + b4 * a43 * c3**2 - sp.Rational(1,12) eq8 = b4 * a43 * a32 * c2 - sp.Rational(1,24) eq9 = c2 - a21 eq10 = c3 - a31 - a32 eq11 = c4 - a41 - a42 - a43 sp.solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11))
(b) Use your method from part (a) to solve the following initial value problem using a step length of \(h = 0.1\).
\[ y' = -3ty, \qquad t \in [0, \tfrac12], \qquad y(0) = 1. \]Present the solution as a table of \(t\) and \(y\) values correct to 6 decimal places.
[10 marks]Solution
\(t\)
\(y\)
0.0
1.000000
0.1
0.985112
0.2
0.941764
0.3
0.873714
0.4
0.786625
0.5
0.687286
import numpy as np def rk4(f, tspan, y0, h): nsteps = int((tspan[1] - tspan[0]) / h) t = np.arange(0, tspan[1] + h, h) y = np.zeros(t.shape) y[0] = y0 for n in range(nsteps): k1 = f(t[n], y[n]) k2 = f(t[n] + 1/3 * h, y[n] + 1/3 * h * k1) k3 = f(t[n] + 2/3 * h, y[n] + h * (-1/3 * k1 + k2)) k4 = f(t[n] + h, y[n] + h * (k1 - k2 + k3)) y[n+1] = y[n] + h / 8 * (k1 + 3 * k2 + 3 * k3 + k4) return t, y def f(t, y): return -3 * t * y h = 0.1 y0 = 1 tspan = [0, 1/2] t, y = rk4(f, tspan, y0, h) print(f"t = {t}") print(f"y = {y}")
(c) The exact solution of the IVP from part (b) is \(y = e^{-(3t^2)/2}\). Compute the global truncation errors for your Runge-Kutta method for \(y(\frac12)\) with step lengths \(h=0.1\) and \(h = 0.05\). Use the global truncation errors to estimate the order of accuracy of your Runge-Kutta method from part (a).
[5 marks]Solution
The numerical solutions using the two step lengths are
\[\begin{split} \begin{align*} y_5(h = 0.1) &= 3.47 \times 10^{-6}, \\ y_5(h = 0.05) &= 1.97 \times 10^{-7}, \end{align*} \end{split}\]so the order estimate is
\[ \text{order} = \frac{\log \left( \frac{3.47 \times 10^{-6}}{1.97 \times 10^{-7}} \right)}{ \log \left( \frac{0.1}{0.05} \right)} = 4.14.\](a) A Runge-Kutta method is expressed in the following Butcher tableau
\[\begin{split} \begin{array}{c|cccc} 0 & 0 \\ \frac13 & \frac13 \\ \frac12 & \frac58 & -\frac18 \\ 1 & \frac32 & \frac32 & -2 \\ \hline & -\frac12 & 3 & -2 & \frac12 \end{array} \end{split}\]Determin the stability function of this Runge-Kutta method, and hence state the order of the method.
The stabilty function of an explicit Runge-Kutta method is
\[ R(z) = 1 + z \sum_{k=0}^{s-1} \mathbf{b}^\mathsf{T}A^k \mathbf{e} z^{k+1}. \]Computing the coefficents
\[\begin{split} \begin{align*} k &= 0: & \mathbf{b}^\mathsf{T}A^0 \mathbf{e} &= (-\tfrac12, 3, -2, \tfrac12) \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = 1, \\ k &= 1: & \mathbf{b}^\mathsf{T}A \mathbf{e} &= (-\tfrac12, 3, -2, \tfrac12) \begin{pmatrix} 0 & 0 & 0 & 0 \\ \frac13 & 0 & 0 & 0 \\ \frac58 & -\frac18 & 0 & 0 \\ \frac32 & \frac32 & -2 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \frac12, \\ k &= 2: & \mathbf{b}^\mathsf{T}A \mathbf{e} &= (-\tfrac12, 3, -2, \tfrac12) \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -\frac{1}{24} & 0 & 0 & 0 \\ -\frac34 & \frac14 & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = -\frac16, \\ k &= 3: & \mathbf{b}^\mathsf{T}A \mathbf{e} &= (-\tfrac12, 3, -2, \tfrac12) \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \frac{1}{12} & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} = \frac{1}{24}, \end{align*} \end{split}\]therefore the stability function is
\[ R(z) = 1 + z + \frac12 z^2 - \frac16 z^3 + \frac{1}{24} z^4. \]\(R(z)\) is the series expansion of \(e^z\) correct to the second-order term, so this Runge-Kutta method is second-order accurate.
(b) An implicit Runge-Kutta method is expressed as the following Butcher tableau
\[\begin{split} \begin{array}{c|cc} \frac14 & \frac14 & 0 \\ \frac34 & \frac12 & \frac14 \\ \hline & \frac12 & \frac12 \end{array} \end{split}\]Determine the stability function of this method and state with reasons whether it is A-stable or not.
Solution
The stability function of an implicit Runge-Kutta method is
\[ R(z) = \frac{\det(I - zA + z \mathbf{eb}^\mathsf{T})}{\det(I - zA)}. \]Therefore
\[\begin{split} \begin{align*} R(z) &= \frac{ \det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac14 & 0 \\ \frac12 & \frac14 \end{pmatrix} + z \begin{pmatrix} \frac12 & \frac12 \\ \frac12 & \frac12 \end{pmatrix} \right) }{ \det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac14 & 0 \\ \frac12 & \frac14 \end{pmatrix} \right) } \\ &= \frac{ \det \begin{pmatrix} 1 + \frac14 z & \frac12 z \\ 0 & 1 + \frac14 z \end{pmatrix} }{ \det \begin{pmatrix} 1 - \frac14 z & 0 \\ -\frac12 z & 1 - \frac14 z \end{pmatrix} } = \frac{1 + \frac12 z + \frac{1}{16}z^2}{1 - \frac12 z + \frac{1}{16}z^2}. \end{align*} \end{split}\]An IRK method is A-stable if
the roots of \(Q(z)\) have positive real parts and
the polynomial \(E(y) = |Q(iy)|^2 - |P(iy)|^2\) is greater than or equal to zero for all \(y \in \mathbb{R}\)
where \(R(z) = P(z) / Q(z)\). The roots of \(Q(z)\) is \(z = \frac14\) so the first condition is satisfied. The polynomial \(E(y)\) is
\[\begin{split} \begin{align*} E(y) &= \left| 1 - \frac12 iy - \frac{1}{16} y^2 \right|^2 - \left| 1 + \frac12 iy - \frac{1}{16} y^2 \right|^2 \\ &= \left( 1 - \frac{1}{16} y^2 \right)^2 + \left( -\frac12 y \right)^2 - \left( 1 - \frac{1}{16} y^2 \right)^2 - \left( \frac12 y \right)^2 \\ &= \frac14 y^2 - \frac14 y^2 = 0. \end{align*} \end{split}\]So the second condition is satisfied and this is an A-stable method.
(c) The initial value problem below
\[ y' = t - 4y, \qquad t \in [0, 1], \qquad y(0) = 1, \]is to be solved using the Runge-Kutta method from part (b) using a step length of \(h = 0.1\).
(i) The stage values of an implicit Runge-Kutta method can be calculated using the following
\[ Y_i = y_n + h \sum_{j=1}^s a_{ij} f(t_n + c_jh, Y_j). \]Calculate the stage values \(Y_1\) and \(Y_2\) for the first step.
Solution
\[\begin{split} \begin{align*} Y_1 &= y_n + h (\tfrac14 (t_n + \tfrac14 h - 4Y_1)), \\ Y_2 &= y_n + h (\tfrac12 (t_n + \tfrac14 h - 4Y_1) + \tfrac14(t_n + \tfrac34 h - 4Y_2)), \end{align*} \end{split}\]which can be rearranged to
\[\begin{split} \begin{align*} (1 + h)Y_1 &= y_n + \tfrac14 ht_n + \tfrac{1}{16}h^2, \\ 2hY_1 + (1 + h) Y_2 &= y_n + \tfrac34 ht_n + \tfrac{5}{16}h^2, \end{align*} \end{split}\]and writing as a matrix equation
\[\begin{split} \begin{align*} \begin{pmatrix} 1 + h & 0 \\ 2h & 1 + h \end{pmatrix} \begin{pmatrix} Y_1 \\ Y_2 \end{pmatrix} = \begin{pmatrix} y_n + \frac14 ht_n + \frac{1}{16}h^2 \\ y_n + \frac34 ht_n + \frac{5}{16} h^2 \end{pmatrix}. \end{align*} \end{split}\]For the first step, \(t_0 = 0\), \(y_0 = 1\) and since \(h = 0.1\) then
\[\begin{split} \begin{align*} \begin{pmatrix} 1.1 & 0 \\ 0.2 & 1.1 \end{pmatrix} \begin{pmatrix} Y_1 \\ Y_2 \end{pmatrix} &= \begin{pmatrix} 1.000625 \\ 1.003125 \end{pmatrix} \\ \begin{pmatrix} Y_1 \\ Y_2 \end{pmatrix} &= \begin{pmatrix} \frac{11}{10} & 0 \\ -\frac{20}{121} & \frac{11}{10} \end{pmatrix} \begin{pmatrix} 1.000625 \\ 1.003125 \end{pmatrix} = \begin{pmatrix} 0.909659 \\ 0.746539 \end{pmatrix}. \end{align*} \end{split}\]Therefore \(Y_1 = 0.909659\) and \(Y_2 = 0.746539\).
(ii) Hence calculate the solution at \(t = 0.1\).
Solution
\[\begin{split} \begin{align*} y_1 &= y_0 + \frac{h}{2} (f(t_0 + \tfrac14 h, Y_1) + f(t_0 + \tfrac34 h, Y_2)) \\ &= 1 + \frac{0.1}{2}((0 + \tfrac14(0.1) - 4(0.909659)) + (0 + \tfrac34 (0.1) - 4(0.746539))) \\ &= 0.673760. \end{align*} \end{split}\](iii) Using MATLAB or Python, compute the solution to this initial value problem over the interval \(t \in [0, 1]\). Present your solution as a table with values stated using 4 decimal places.
Solution
\(n\)
\(t_n\)
\(y_n\)
0
0.0
1.0000
1
0.1
0.6738
2
0.2
0.4636
3
0.3
0.3312
4
0.4
0.2509
5
0.5
0.2053
6
0.6
0.1831
7
0.7
0.1765
8
0.8
0.1803
9
0.9
0.1912
10
1.0
0.2067
import numpy as np def f(t, y): return t - 4 * y y0 = 1 h = 0.1 nsteps = int((1 - 0) / h) t = np.arange(0, 1 + h, h) y = np.zeros(t.shape) y[0] = y0 A = np.array([[ 1 + h, 0 ], [ 2 * h, 1 + h ]]) Ainv = np.linalg.inv(A) print("| n | tn | yn |") print("|----|-----|--------|") print(f"| {0:2d} | {t[0]:0.1f} | {y[0]:0.4f} |") for n in range(nsteps): b = np.array([y[n] + 1/4 * h * t[n] + 1/16 * h**2, y[n] + 3/4 * h * t[n] + 5/16 * h**2]) Y = Ainv @ b y[n+1] = y[n] + h / 2 * (f(t[n] + 1/4 * h, Y[0]) + f(t[n] + 3/4 * h, Y[1])) t[n+1] = t[n] + h print(f"| {n+1:2d} | {t[n+1]:0.1f} | {y[n+1]:0.4f} |")
(a) The \(B(k)\) and \(C(k)\) order conditions for a Runge-Kutta method are
\[\begin{split} \begin{align*} B(k):&& \sum_{i=1}^s b_i c_i^{j-1} &= \frac{1}{j}, & j &= 1 \ldots k, \\ C(k):&& \sum_{j=1}^s a_{ij}c_j^{\ell - 1} &= \frac{1}{\ell} c_i^\ell, & i &= 1 \ldots s, & \ell &= 1 \ldots k, \end{align*} \end{split}\]and the \(n\)-order Legendre polynomial is
\[ \begin{align*} P_n(x) = \sum_{i=0}^n {n \choose i} {n + i \choose i} (x - 1)^i, \end{align*} \]where \({n \choose i}\) is the Binomial coefficient.
An \(s\)-stage Radau IIA method has order \(2s - 1\), the \(c_i\) values are the roots of \(P_s(x) - P_{s-1}(x)\) and the \(a_{ij}\) and \(b_i\) coefficients are chosen to satisfy the \(B(2s)\) and \(C(s)\) order conditions respectively.
Derive a third-order Radau IIA method expressing your solution as a Butcher tableau.
[10 marks]Solution
A third-order Radau IIA method has \(s = 2\) stages. The values of \(c_1\) and \(c_2\) are the roots of \(P_2(x) - P_1(x) = 0\)
\[\begin{split} \begin{align*} 0 &= P_2(x) - P_1(x) \\ &= 1 + {2 \choose 1}{3 \choose 1}(x - 1) + {2 \choose 2}{4 \choose 2}(x - 1)^2 - 1 - {1 \choose 1}{2 \choose 1} (x - 1) \\ &= 6x^2 - 8x + 2 \\ &= (3x - 1)(x - 1) \end{align*} \end{split}\]so \(c_1 = \frac{1}{3}\) and \(c_2 = 1\). The \(b_1\) and \(b_2\) values satisfy the \(B(2)\) order conditions
\[\begin{split} \begin{align*} b_1 + b_2 &= 1, \\ \frac{1}{3}b_1 + b_2 &= \frac{1}{2}, \end{align*} \end{split}\]so \(b_1 = \frac{3}{4}\) and \(b_2 = \frac{1}{4}\). The \(A\) values satisfy the \(C(1)\) order conditions
\[\begin{split} \begin{align*} a_{11} + a_{12} &= \frac{1}{3}, \\ a_{21} + a_{22} &= 1, \\ \frac{1}{3}a_{11} + a_{12} &= \frac{1}{18}, \\ \frac{1}{3}a_{21} + a_{22} &= \frac{1}{2}. \end{align*} \end{split}\]Subtracting the third equation from the first gives \(a_{11} = \frac{5}{12}\) so \(a_{12} = -\frac{1}{12}\). Subtracting the fourth equation from the second gives \(a_{21} = \frac{3}{4}\) and \(a_{22} = \frac{1}{4}\). So the Butcher tableau for the third-order Radau IIA method is
\[\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}\](b) A 2-stage SDIRK method is given below
\[\begin{split} \begin{align*} \begin{array}{c|cc} \frac{1}{8} & \frac{1}{8} & 0 \\ \frac{7}{8} & \frac{3}{4} & \frac{1}{8} \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{align*} \end{split}\]Write the stability function of this method and determine whether it is A-stable or not. You may find the following formulae useful:
\[\begin{split} \begin{align*} R(z) &= \frac{\det(I - zA + z\mathbf{eb}^\mathsf{T})}{\det(I - zA)}, \\ E(y) &= |Q(iy)|^2 - |P(iy)|^2. \end{align*} \end{split}\][10 marks]Solution
The stability function for this DIRK method is
\[\begin{split} \begin{align*} R(z) &= \frac{\det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac18 & 0 \\ \frac34 & \frac18 \end{pmatrix} + z \begin{pmatrix} \frac12 & \frac12 \\ \frac12 & \frac12 \end{pmatrix} \right) }{\det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac18 & 0 \\ \frac34 & \frac18 \end{pmatrix} \right) } \\ &= \frac{\det \begin{pmatrix} 1 + \frac{3}{8}z & \frac{1}{2}z \\ -\frac14 z & 1 + \frac{3}{8}z \end{pmatrix} }{\det \begin{pmatrix} 1 - \frac18 z & 0 \\ -\frac34 z & 1 - \frac18 z \end{pmatrix} } = \frac{1 + \frac34 z + \frac{17}{64}z^2}{1 - \frac14 z + \frac{1}{64}z^2} \end{align*} \end{split}\]The roots of the denominator are
\[ \begin{align*} z &= \frac{\frac14 \pm \sqrt{\frac{1}{16} - \frac{1}{16}}}{\frac{1}{32}} = 8, \end{align*} \]which is positive so the first condition for A-stability is satisfied.
\[\begin{split} \begin{align*} E(y) &= |Q(iy)|^2 - |P(iy)|^2 \\ &= \left| 1 - \frac14 iy - \frac{1}{64}y^2 \right|^2 - \left| 1 + \frac34 iy - \frac{17}{64}y^2 \right|^2 \\ &= \left( 1 - \frac{1}{64}y^2 \right)^2 + \left( -\frac{1}{4}y \right)^2 - \left( 1 - \frac{17}{64}y^2 \right)^2 - \left( \frac{3}{4}y \right)^2 \\ &= 1 - \frac{1}{32}y^2 + \frac{1}{4096}y^4 + \frac{1}{16}y^2 - 1 + \frac{17}{32}y^2 - \frac{289}{4096}y^4 - \frac{9}{16}y^2 \\ &= - \frac{9}{128}y^4 \end{align*} \end{split}\]Since \(E(y) < 0\) for all \(y \in \mathbb{R}\) so this DIRK method is non A-stable.
(c) Explain what is meant by a stiff system of equations. What methods would be most suitable for solving stiff and non-stiff problems?
[5 marks]Solution
(a) The \(B(k)\) and \(C(k)\) order conditions for an \(s\)-stage Runge-Kutta method are
\[\begin{split} \begin{align*} B(k): && \sum_{i=1}^s b_ic_i^{j - 1} &= \frac{1}{j}, & j &= 1, \ldots, k, \\ C(k): && \sum_{j=1}^s a_{ij}c_j^{\ell - 1} &= \frac{1}{\ell} c_i^\ell, & i &= 1, \ldots, s, & \ell &= 1, \ldots, k, \end{align*} \end{split}\]and the Legendre polynomial is
\[ P_n(x) = \sum_{i=0}^n {n \choose i}{n+i \choose i} (x - 1)^i, \]where \({n \choose i}\) is the Binomial coefficient.
Given that an \(s\)-stage Gauss-Legendre implicit Runge-Kutta method has order \(2s\) where the \(c_i\) coefficients are the roots of the Legendre polynomial \(P_s(x)\) and the \(b_i\) and \(a_{ij}\) coefficients satisfy the \(B(s)\) and \(C(s)\) order conditions, derive a fourth-order Gauss-Legendre method.
[10 marks]Solution
A fourth-order Gauss-Legendre method has \(s = 2\) stages. The \(c_i\) values are the roots of \(P_2(x)\)
\[\begin{split} \begin{align*} P_3(x) &= 1 + {2 \choose 1}{3 \choose 1}(x - 1) + {2 \choose 2}{4 \choose 2}(x - 1)^2 \\ &= 1 + 6(x - 1) + 6 (x^2 - 2x + 1) \\ &= 6x^2 - 6x + 1, \end{align*} \end{split}\]using Python or MATLAB to compute the roots of \(P_3(x)\) results in
\[ c_1 = \frac12 - \frac{\sqrt{3}}{6}, \qquad c_2 = \frac12 + \frac{\sqrt{3}}{6}.\]The \(b_i\) and \(a_{ij}\) satisfy the \(B(2)\) and \(c(2)\) order conditions
\[\begin{split} \begin{align*} b_1 + b_2 &= 1, \\ b_1c_1 + b_2c_2 &= \frac12, \\ a_{11} + a_{12} &= c_1, \\ a_{21} + a_{22} &= c_2, \\ a_{11}c_1 + a_{12}c_2 &= \frac12 c_1^2, \\ a_{21}c_1 + a_{22}c_2 &= \frac12 c_2^2. \end{align*} \end{split}\]Solving using Python or MATLAB gives
\[\begin{split} \begin{array}{c|cc} \frac12 - \frac{\sqrt{3}}{6} & \frac14 & \frac14 - \frac{\sqrt{3}}{6} \\ \frac12 + \frac{\sqrt{3}}{6} & \frac14 + \frac{\sqrt{3}}{6} & \frac14 \\ \hline & \frac12 & \frac12 \end{array} \end{split}\]x = sp.symbols('x') c = sp.solve(6 * x**2 - 6 * x + 1) c1 = c[0] c2 = c[1] print(c1, c2) b1, b2, a11, a12, a21, a22 = sp.symbols("b1, b2, a11, a12, a21, a22") eq1 = b1 + b2 - 1 eq2 = b1 * c1 + b2 * c2 - sp.Rational(1,2) eq3 = a11 + a12 - c1 eq4 = a21 + a22 - c2 eq5 = a11 * c1 + a12 * c2 - c1**2 / 2 eq6 = a21 * c1 + a22 * c2 - c2**2 / 2 sp.solve((eq1, eq2, eq3, eq4, eq5, eq6))
(b) Another implicit method has the following Butcher tableau:
\[\begin{split} \begin{array}{c|cc} \frac12 & \frac12 & 0 \\ \frac32 & -\frac12 & 2 \\ \hline & -\frac12 & \frac32 \end{array} \end{split}\]State with reasons whether it is A-stable or not. You may find the following useful:
\[\begin{split} \begin{align*} R(z) &= \frac{\det(I - zA + z\mathbf{eb}^\mathsf{T})}{\det(I - zA)}, \\ E(y) &= |Q(iy)|^2 - |P(iy)|^2, \end{align*} \end{split}\]where \(R(z) = P(z) / Q(z)\) and \(i^2 = -1\).
[10 marks]Solution
The stability function is
\[\begin{split} \begin{align*} R(z) &= \frac{ \det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac12 & 0 \\ -\frac12 & 2 \end{pmatrix} + z \begin{pmatrix} -\frac12 & \frac32 \\ -\frac12 & \frac32 \end{pmatrix} \right) }{ \det \left( \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - z \begin{pmatrix} \frac12 & 0 \\ -\frac12 & 2 \end{pmatrix} \right) } \\ &= \frac{\det \begin{pmatrix} 1 - z & \frac32 z \\ 0 & 1 - \frac12 z\end{pmatrix} }{ \det \begin{pmatrix} 1 - \frac12 z & 0 \\ \frac12 z & 1 - 2z \end{pmatrix} } = \frac{1 - \frac32 z + \frac12 z^2}{1 - \frac52 z + z^2}. \end{align*} \end{split}\]A method is A-stable if
the roots of \(Q(z)\) have positive real parts and
\(E(y) \geq 0\) for all \(y \in \mathbb{R}\).
Checking the roots of \(Q(z)\)
\[ \begin{align*} z &= \frac{\frac52 \pm \sqrt{\frac{25}{4} - 4}}{2} = \frac54 \pm \frac34=\frac12, 2 \end{align*} \]So the first condition is satisfied. Checking \(E(y) \ge 0\)
\[\begin{split} \begin{align*} E(y) &= \left| 1 - \frac52 iy - y^2 \right|^2 - \left| 1 - \frac32 iy - \frac12 y^2 \right|^2 \\ &= (1 - y^2)^2 + \left( -\frac52 y \right)^2 - \left( 1 - \frac12 y^2 \right)^2 - \left( -\frac32 y \right)^2 \\ &= 1 - 2y^2 + y^4 + \frac{25}{4}y^2 - 1 + y^2 - \frac14 y^4 - \frac94 y^2 \\ &= 3y^2 + \frac34 y^4, \end{align*} \end{split}\]so \(E(y) > 0\) for all \(y \in \mathbb{R}\) so the second condition is satisfied and this is an A-stable method.
(c) What are the advantages and disadvantages of using an A-stable implicit Runge-Kutta method over an explicit Runge-Kutta method. In what situation would you choose to implement one over the other?
[5 marks]Solution
Matrix Decomposition Methods#
(a) Use LUP decomposition to solve the following system of equations
\[\begin{split} \begin{align*} 2x_1 + 2x_2 + x_3 &= 11, \\ 4x_1 + x_2 + 3x_3 &= 46, \\ -3x_1 - 5x_2 + 4x_3 &= 14. \end{align*} \end{split}\][10 marks]Solution
Performing partial pivoting on the coefficient matrix
\[\begin{split} \begin{align*} & \begin{pmatrix} 2 & 2 & 1 \\ 4 & 1 & 3 \\ -3 & -5 & 4 \end{pmatrix} \begin{matrix} R_1 \leftrightarrow R_2 \\ \phantom{x} \\ \phantom{x} \end{matrix} & \longrightarrow & \begin{pmatrix} 4 & 1 & 3 \\ 2 & 2 & 1 \\ -3 & -5 & 4 \end{pmatrix} \begin{matrix} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{matrix} & \longrightarrow & \begin{pmatrix} 4 & 1 & 3 \\ -3 & -5 & 4 \\ 2 & 2 & 1 \end{pmatrix} = PA. \end{align*} \end{split}\]Performing the same row swaps on \(I_3\)
\[\begin{split} \begin{align*} & \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{matrix} R_1 \leftrightarrow R_2 \\ \phantom{x} \\ \phantom{x} \end{matrix} & \longrightarrow & \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{matrix} \\ R_2 \leftrightarrow R_3 \\ \phantom{x} \end{matrix} & \longrightarrow & \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} = P. \end{align*} \end{split}\]Computing the LU decomposition of \(PA\)
\[\begin{split} \begin{align*} j &= 1: & u_{11} &= a_{11} = 4, \\ && \ell_{21} &= \frac{1}{u_{11}}a_{21} = \frac14(-3) = -\frac34, \\ && \ell_{31} &= \frac{1}{u_{11}}a_{31} = \frac14(2) = \frac12, \\ j &= 2: & u_{12} &= a_{12} = 1, \\ && u_{22} &= a_{22} - \ell_{21}u_{12} = -5 - \left(-\frac34\right)(1) = -\frac{17}{4}, \\ && \ell_{32} &= \frac{1}{u_{22}}(a_{32} - \ell_{31}u_{12}) = -\frac{4}{17} \left( 2 - \frac12(1) \right) = -\frac{6}{17}, \\ j &= 3: & u_{13} &= a_{13} = 3, \\ && u_{23} &= a_{23} - \ell_{21}u_{13} = 4 - \left( -\frac34 \right)(3) = \frac{25}{4}, \\ && u_{33} &= a_{33} - \ell_{31}u_{13} - \ell_{32}u_{23} = 1 - \frac12(3) - \left( -\frac{25}{4} \right) = \frac{29}{17}, \end{align*} \end{split}\]therefore
\[\begin{split} \begin{align*} L &= \begin{pmatrix} 1 & 0 & 0 \\ -\frac34 & 1 & 0 \\ \frac12 & -\frac{6}{17} & 1 \end{pmatrix}, & U &= \begin{pmatrix} 4 & 1 & 3 \\ 0 & -\frac{17}{4} & \frac{25}{4} \\ 0 & 0 & \frac{29}{17} \end{pmatrix}, & P &= \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}. \end{align*} \end{split}\]Solving \(L\mathbf{y} = P \mathbf{b}\)
\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 & 0 \\ -\frac34 & 1 & 0 \\ \frac12 & -\frac{6}{17} & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} &= \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \begin{pmatrix} 11 \\ 46 \\ 14 \end{pmatrix} = \begin{pmatrix} 46 \\ 14 \\ 11 \end{pmatrix} \quad \implies \quad \begin{aligned} y_1 &= 46 \\ y_2 &= 14 - (-\tfrac34)(46) = \tfrac{97}{2}, \\ y_3 &= 11 - \tfrac12(46) - (-\tfrac{6}{17}) = \tfrac{87}{17}. \end{aligned} \end{align*} \end{split}\]Solving \(U \mathbf{x} = \mathbf{y}\)
\[\begin{split} \begin{align*} \begin{pmatrix} 4 & 1 & 3 \\ 0 & -\frac{17}{4} & \frac{25}{4} \\ 0 & 0 & \frac{29}{17} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} &= \begin{pmatrix} 46 \\ \frac{97}{2} \\ \frac{87}{17} \end{pmatrix} \quad \implies \quad \begin{aligned} x_3 &= \tfrac{17}{29} (\tfrac{87}{17}) = 3, \\ x_2 &= -\tfrac{17}{4} (\tfrac{97}{2} - \tfrac{25}{4}(3)) = -7, \\ x_1 &= \tfrac{1}{4} (46 - (-7) - 3(3)) = 11. \end{aligned} \end{align*} \end{split}\](b) The Gram-Schmidt method for computing the QR decomposition of the matrix \(A\) is:
for \(j = 1 \ldots n\)
\(r_{ij} = \mathbf{q}_i \cdot \mathbf{a}_j, \qquad i = 1, \ldots, j - 1\)
\(\mathbf{u}_j = \mathbf{a}_i - \displaystyle\sum_{i=1}^{j - 1} r_{ij} \mathbf{q}_i\)
\(r_{jj} = \| \mathbf{u}_j \|\)
\(\mathbf{q}_j = \dfrac{\mathbf{u}_i}{r_{jj}}\)
Calculate the QR decomposition of the following matrix using the Gram-Schmidt process and verify that your decomposition is an orthogonal matrix.
\[\begin{split} \begin{align*} A = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \\ 1 & 3 & 4 \end{pmatrix} \end{align*} \end{split}\][10 marks]Solution
\[\begin{split} \begin{align*} j &= 1; & \mathbf{u}_1 &= \mathbf{a}_1 = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \\ && r_{11} &= \| \mathbf{u}_1 \| = \sqrt{3}, \\ && \mathbf{q}_1 &= \frac{\mathbf{u}_1}{r_{11}} = \frac{1}{\sqrt{3}} \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} = \begin{pmatrix} \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \end{pmatrix}, \\ j &= 2: & r_{12} &= \mathbf{q}_1 \cdot \mathbf{a}_2 = \begin{pmatrix} \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \end{pmatrix} \cdot \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} = 2\sqrt{3}, \\ && \mathbf{u}_2 &= \mathbf{a}_2 - r_{12} \mathbf{q}_1 = \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} - 2 \sqrt{3} \begin{pmatrix} \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \end{pmatrix} = \begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix}, \\ && r_{22} &= \| \mathbf{u}_2 \| = \sqrt{2}, \\ && \mathbf{q}_2 &= \frac{\mathbf{u}_2}{r_{22}} = \frac{\sqrt{2}}{2} \begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix} = \begin{pmatrix} -\frac{\sqrt{2}}{2} \\ 0 \\ \frac{\sqrt{2}}{2}\end{pmatrix}, \\ j &=3 & r_{13} &= \mathbf{q}_1 \cdot \mathbf{a}_3 = \begin{pmatrix} \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \end{pmatrix} \cdot \begin{pmatrix} 1 \\ 3 \\ 4 \end{pmatrix} = \frac{8\sqrt{3}}{3}, \\ && r_{23} &= \mathbf{q}_2 \cdot \mathbf{a}_3 = \begin{pmatrix} -\frac{\sqrt{2}}{2} \\ 0 \\ \frac{\sqrt{2}}{2}\end{pmatrix} \cdot \begin{pmatrix} 1 \\ 3 \\ 4 \end{pmatrix} = \frac{3\sqrt{2}}{2}, \\ && \mathbf{u}_3 &= \mathbf{a}_3 - r_{13} \mathbf{q}_1 - r_{23} \mathbf{q}_2 = \begin{pmatrix} 1 \\ 3 \\ 4 \end{pmatrix} - \frac{8\sqrt{3}}{3} \begin{pmatrix} \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \\ \frac{\sqrt{3}}{3} \end{pmatrix} - \frac{3\sqrt{2}}{2} \begin{pmatrix} -\frac{\sqrt{2}}{2} \\ 0 \\ \frac{\sqrt{2}}{2}\end{pmatrix} = \begin{pmatrix} -\frac16 \\ \frac13 \\ -\frac16 \end{pmatrix}, \\ && r_{33} &= \| \mathbf{u}_3 \| = \frac{\sqrt{6}}{6}, \\ && \mathbf{q}_3 &= \frac{\mathbf{u}_3}{r_{33}} = \sqrt{6} \begin{pmatrix} -\frac16 \\ \frac13 \\ -\frac16 \end{pmatrix} = \begin{pmatrix} -\frac{\sqrt{6}}{6} \\ \frac{\sqrt{6}}{3} \\ -\frac{\sqrt{6}}{6} \end{pmatrix}, \end{align*} \end{split}\]therefore
\[\begin{split} \begin{align*} Q &= \begin{pmatrix} \frac{\sqrt{3}}{3} & -\frac{\sqrt{2}}{2} & -\frac{\sqrt{6}}{6} \\ \frac{\sqrt{3}}{3} & 0 & \frac{\sqrt{6}}{3} \\ \frac{\sqrt{3}}{3} & \frac{\sqrt{2}}{2} & -\frac{\sqrt{6}}{6} \end{pmatrix}, & R &= \begin{pmatrix} \sqrt{3} & 2 \sqrt{3} & \frac{8\sqrt{3}}{3} \\ 0 & \sqrt{2} & \frac{3\sqrt{2}}{2} \\ 0 & 0 & \frac{\sqrt{6}}{6} \end{pmatrix}. \end{align*} \end{split}\](c) Explain how LU decomposition used in the calculation of the stage values of an implicit Runge-Kutta method. What advantages does it have for this application over other methods for solving systems of linear equations?
[5 marks]Solution
LU decomposition is commonly used to compute the solution to the linear system required by Newton’s method to solve for the stages values. Since we assume the Jacobian matrix is constant for all Newton iterations, the LU decomposition only needs to be computed once per step, whereas other solution methods for linear systems such as Gaussian elimination would need to be re-computed for each iteration.