LSOIBT, Implicit Ordinary Differential Equations System Block Tridiagonal Matrices

1 - Description of program or function: LSOIBT is a set of general-purpose FORTRAN routines solver for the initial value problem for ordinary differential equation systems. It is suitable for both stiff and non-stiff systems. LSOIBT treat systems in the linearly implicit form A(t,y) dy/dt = g(t,y), A = a square matrix, i.e. with the derivative dy/dt implicit, but linearly so. It allows A to be singular, in which case the system is a differential-algebraic equation (DAE) system. In that case, the user must be very careful to supply a well-posed problem with consistent initial conditions. LSOIBT, written jointly with C. S. Kenney, solves linearly implicit systems in which the matrices involved are all assumed to be block-tridiagonal. Linear systems are solved by the LU method. The LSOIBT source is commented extensively to facilitate modification. Both a single-precision version and a double-precision version are available. 2 - Methods: This program solves a semi-discretized form of the following system of three PDEs (each similar to a Burgers equation): u(i)t = -(u(1)+u(2)+u(3)) u(i)x + eta(i) u(i)xx (i=1,2,3), on the interval -1.le. x.le. 1, and with time t.ge. 0. The diffusion coefficients are eta(*) =.1,.02,.01. The boundary conditions are u(i) = 0 at x = -1 and x = 1 for all i. The initial profile for each u(i) is a square wave: u(i) = 0 on 1/2.lt. abs(x).le. 1, u(i) = amp(i)/2 on abs(x) = 1/2, u(i) = amp(i) on 0.le. abs(x).lt. 1/2, where the amplitudes are amp(*) =.2,.3,.5. A simplified Galerkin treatment of the spatial variable x is used, with piecewise linear basis functions on a uniform mesh of 100 intervals. The result is a system of ODEs in the discrete values u(i,k) approximating u(i) (i=1,2,3) at the interior points (k = 1,...,99). The ODEs are: (u'(i,k-1) + 4 u'(i,k) + u'(i,k+1))/6 = -(1/6dx) (c(k-1)dul(i) + 2c(k)(dul(i)+dur(i)) + c(k+1)dur(i)) + (eta(i)/dx2) (dur(i) - dul(i)) (i=1,2,3, k=1,...,99), where c(j) = u(1,j)+u(2,j)+u(3,j), dx =.02 = the interval size, dul(i) = u(i,k) - u(i,k-1), dur(i) = u(i,k+1) - u(i,k). Terms involving boundary values (subscripts 0 or 100) are dropped from the equations for k = 1 and k = 99 above