Krylov approximations for matrix square roots in stiff boundary value problems

Recently, we proposed an algebraic difference scheme, with extended stability properties, for linear boundary value problems involving stiff differential equations of first order. Here, an efficient approximation scheme is presented for matrix square roots, which provides the stabilization of that scheme in case of stiffness. It combines the use of low-rank matrix approximations from projections onto Krylov subspaces with an accelerated sign iteration for the matrix square root. The Krylov approximation, being accurate in eigenspaces with large eigenvalues, preserves the stability of the scheme, and the 0(n3) square root computation need be performed only in lower dimension. Operation counts and numerical results show that the effort for the numerical scheme is essentially proportional to the number of stiff components, but not to the norm of the coefficient matrix. Approximation properties of low-rank Krylov matrices, which may be of independent interest, are analyzed. 1. THE SQRT ONE-STEP DIFFERENCE SCHEME A standard approximation for the differential equation in the linear boundary value problem (BVP) u'(x) = A(x) u(x)+ g (x) X E [0 1], U(X) ER, Bou(O) + B u(l) = 8, is the trapezoidal rule on a suitable grid, 0 = xo 1 , this restriction leads to unacceptable small stepsizes in regions where the solution is smooth. This situation is usually referred to as "stiffness." In certain cases, however, a relaxation Received March 27, 1990; revised December 5, 1990. 1991 Mathematics Subject Classification. Primary 65L10, 65L20, 65F30.