The use of sparse matrix technique in the numerical integration of stiff systems of linear ordinary differential equations

Abstract

Consider the system of the ordinary differential equations of first order: y′= A( t) y+ b( t), y ϵ R s , t ϵ [ a,b] ⊂ R, y(a)= y O , y O given. Denote by λ i (( t) , i=1(1) s, t ϵ [ a,b], the eigenvalues of matrix A(t). Let υ i = Re(λ i ( t)), ν i (t)= Im(λ i (t)), λ=max (|λ i (t)|), itν=max(|λ i (t)|), i=1(1) s, t ϵ [ a,b]. Assume that: (1) ν is large, (2) λ i ( t) α0, i=1(1) s, t ϵ[ a,b], (3) ϵ ⪡v. Some problems arising in the spectroscopic resonance theory satisfy the above coditions. A modified Nθrsett 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(t) is called sparse in this case) and that matrix A(t) is large. An attempt to exploit the sparsity of matrix A( t) during the numerical solution of the system of ordinary differential equations is carried out in this paper. A detailed analysis of the impletmentation of some sparse matrix techniques (based in a Gustavson scheme and a generalized Markowitz pivotal strategy) in the integration of linear systems of ordinary differential equation is presented. The possibility ofimproving 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 recommendations are drawn.