Assignment 3 - Implicit Numerical Solvers

Assignment 3 - Implicit Numerical Solvers#

In this exercise we will look at implicit solvers, which can roughly be identified as solvers where the next step \(y_{n+1}\) is dependent on itself. The simplest of these is the backward Euler, also called the implicit Euler method. We can look at the forward and the backward Euler methods side by side, given as:

(112)#\[ \begin{align}\begin{aligned}\text{Forward Euler: } y_{n+1} = y_n + h f(t_n, y_n)\\\text{Backward Euler: } y_{n+1} = y_n + h f(t_{n+1}, y_{n+1})\end{aligned}\end{align} \]

Notice that \(y_{n+1}\) appears on both sides of the equality, this is the essence of what makes it implicit.

Hint

Have you read the code examples in Numerical Methods for ODEs?

Problem 1 - Newton solver#

When dealing with implicit solvers, we have to do something to shake the solution out. If we have a linear system, this generally means that we have to invert a matrix. Consider for instance the linear system \(\dot{y} = Ay\) with the implicit Euler method. This can be written in the matrix-equation form

(113)#\[(I + hA) y_{n+1} = y_n\]

which can be solved by finding the inverse of the matrix \(I - hA\). However, for a non-linear system, where we really need to evaluate some right-hand side function \(f(t_{n+1}, y_{n+1})\), we’re not so lucky. But all is not lost, Isaac Newton and Joseph Raphson found a nice numerical algorithm that can be used to solve even non-linear systems. There are many other methods that can do this, that is, other fixed point methods and root finding algorithms, but Newton-Raphson is a simple but effective one, and usually sufficient and performant for the systems we investigate.

Tasks

  1. Implement the Newton-Raphson method (also simply known as Newton’s method) for an arbitrary function \(f(t, y)\). The inputs to the function should include the initial guess \((t_0, x_0)\), the function \(f(t, x)\) and its Jacobian \(J_f(t, x)\).

    Hint

    In order to verify that you have implemented the method correctly, you can try to find the roots of a polynomial you know the answer to. For instance, the polynomial

    (114)#\[(x_1 + 1)^3 + (x_2 + 1)^3 + (x_3 + 1)^3\]

    should have all its roots at \(x_1 = x_2 = x_3 = -1\). Since it is a polynomial, the analytical evaluation of the Jacobian should be straightforward to compute.

Problem 2 - Implicit Euler and Midpoint#

With a working Newton solver, we can continue with implementing implicit solvers. The two we are going to take a look at are the implicit Euler method and the implicit Midpoint method. These methods have Butcher tables

(115)#\[\begin{split}\text{Implicit Euler:}\quad \begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array} \\ \text{Implicit Midpoint:}\quad \begin{array}{c|c} \frac12 & \frac12\\ \hline & 1 \end{array}\end{split}\]

When you test your implementation of these solvers, use the harmonic oscillator as a test-case. The harmonic oscillator is a simplification of the damped non-linear oscillator from the previous assignment, where we remove the damping term and use the small-angle approximation (\(\sin(\theta) = \theta\)) to remove the non-linearity. This leaves us with the equation

(116)#\[\ddot{\theta} + \omega^2 \theta = 0\]

With initial conditions \(\theta(0) = 1, \dot{\theta}(0) = 0\), we get that the analytical solution is simply

(117)#\[y(t) = \cos(\omega t)\]

This allows us to use the analytical solution directly when determining the error and convergence order of the implemented methods, even when the simulation time is very long. To make things even simpler, we will set \(\omega = 1\).

Note

You may have observed that this is actually a linear system, and thus we could “just invert the matrix” to get the answer and not bother with the non-linear solver. This is correct, but developing a general non-linear solver will be useful for future assignments where things may not longer be the case.

Tasks

  1. Implement implicit Euler using the Newton solver you programmed in the previous problem.

  2. Implement implicit Midpoint using the Newton solver you programmed in the previous problem.

    Hint

    The half-step in the implicit Midpoint can be a bit tricky to figure out. It turns out it is possible to think of the implicit midpoint as an implicit Euler half-step from \(y_n\) to \(y_{n+\frac12}\), and then an explicit Euler half-step from \(y_{n+\frac12}\) to \(y_{n+1}\). This is really where the midpoint method gets its name, since \(y_{n+\frac12}\) is the midpoint.

  3. Find the experimental order of convergence for both methods. Is this unexpected considering that both methods are one-stage methods?

  4. Simulate the harmonic oscillator with a fairly large timestep (maybe \(h = 0.1\)) and for a “long” time (maybe \(t_f = 10\)). Do this for both the explicit Euler, the implicit Euler and the implicit Midpoint method, then compare each solution with the analytical solution. What do the different methods do? Are they able to simulate the system for those solver parameters?

Problem 3 - Stability and stiff systems#

In the previous assignment, we looked at the stability of explicit methods. With a similar treatment of implicit methods, we hope it will be clear why one might be interested in implicit methods. Not only that, but maybe you will get an intuition for what separates A- and L-stability.

In this problem, we will look at a known stiff system described by

(118)#\[\begin{split}\dot{x}_1 = - 0.5 x_1 + 20 x_2 \\ \dot{x}_2 = - 20 x_2\end{split}\]

