Matrix representation of Semi-Implicit Method

Higher order equations, or series of ODEs are more easily solved using the Semi-Implicit Method if represented as a matrix form. Thus,

 

Letters in bold in the above formulation are matrix representations, and I is the identity matrix. Having represented A as a matrix, the Backward Euler can now be formulated as:

 

The matrix formulation of the Backward Euler negates an iterative process, but introduces a matrix inversion. One matrix inversion suffices for the set of ODE’s at each time step, since the partial derivatives are constant. For a set of n ODE’s, the matrix A is of size nXn, and the variable matrix yn is nX1. A matrix multiplication thus needs to be performed to determine y(n+1).

Press et al (Press, Teukolsky, Vetterling & Flannery, Numerical Recipes, The Art of Scientific Computing) recommend LU-decompositions for the matrix inversion process, followed by the LU-back substitution. Subroutines LUDCMP and LUBKSB can be implemented for this process.

The algorithm for the semi-implicit process can be summarized as:

1)    % Define step size (h), initial function y(0), initial time t0, final time tf, eps
2)    nosteps = (tf – t0)/h
3)    for i ← 1 to nosteps

-       Calculate matrix A
-       Invoke subroutine LUCDMP for A
-       B ← LUDCMP (A)
-       Invoke LUBKSB to compute h*B*f(yn)
-       f(yn) ← LUBKSB(h*B*f(yn))
-       y[i] ← y[i-1] + f(yn)
-       t0 ← t0 + h

4)    % Display results