Reduced storage matrix methods in stiff ODE systems

Stiff initial-value ODE methods (e.g. BDF) normally require the Jacobian matrix as part of a Newton-like iteration within each implicit time step, and require it to be formed and stored explicitly, whether the linear algebraic method used is direct or iterative. But, a Krylov subspace iterative linear system method, which involves the system matrix only in operator form, can be made part of an inexact Newton method within such a stiff ODE method in a matrix-free manner, requiring no explicit Jacobian matrix storage. Such combinations, using BDF methods, have been implemented with Arnoldi iteration, GMRES (the Generalized Minimum RESidual method), and CG (the Conjugate Gradient Method). Various practical matters (scaling, starting, stopping, etc.) are dealt with in the stiff ODE context. In the context of general nonlinear algebraic systems, we provide some theoretical foundation for the combined Newton-Krylov method by giving convergence results that include errors due to the difference quotient approximation to the linear operator. Earlier tests showed matrix-free methods to be quite effective, at least when the spectrum of the kproblem Jacobian is rather tightly clustered. To improve their robustness, we have added preconditioning, in an experimental solver called LSODPK, which we tested on ODE systems that arise from time-dependent PDE systems by the method of lines. Preconditioner matrices can be formed from the interaction or reaction terms, from the spatial transport terms, or both (as in operator splitting). The additional matrix storage can be reduced greatly by grouping the diagonal blocks in a natural way. The methods appear to be quite effective in improving both speed and storage economy over traditional stiff methods, and over matrix-free methods without preconditioning.