with initial conditions \(x_1(0) = 0, x_2(0) = 1\).

One way to observe that it is stiff is that the behavior of the system will be very different when \(x_2\) is large and when it is small. It is the case that \(x_2\) will decay very quickly due to the coefficient \(-20\), but during the decay \(x_1\) will increase roughly as fast as \(x_2\) decays. When \(x_2\) no longer contributes much, when it has become small, \(x_1\) will also decay, but much slower due to the coefficient \(-0.5\).

Tasks

  1. Define stiff systems in your own words.

    Hint

    There is no precise mathematical definition of stiffness, and most attempts at defining it defer to general statements about when the “phenomenon of stiffness” appears in a numerical solution. Reflect over what you think it means for a system to “be stiff”.

  2. Find the stability functions for the implicit Euler and the implicit Midpoint method.

  3. Plot the stability functions for both methods and conclude whether they are A-stable. Also, are any of the methods L-stable?

  4. Use both the explicit Euler, implicit Euler and the implicit Midpoint methods to numerically solve the stiff system Increase the stepsize from \(h = 0.01\) to \(h = 0.1\) and describe what you observe. Include a plot of \(x_1\) over time for each of the methods for \(h = 0.1\).

  5. What do you think is the main concerns one should address when choosing a numerical solver?

    Hint

    Some keywords may be accuracy, stability, computational cost and solver properties.

  6. (Optional) If you implemented an adaptive controller for the timestep in the previous assignment, try to use it to simulate the system and see what timesteps it selects. Are the selected timesteps reasonable?

Problem 4 (Optional) - General IRK methods#

By now, we have dealt with a large number of Runge-Kutta methods, all of which can be represented with a Butcher table. Maybe the time is right to implement a general routine that can take any Butcher table and perform the corresponding RK method? In general, an RK method will then look like

(119)#\[\begin{split}\begin{array}{c|c} \vec{c} & A \\ \hline & \vec{b} \end{array}\end{split}\]

where \(A\) defines the linear combination of the stages inside the ODE right-hand side evaluations, \(b\) is the final linear combination of stages to get the next step and \(c\) is the vector of appropriate time offsets. Until now, the time offsets have not been very important, since our systems have been autonomous, but often in real life, the forces we’re subjected to change not only in space, but also time, so it is important.

Let us consider the system

(120)#\[\begin{split}\dot{x}_1 = -\lambda_1 x_1 + A_1 \cos(\omega_1 t), \\ \dot{x_2} = -\lambda_2 x_2 + A_2 \cos(\omega_2 t).\end{split}\]

It is not a system of homogeneous ODEs since it has the cosine driving term, which makes it non-autonomous. The ODEs are decoupled, so they can be solved one at a time using an integrating factor, which gives the solution

(121)#\[x_i(t) = \frac{A_i (\lambda_i \cos(\omega_i t) + \omega_i \sin(\omega_i t))}{\lambda_i^2 + \omega_i^2} + e^{-\lambda_i t} \left( y_0 - \frac{A_i\lambda_i}{\lambda_i^2 + \omega_i^2} \right)\]

To test that the general IRK implementation is correct, you may of course compare the performance between your hard-coded previous implementations to the general scheme, but you may also test it on an RK method you haven’t implemented yet: the two-stage Gauss-Legendre collocation method. This method is given by the vectors and matrices:

(122)#\[\begin{split}A = \begin{pmatrix} \frac14 & \frac14 - \frac{\sqrt3}{6} \\ \frac14 + \frac{\sqrt3}{6} & \frac 14 \end{pmatrix}, \quad \vec{b} = \begin{pmatrix} \frac12 \\ \frac12 \end{pmatrix}, \quad \vec{c} = \begin{pmatrix} \frac12 - \frac{\sqrt3}{6} \\ \frac12 + \frac{\sqrt3}{6} \end{pmatrix}.\end{split}\]

Since it is a collocation method, it should have twice the convergence order one would expect. As a two-stage method, it should then be a fourth order accurate method, which you can numerically verify by computing the experimental order of convergence.

Tasks

  1. Implement a function that applies a general IRK method to a system defined by its right-hand side function \(f(t, x)\) and the Jacobian of this function \(J_f(t, x)\). The general IRK function should handle state vectors of arbitrary dimensions and an arbitrary number of implicit or explicit stages.

  2. Both explicit and implicit RK methods can be described by Butcher tables. However, it is wasteful to use the Newton on a known explicit scheme, since it can be solved with a forward pass through the stages. Extend your general code to detect whether the Newton solver is required or not.

    Hint

    One only needs to investigate the \(A\)-matrix do determine whether the RK method is implicit or explicit. What is the condition \(A\) must satisfy?

  3. There is a third “type” of Runge-Kutta methods hidden in this general approach, the so-called diagonally-implicit Runge-Kutta methods, or DIRK among friends. The major innovation here is that while the DIRK is implicit, any stage \(p\) is only dependent on the \(p-1\) stages before it. This means that we do not need to solve for all stages simultaneous, but can solve one at a time, until we get to the last one. Extend your general IRK function to handle this case, as well, possibly saving some precious computational time.

    Hint

    There is also a special condition the \(A\)-matrix must satisfy here if the scheme is to be a DIRK. What is it?

Contents