This file contains an explanation of the difference between implicit and explicit time integration schemes. The content is intended for those who want to learn a bit more than what the ForwardandBackwardEulerExplorer GUI can offer.This file contains an explanation of the difference between implicit and explicit time integration schemes. The content is intended for those who want to learn a bit more than what the ForwardandBackwardEulerExplorer GUI can offer.
clear all;
Start with a simple linear ODE,

From ODE stability analysis we know that our eigenvalues need to be in the left half plane. In this case we only have one eigenvalue,
. So for stability we need
.
Now we can analyzed the staiblity for finite
. Given the general for of an ODE

the Forward Euler method is expressed as:

Applying Forward Euler to
gives,

The solution has the following form,

where
is the amplification factor such that,

We want to determine under what conditions
since this would mean that
would grow unbounded as
. One way to do this is to solve for the stability boundary for which
. To do this, let
(since
) where
. Making this substitution into the amplification factor we get:

For stability using Forward Euler we need

theta = linspace(0,2*pi,100); lambda_delta_t = exp(i*theta)-1; plot(real(lambda_delta_t),imag(lambda_delta_t)) grid on; axis([-3 3 -3 3]) text(-1.0,.3,'stable','HorizontalAlignment','center'); xlabel('real(\lambda\Delta t)') ylabel('imag(\lambda\Delta t)') title('Stability Region for Forward Euler')
For a given problem, i.e. with a given
, the timestep,
, must be chosen so that the algorithm remains stable for
.
As an example let's take
,
. This ODE has an exact solution of
which we can compare to the results of Forward Euler for
and
.
v0 = 1; lambda = -2; tex = linspace(0,20); uex = exp(-2*tex); dt = [0.1 0.9 1.1]; for ind=1:3 N = ceil(20/dt(ind)); t{ind} = dt(ind)*[0:N]; v{ind} = (1+lambda*dt(ind)).^[0:N] * v0; end subplot(3,1,1) plot(tex,uex,t{1},v{1}) xlabel('t'); ylabel('v'); xlim([0 20]); text(15,0.5,'\Delta t = 0.1') subplot(3,1,2) plot(tex,uex,t{2},v{2}) xlabel('t'); ylabel('v'); xlim([0 20]); text(15,0.5,'\Delta t = 0.9') subplot(3,1,3) plot(tex,uex,t{3},v{3}) xlabel('t'); ylabel('v'); xlim([0 20]); text(15,25,'\Delta t = 1.1')
If we have a stiff system with large negative real eigenvalues, using an explicit time integration scheme can be very inefficent as we will need to use very small time steps. A more efficient approach to numerically integrate a stiff problem would be to use a metho with eigenvalue stability for large negative real eigenvalues. Implicit methods often have good stability along the negative real axis. The simplice implicit metho is the Backward Euler Method,

The amplification factor for Backword Euler is

and for stability using Backward Euler we need

which is shown below
clf; theta = linspace(0,2*pi,100); lambda_delta_t = exp(i*theta)+1; plot(real(lambda_delta_t),imag(lambda_delta_t)) grid on; axis([-3 3 -3 3]) text(1.0,.3,'unstable','HorizontalAlignment','center'); xlabel('real(\lambda\Delta t)') ylabel('imag(\lambda\Delta t)') title('Stability Region for Backward Euler')
Returning to the example problem,
,
. We can compare to the results of Backward Euler for
and
.
v0 = 1; lambda = -2; tex = linspace(0,20); uex = exp(-2*tex); dt = [0.1 0.9 1.1]; for ind=1:3 N = ceil(20/dt(ind)); t{ind} = dt(ind)*[0:N]; v{ind} = (1/(1-lambda*dt(ind))).^[0:N] * v0; end subplot(3,1,1) plot(tex,uex,t{1},v{1}) xlabel('t'); ylabel('v'); xlim([0 20]); text(15,0.5,'\Delta t = 0.1') subplot(3,1,2) plot(tex,uex,t{2},v{2}) xlabel('t'); ylabel('v'); xlim([0 20]); text(15,0.5,'\Delta t = 0.9') subplot(3,1,3) plot(tex,uex,t{3},v{3}) xlabel('t'); ylabel('v'); xlim([0 20]); text(15,0.5,'\Delta t = 1.1')
So implicit methods are WAY better. Why don't we use them all the time?
Consider the Forward Euler method applied to
where
is a
matrix,

In this explicit algorithm, the largest computational cost is the matrix vector multiply,
, which is an
operation. Now for the implicit Backward Euler method,

Re-arraging to solve for
gives:


Thus to vind
for Backward Euler requires the solution of a
system of equations which is an
cost. As a result, for large systems, the cost of the
linear solution may begin to outweigh the benefits of the larger timesteps that are possible when using implicit methods.
