DASSL and DASPK

DASSL uses backward differentiation formula (BDF) methods [6] to solve a system of DAEs or ODEs. The methods are variable step-size variable order. The system of equations in DASSL is written in an implicit ODE form like

F(t,y,y') = 0,

where y' denotes the time derivatives of y. The BDF methods used in DASSL require the solution of a large system of nonlinear equations

$$F(t_n, y_n, \alpha_n y_n+\beta_n) = 0$$

on each time step. Here, $\alpha_n$ and $\beta_n$ are scalars which depend on the method and stepsize. In DASSL, this system is solved by a modified Newton iteration. Each iteration of the Newton method requires the solution of a linear system

Ay_n^(k+1) = b_n^(k),

where the matrix A is given by

$$A = \alpha_n \frac{\partial F}{\partial y'} + \frac{\partial F}{\partial y}.$$

The one-dimensional PDEs studied in this Chapter generates a matrix which is block tridiagonal. In DASSL, this linear system is solved via a banded direct solver. Because the CPU cost to solve this linear system is proportional to the bandwidth of the matrix, this solver is quite efficient if the bandwidth of the matrix is relatively small. Different moving mesh strategies result in different bandwidths, which is a very important factor in considering the efficiency of the method. The reader can refer to [13] for more implementation details.

In the past few years, a new DAE solver DASPK has been developed by Brown, Hindmarsh and Petzold [3]. DASPK used the preconditioned incomplete GMRES method to solve the linear system at each Newton iteration. A useful property in this context of the GMRES iteration is that it does not require the matrix explicitly, but only requires the matrix-vector product Av for a given vector v. In DASPK, Av is approximated by the difference

$$Av \simeq \d{F(t, y+\delta v, \alpha(y+\delta v)+\beta)-F(t,y,\alpha y+\beta)} {\delta}.$$

The basic method in DASPK can often be made much more robust and efficient with the use of a preconditioner. In practice, a block diagonal preconditioner like block-SSOR or block-Jacobi are often very effective for reaction-diffusion equations. For general equations, an incomplete LU (ILU) factorization preconditioner is often used. For two-dimensional partial differential equations, the bandwidth of the matrix involved in the linear system at each Newton iteration is quite large. Thus the banded solver used in standard DASSL is quite inefficient. DASPK with a simple preconditioner as described above is very effective in this case.