1.8 Multi-Step Methods

1.8.3 Backwards Differentiation Methods

Measurable Outcome 1.12, Measurable Outcome 1.15, Measurable Outcome 1.17, Measurable Outcome 1.18, Measurable Outcome 1.19

Backwards differentiation methods are some of the best multi-step methods for stiff problems. The backwards differentiation formulae are of the form,

The coefficients for the first through fourth order methods are given in the table below.

The stability boundary for these methods are shown below. As can be seen, all of these methods are stable everywhere on the negative real axis, and are mostly stable in the left-half plane in general. Thus, backwards differentiation work well for stiff problems in which stong damping is present.

MATLAB^®'s ODE Integrators

MATLAB has a a set of tools for integration of ODE's. We will briefly look at two of them: ode45 and ode15s. ode45 is designed to solve problems that are not stiff while ode15s is intended for stiff problems. ode45 is based on a four and five-stage Runge-Kutta integration (discussed in Section 1.9), while ode15s is based on a range of highly stable implicit integration formulas (one option when using ode15s is to use the backwards differentiation formulas). As a short illustration on how these MATLAB ODE integrators are implemented, the following script solves the one-dimensional diffusion problem from Section 1.7.1 using either ode45 or ode15s. The specific problem we consider here is a bar which is initially at a temperature \(T_{init} = 400 K\) and at \(t=0\), the temperature at the left and right ends is suddenly raised to \(800 K\) and \(1000 K\), respectively.

% MATLAB script: dif1d_main.m
%
% This code solve the one-dimensional heat diffusion equation
% for the problem of a bar which is initially at T=Tinit and
% suddenly the temperatures at the left and right change to
% Tleft and Tright.
%
% Upon discretization in space by a finite difference method,
% the result is a system of ODE's of the form,
%
% u_t = Au + b
%
% The code calculates A and b.  Then, uses one of MATLAB's
% ODE integrators, either ode45 (which is based on a Runge-Kutta
% method and is not designed for stiff problems) or ode15s (which
% is based on an implicit method and is designed for stiff problems).
%

clear all; close all;

sflag = input('Use stiff integrator? (1=yes, [default=no]): ');

% Set non-dimensional thermal coefficient
k   = 1.0; % this is really k/(rho*cp)

% Set length of bar
L = 1.0; % non-dimensional

% Set initial temperature
Tinit = 400;

% Set left and right temperatures for t>0
Tleft = 800;
Tright = 1000;

% Set up grid size
Nx = input(['Enter number of divisions in x-direction: [default=' ...
            '51]']);
if (isempty(Nx)),
  Nx = 51;
end

h  = L/Nx;
x  = linspace(0,L,Nx+1);

% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
Tmax = 0.5;

% Initialize a sparse matrix to hold stiffness & identity matrix
A = spalloc(Nx-1,Nx-1,3*(Nx-1));
I = speye(Nx-1);

% Calculate stiffness matrix

for ii = 1:Nx-1,

  if (ii > 1),
    A(ii,ii-1) = k/h^2;
  end

  if (ii < Nx-1),
    A(ii,ii+1) = k/h^2;
  end

  A(ii,ii) = -2*k/h^2;

end

% Set forcing vector
b = zeros(Nx-1,1);
b(1)    = k*Tleft/h^2;
b(Nx-1) = k*Tright/h^2;

% Set initial vector
v0 = Tinit*ones(1,Nx-1);

if (sflag == 1),

  % Call ODE15s
  options = odeset('Jacobian',A);
  [t,v] = ode15s(@dif1d_fun,[0 Tmax],v0,options,A,b);

else

  % Call ODE45
  [t,v] = ode45(@dif1d_fun,[0 Tmax],v0,[],A,b);

end

% Get midpoint value of T and plot vs. time
Tmid = v(:,floor(Nx/2));
plot(t,Tmid);
xlabel('t');
ylabel('T at midpoint');

As can be seen, this script pre-computes the linear system \(A\) and the column vector \(b\) since the forcing function for the one-dimensional diffusion problem can be written as the linear function, \(f = Av + b\). Then, when calling either ODE integrator, the function which returns \(f\) is the first argument in the call and is named, dif1d_fun. This function is given below:

% MATLAB function: dif1d_fun.m
%
% This routine returns the forcing term for
% a one-dimensional heat diffusion problem
% that has been discretized by finite differences.
% Note that the matrix A and the vector b are pre-computed
% in the main driver routine, dif1d_main.m, and passed
% to this function.  Then, this function simply returns
% f(v) = A*v + b.  So, in reality, this function is
% not specific to 1-d diffusion.

function [f] = dif1d_fun(t, v, A, b)

f = A*v + b;

As can be seen from dif1d_fun, \(A\) and \(b\) have been passed into the function and thus the calculation of \(f\) simply requires the multiplication of \(v\) by \(A\) and the addition of \(b\).

The major difference between the implementation of the ODE integrators in MATLAB and our discussions is that MATLAB's implementations are adaptive. Specifically, MATLAB's integrators estimate the error at each iteration and then adjust the timestep to either improve the accuracy (i.e. by decreasing the timestep) or efficiency (i.e. by increasing the timestep).

The results for the stiff integrator, ode15s are shown in Figure 1.22.

These results look as expected (note: in integrating from \(t=0\) to \(t=0.5\), a total of 64 timesteps were taken).

The results for the non-stiff integrator are shown in Figure 1.23 and in a zoomed view in Figure 1.24. The presence of small scale oscillations can be clearly observed in the ode45 results. These oscillations are a result of the large negative eigenvalues which require small \({\Delta t}\) to maintain stability. Since the ode45 method is adaptive, the timestep automatically decreases to maintain stability, but the oscillatory results clearly show that the stability is barely achieved. Also, as a measure of the relative inefficiency of the ode45 integrator for this stiff problem, note that 6273 timesteps were required to integrate from \(t=0\) to \(t=0.5\).

One final concern regarding the efficiency of the stiff integrator ode15s. In order for this method to work in an efficient manner for large systems of equations such as in this example, it is very important that the Jacobian matrix, \({\partial f}/{\partial u}\) be provided to MATLAB. If this is not done, then ode15s will construct an approximation to this derivative matrix using finite differences and for large systems, this will become a significant cost. In the script dif1d_main, the Jacobian is communicated to the ode15s integrator using the odeset routine. Note: ode45 is an explicit method and does not need the Jacobian so it is not provided in that case.