Every time we encounter a problem where things change on two vastly different scales, the equations get stiff.
But even innocent ODEs of second order can easily lead to stiff systems of first order ODEs.
Consider the following example:
(4.79) |
x(t) = A e^{a t} | (4.80) |
(4.81) |
x(t) = A e^{10 t} + B e^{-10 t} | (4.82) |
Even if you set your initial conditions so that A = 0, e.g.,
(4.83) |
We had seen it already when we talked about the instability of the explicit Euler method.
Here is how you can stabilize the solution method regardless of the length of your time step, sic!
Consider the following very simple equation:
(4.84) |
(4.85) |
(4.86) |
(4.87) |
This solution method is said to be unconditionally stable.
The reasoning can be easily expanded to a system of ODEs
as follows. Consider the following:
(4.88) |
(4.89) |
(4.90) |
We have discussed this step in P573.
If all are positive or zero (i.e., C is said to be positive definite) then the method is stable for any step size .
In order to solve the equation we can invert at the first step, and then, assuming that is fixed, simply keep reusing it, because C is constant.
If C = C(t) but still definite positive for every t, or if we're going to change , we may have to invert at every step, or to invoke one of our eigen-value procedures to do the job - and then rotate the solution as need be.
In general
(4.91) |
(4.92) |
If the time step
is a short one, then
(4.93) |
(4.94) |
(4.95) |
The implicit and semi-implicit methods discussed so far are, like the Euler method, first-order accurate only. Higher order implicit methods exist too, but they are prohibitively costly and thus seldom used. It is usually cheaper to shorten the time step while enjoying the stability of a first-order implicit or semi-implicit method at the same time. For a longer time step a semi-implicit method may break to begin with, so the resulting cost would be additionally compounded by having to go back to a fully non-linear case.
Semi-implicit methods are used commonly whenever radiative transfer has to be combined with CFD computations. The reason for this is a very short time scale associated with radiative transfer phenomena compared to CFD time scales. Weather prediction codes and ocean-atmosphere circulation codes are the place where you'll encounter such techniques. Another area is simulation of astrophysical systems, such as accretion disks, supernova explosions, and collisions between neutron stars.
Implicit and semi-implicit methods applied to systems of ODEs and PDEs (the latter can be reduced to the former) require the use of matrix inversion methods, or, at the very least, the use of linear equation solvers. This is what we are going to focus on in the next chapter.