Commtrrs & Chemistry Vol. 4, pp. I-12 Pagmon Press Ltd.. IWO. Printed in Oreat Britain zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCB THE USE OF SPARSE MATRIX TECHNIQUE IN THE NUMERICAL INTEGRATION OF STIFF SYSTEMS OF LINEAR ORDINARY DIFFERENTIAL EQUATIONS KJELD SCHAUMBURG Department of Chemical Physics, University of Copenhagen, The H.C. (arsted Institute, 5, Universiietsparken, DK 2100 Copenhagen 0, Denmark JERZY WASNIEWSKI The Regional Computing Center at the University of Copenhagen 5, Vermundsgade, DK-2100 Copenhagen 0, Denmark and ZAHARI ZLATEV Department of Mathematics and Statistics, The Royal University for Veterinary and Agriculture 4O, Thorvald- sensvej, DK-I871 Copenhagen V, Denmark (Received 23 July 1979) Abstract-Consider the system of ordinary differential equations of first order: y’ = A(f)y + b(t), y E R”, t E [a, b] c R, y(n) = yo, y. given. Denote by h,(f), i = I(I)s, I E [(I, b], the eigenvalues of matrix A(t). Let 140 = ReMOk 40 = IN%(O), A = max(lW p = max (MOI), Y = ma&(t)l), zyxwvutsrqponmlkjihgfedcbaZYX i = l(l)s, t E [n, b]. Assume that: (I) A is large, (2) &I) ~0, i = l(l)s, t E [a, h], (3) B-Z v. Some problems arising in the spectroscopic resonance theory satisfy the above conditions. A modified Nprsett method has successfully been used to solve numerically the above system in a previous paper (Schaumburg et al., 1979). The same method is used in this paper also. However, it is now assumed that many elements in matrix A(t) are equal to zero (matrix A(I) is called sparse in this case) and that matrix A(t) is large. An attempt to exploit the sparsity of matrix A(l),during the numerical solution of the system of ordinary differential equations is carried out in this paper. A detailed analysis of the implementation of some sparse matrix techniques (based on a Gustavson storage scheme and a generalized Markowitz pivotal strategy) in the integration of linear systems of ordinary differential equations is presented. The possibility of improving the results by the use of iterative refinement and large values of a special parameter (called a drop-tolerance) is also discussed. Five different algorithms are compared and some conclusions and recom- mendations are drawn. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 1. STATEMENT OF THE PROBLEM Consider the first order linear system of ordinary differential equations (this problem is called LIVP 1 in Schaumburg et al., 1979) Y’ = A(Oy + b(t), y(a) = Y,, hgiven, (1.1) where zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA t E [n, b] C R, y E R”, b(t) E R”, (1.2) s is a positive integer; A(t) E R’“‘, i.e. A(t)={a,(t)};,z,, aif E R. (1.3) Denote by Ai( i = l( l)s, t E [a, 61, the eigenvalues of matrix A(t). Let pi(r)= Re{hi(f)} and vi(t)= Im{(ni(t)}, i = I(l)s, t E [a, b]. Consider the following three numbers A = max (IAi(t) i = I(l)& t E la, b]: (1.4) CL= max(jIli(t)/}, i = l(l)s, 1 E [a, b]; (1.5) v =max fIw(r)ll, i = I(l)s, r E [a, b]. (1.6) 1 Assume that: A % M, where M is large; (1.7) /A;(?)lO, i = l(l)& t E [a, b]; (1.8) p 4 ZJ. (1.9) The problems arising in the spectroscopic resonance theory (see, e.g. Bildspre et al., 1976; Jacobsen et al., 1976; Jacobsen & Schaumbure. 1976) are of tvne (1.1) where (1.7H1.9) are satisfied. here weA 1 shbuld emphasize that all our conclusions were drawn from experiments with spectroscopic LIVP’s 1, however, we believe that these conclusions are true for any LIVP 1 satisfying (1.7)-( 1.9). The performance of different codes, available at RECKU (The Regional Computing Centre at the Uni- versity of Copenhagen and NEUCC (The Northern Europe University Computing Centre), in the solution of LIVP’s I arising from the spectroscopic resonance theory has been investigated by Schaumburg & Was- niewski (1978). They have found that the Ncirsett code SIRKUS, (based bn semiexplicit two&e Runge- Kutta formulae (see Ngrsett, 1974), is the most efficient code among the codes tested, see more details in Schaumburg & Wasniewski (1978).