The Thomas Algorithm

The standard method of solving a set of simultaneous equations with a computer is to cast the set as a matrix equation. For example, the set 

To find the set of solutions, x, we pre-multiply both sides of the equation 

For this simple example with only three unknowns, we could calculate the product A b by hand however, for systems with a large number of unknowns, this is not practical. For a dense matrix A, one with mostly non-zero values, the normal way to calculate A b is the general Gaussian elimination algorithm [4]. This algorithm takes a number of mathematical operations that is proportional to n where n is the number of unknowns. 

The equation set we consider here (3.36) has n equations so A will be an n X n matrix. As n may be arbitrarily large, the cost of calculating A by Gaussian elimination may be very high — it will take a long time to compute. Fortunately, the equation set has a very useful property of which we can take advantage. For each spatial point i, the equation is of the form 

p and 7 do not vary with i, except at the electrode surface boundary, 

---

The Thomas algorithm begins by a for ward elimination, row by row starting down from the top row (j=l, the condenser stage), to give the following replacements shown in Fig. 13-50Z . For row T. 

The implicit Crank-Nicholson integration method was used to solve the equation. Radial temperature and concentrations were calculated using the Thomas algorithm (Lapidus 1962, Carnahan et al,1969). This program allowed the use of either ideal or non-ideal gas laws. For cases using real gas assumptions, heat capacity and heat of reactions were made temperature dependent. 

When all of the coefficients are known, this can be solved for the concentrations of component i in every stage. A straightforward method for solving a tridiagonal matrix is known as the Thomas algorithm to which references are made in Sec. 13.10, Basis for Computer Evaluation of Multicomponent Separations Specifications. ... 

For implicit schemes, we will obtain a system of linear algebraic equations that must be solved. As mentioned in Example 8.1, one-dimensional diffusion problems generate tri-diagonal matrices, that can be solved for using the Thomas algorithm or other fast matrix routines. Equation (8.83) can be written as... 

The solution methodology of the determinants is similar to that of the well-known Thomas algorithm used for the numerical solution of a differential equation with the finite-difference method [50]. An essential difference from the Thomas algorithm is that the first step ofthe algorithm here is a so-called backward process. This means that the calculation of T starts from the last sublayer, that is, from the Mth sublayer ofthe determinant and it is continued down to the 1st sublayer. Thus, the value of Ti is obtained directly, in the fist calculation step. Then, applying the known value of Ti, the value of Pi can be obtained by means of the fist boundary condition at X= 0, namely ... 

In general case Eqs. (4.60) and (4.61) present infinite sets of the five-term (pentadiagonal) recurrence relations with respect to the index l. In certain special cases (t - 0 or a - 0), they reduce to three-term (tridiagonal) recurrence relations. In this section the sweep procedure for solving such relations is described. This method, also known as the Thomas algorithm, is widely used for recurrence relations entailed by the finite-difference approximation in the solution of differential equations (e.g., see Ref. 61). In our case, however, the recurrence relation follows from the exact expansion (4.60) of the distribution function in the basis of orthogonal spherical functions and free of any seal of proximity, inherent to finite-difference method. Moreover, in our case, as explained below, the sweep method provides the numerical representation of the exact solution of the recurrence relations. 

Tridiagonal systems result from Eq. (10.24a). These can be solved efficiently using the Thomas algorithm, as discussed in the last section of this chapter. 

For tridiagonal matrices, the decomposition of the matrix into a product of a lower and an upper diagonal matrix leads to an efficient algorithm known as the Thomas algorithm. For a system of the form... 

The Thomas algorithm will always converge if the tridiagonal matrix is diagonally dominant. In other words, the matrix is such that... 

A subroutine in FORTRAN code is written below for the Thomas algorithm. 

Up to this point, the treatments have involved reactions for which the discrete form of the reaction-diffusion equations involve only terms in concentration of the species to which the discrete equation applies. That is, if there were two substances involved, O and R as above, then the discrete equation at a point i had terms only in C 0 for species O, and only C R for species R. This made it possible to use the Thomas algorithm to reduce a system like (6.27) to (6.28), treating the two species systems separately. They then get coupled through the boundary conditions. 

When homogeneous reactions take place, it often happens that some of the discrete equations contain terms in concentration for more than the one species, and it is then not generally possible to use the Thomas algorithm to reduce the systems. These systems are said to be coupled. An example will illustrate this situation. 

The above system, although leading to a quadradiagonal system of equations, can still be solved by a smallish extension of the Thomas algorithm [153]. Consider the last two equations of (8.33) and rewrite them, putting the bulk concentration terms on the right-hand side ... 

Thus far, this looks just like the Thomas algorithm for the tridiagonal system, as described above in Sect. 8.3. Prom here on, however, the processes diverge. We need to keep both substitutions for C N and C N 1 and use them in the third-last equation, which contains both. This process is continued backwards, reducing all equations with four unknowns to new ones with just two unknowns. The expressions resulting from this are the following ... 

Coupled equations are those in which some or all of the dynamic equations have terms in more than one of the variables (concentrations). This leads, upou discretisation, to systems of discrete equations that cannot usually be solved using the Thomas algorithm because, no matter how one orders the concentration vectors, the systems correspond to matrix equations that are... 

Standard numerical methods may be used to solve the linear system a number of these are listed in Table 6. Some of these, such as the Thomas algorithm and FIFD, take advantage of the structure of [Ai]. The linear solver is the critical engine behind the simulation. To illustrate this, consider that simulation of more challenging electrode geometries, such as the microband... 

A submatrix is formed at each point on the grid, relating the concentration of each species to the others (kinetically). The material balance equation for all the species may be written with this submatrix down the diagonal - resulting in a block tridiagonal matrix. This may be solved using a matrix version of the Thomas algorithm which requires each submatrix to be inverted (by LU factorization). 

The alternating direction implicit (ADI) method (Peaceman and Rachford, 1955) is a partially implicit method. The equation is rearranged so that one coordinate may be solved implicitly using the Thomas algorithm whilst the others are treated explicitly. If this is done alternately, each coordinate has a share of the implicit iterations and the efficiency (Gavaghan and Rollett, 1990) as well as the stability is improved. The method was used by Heinze for microdisc simulations (Heinze, 1981 Heinze and Storzbach, 1986) and has subsequently been adopted by others (Taylor et al, 1990 Fisher et al., 1997). 

Compute a new set of values of the TJ tear variables by solving simultaneously the set of N energy-balance equations (13-72), which are nonlinear in the temperatures that determine the enthalpy values. When linearized by a Newton iterative procedure, a tridiagonal-matrix equation that is solved by the Thomas algorithm is obtained. If we set gj equal to Eq. (13-72), i.e., its residual, the linearized equations to be solved simultaneously are... 

---