Solution of ODE’s using the State-Space Method

The state-space representation of

   provides a no-fuss approach to converting problems with 2^nd order ODE’s. The derivative equations take the form

 

The derivative definition thus requires a matrix inversion for A. 

 

Matrix inversion can be performed using techniques such as Gauss-Jordan. Subroutine gaussj from Numerical Recipes can be effortlessly used for the inversion. Rather than perform the matrix inversion, Press et al recommend a LU-decomposition followed by an LU-back substitution to solve for the differential form. They suggest that using the LU-decomposition is “a factor of 3 better than subroutine gaussj”.

Following the recommendation of Press et al, the preferred way to solve for the differential form is to assume it to be a set of linear equations of the form

 

The derivative form is then obtained by performing an LU-decomposition of Q, followed by the LU-backsubstitution of the LU version of Q with P. Subroutines ludcmp and lubksb work superbly efficient for this case.

The Semi-Implicit Euler's Method

As noted earlier, the solution of ODE’s using implicit methods (such as the Backward Euler’s method) requires an iterative procedure at each time step. Let’s take a second look at the Backward Euler’s method:

 

 Atkinson et al (Atkinson, Han & Stewart, Numerical Solution of ODE’s) recommend a two-step approach to solving the above equation in order to avoid iterations. The two-step process can be summarized as follows:

 

 

Atkinson et al also indicate that by solving it as a two-step process, the method is no longer absolutely stable. That is, while the original Backward-Euler’s method would converge at a solution for any step size, the two-step process would not.

Methods that are used to solve for implicit formulations without an iterative process are termed as “Semi-Implicit” methods. The computational effort for iterations at each time step is avoided in the semi-implicit method, while at the same time “reasonable” stability can be guaranteed. Press et al (Press, Teukolsky, Vetterling & Flannery, Numerical Recipes, The Art of Scientific Computing) state “It is not guaranteed to be stable, but it usually is…”.