The use of approximate factorization in stiff ODE solvers (Q1298665)

The authors consider an `approximate factorization' method as a replacement for the Newton method for dealing with the nonlinear systems which arise in the solution of stiff ordinary differential equations (ODEs) by implicit methods. Essentially, at each step of the Newton iteration, the Jacobian matrix is replaced by an approximate factorization for which solution is much easier, using an idea which is closely related to the idea of operator splitting (as in, for example, the Godunov-Ryabenki scheme for time-dependent partial differential equations in more than one space dimension). At each step of the now modified Newton scheme, therefore, an inner iteration is set up using the approximated linear equations. The convergence, at each Newton step, of the inner iterations to the exact Newton iterate for that step, is proved under a restriction on the stepsize \(h\) in the initial value solver. This stepsize restriction is typically less severe than the stepsize restriction required by a fixed point iteration. Because it avoids the cost of a full LU factorization, of the Jacobian, however, the new approach is cheaper than a typical Newton iteration. The authors have combined the new scheme into a unified iteration process which is implemented in a code designed to deal automatically with problems which have varying degrees of stiffness over the range of integration. Numerical results on a number of benchmark problems demonstrate the efficacy of the new approach.