MATLAB Examples

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; 

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.