Three naive methods for solving systems of ordinary differential equations
by Elias Hernandis • Published April 2, 2020 • Tagged maths, differential equations, numerical methods
We review three naive iterative methods for solving systems of ordinary differential equations which are widely discussed in undergraduate level scientific computing courses. Implementations are given in MATLAB code.
There are several strategies for solving the systems of ODEs. Some systems, such as the linear ones, can relatively easily be solved analytically. Others, mandate that we use numerical methods. We shall first introduce a method based on an analytic result that will allow us to validate the other methods, at least on the linear case. Then we move on to other methods which support arbitrary systems of ODEs.
This is a work in progress post. This means I have probably not read it even twice. It also means that maybe it will never be completed.
The exponential matrix method
This method is based on the following result from ODE theory. Consider a linear system of ordinary differential equations of the form
where is a function, and is a real-valued matrix. Then the solution to the previous initial value problem is given by the matrix exponential:
It is very convenient that MATLAB provides the function
expm which finds
the value of the matrix exponential numerically. Hence, implementing this
method only requires defining the initial condition vector , the system
matrix and calculating the matrix exponential.
Consider the system
The corresponding MATLAB translation and numerical solution can be computed by
A = [1 2 3; 4 5 6; 7 8 9] u0 = sin(1:3) u = @(t) expm(t .* A) * u0'
A handy way of working with the solution is to define it as a MATLAB function
handle. A function handle is the way to define a mathematical function in
MATLAB. It allows us to evaluate an expression that depends on one or more
variables without having to retype the expression each time. In the previous
example, we can find the value of at say time by evaluating the
In the context of the method of lines, the matrix is some finite difference matrix obtained by approximating the space derivative in the first step of the method. Further examples are given in the application section.
In terms of stability, this method does not present any problems although some systems may yield solutions with very big numbers which overflow MATLAB's floating point representation of numerical values. Since we will only consider small time intervals, this will not be a problem. On the other hand, a downside of this method is that it can only be used when the ODE system at hand is linear.
The Euler-Forward method
Another way of solving systems of ODEs is to apply the so-called Euler forward method. The main idea behind this method is to approximate the derivative by a finite difference
In the context of (linear) systems of ODEs we have that
where is a column vector made from the approximations obtained in the first part of the method of lines. Since there are no restrictions on what we can have on the right hand side of the equality so this very simple method works for non-linear systems as well. To put it more formally, given an initial value problem of the form
the iteration of the Euler forward method is given by
Notice how we have arranged the iteration: given we have an explicit formula for obtaining the next step . This makes implementing this method extremely simple and also gives it an alternative name, the explicit method.
However, this simplicity comes at a price. The EF method is only stable under some circumstances. We will not dive into the general stability considerations of the EF method in this report.
The Euler-Backward method
The stability problems of the EF method lead us to this slightly more advanced version: the Euler backward method. Again, it is based on approximating the derivative by a finite difference, only this time we choose the index for the right-hand side. In the context of a system of the form described in the previous section, the iteration is given by
Contrarily to the EF case, the next step in the iteration appears implicitly in the iteration. This gives this method its alternative name, the implicit method. This means that computing an iteration will be much more complicated: depending on we may even have to solve a non-linear system of equations. If we restrict ourselves to the case of linear systems, then and thus the iteration becomes
i.e. a system of linear equations. Matlab makes this very easy to implement
using the backslash (
"\") operator. Solving a linear system is not as
trivial as evaluating an expression as we did in EF, but this method has the
benefit that it is always stable for the applications we consider.
The IMEX method
When a model is described by a non-linear system of ODEs, it is almost impossible to apply the EB method since it would require a system of non-linear equations for every step of the iteration. In this case, EF is still a possibility, but, on the other hand, we still would like to have some stability warranties. We introduce the IMEX method. The IMEX method owes its name to the fact that it is a combination of the EB (or IMplicit) and the EF (or EXplicit) methods.
To keep things simple we shall focus on ODE systems of the particular form
where we can think of the right-hand side as having a linear part and a potentially non-linear part . The IMEX method uses EB for the "easy" linear part and EF for the "hard" non-linear part. Hence the iteration is given by
This iteration still has an implicit form, but at least we know the right-hand side from the previous step and hence the system is always linear and can be solved with the MATLAB backslash operator. Stability consideration in for this method are much, much more elaborate and are not discussed here.