UCES

"Efficient Numerical Solution of Tridiagonal Algebraic Systems"

Course:   Methods and Analysis in UCES

Chapter:   3.0   Numerical Solution of Boundary Value Problems
Section:   3.2   Efficient Numerical Solution of Tridiagonal Algebraic Systems
Introduction:   operation counts and memory

Applied Area:   heat transfer in a fin, continuum model

Model:   finite difference or discrete model

Method:   tridiagonal algorithm

Implementation:   Maple, Matlab

Assessment:   1D versus 2D diffusion, existence, stability

Homework:

1. Math:   Calculus 2, algebraic systems
2. Computing:   Maple
3. Science:   Physics

Objectives:

1. Math:   Analysis of tridiagonal matrices
2. Computing:   Efficient solution of tridiagonal problems
3. Science:   Heat transfer in a fin, effect of convective boundary

General Information: A careful derivation and analysis of the tridiagonal algorithm is given. The 1D model of a cooling fin is presented. This model will later be generalized to the 2D and 3D diffusion models.

Contact Person:   R. E. White, NCSU, white@math.ncsu.edu

Revision Date:   8-3-94

Copyright ©1994, R. E. White. All rights reserved.

Return to Table of Contents

3.2   Efficient Numerical Solution of Tridiagonal Algebraic Systems

In the thin wire problem we derived a tridiagonal matrix which was used in the finite difference approximation of the differential equation. It is very common to obtain either similar tridiagonal matrices or more complicated matrices which have blocks of tridiagonal matrices. This lecture is concerned with a very efficient version of the Gaussian elimination algorithm for the solution of tridiagonal algebraic systems. The full version of a Gaussian elimination algorithm for n unknowns requires order n³ /3 operations and order n² storage locations. By taking advantage of the number of zeros and their location, the Gaussian elimination algorithm for tridiagonal systems can be reduced to order 5n operations and order 6n storage locations!

3.2.2.   Applied Area.

Consider a mass which must be cooled by transmitting heat from the mass to the surrounding region. Examples include an electrical amplifier, a transformer on a power line, or a gasoline engine. One way to do this is to attach fins to this mass so that the surface area that transmits the heat will be large. We wish to be able to model this heat flow so that one can determine whether or not a particular configuration will sufficiently cool the mass.

In order to start the modeling process we will make some assumptions that will simplify the model. In later lectures we will return to this model and reconsider some of these assumptions. First, assume no time dependence and that the fin is long so that the temperature is only a function of the distance from the mass to be cooled. Thus, there is diffusion in only one direction.

Figure: Cooling Fin

Second, assume the heat lost through the surface of the fin is similar to Newton's law of cooling:

heat loss through the   lateral surface of the fin	= (area) (time interval) c(usur - u)
	= h(2W + 2T) t c(usur - u)

where u_sur is the surrounding temperature and the c reflects the ability of the fin's surface to transmit heat to the surrounding region.

Third, assume heat diffuses in the x direction according to Fourier's heat law where K is the thermal conductivity:

For interior volume elements with x+h < L = 1 (the tip of the fin),

0	(heat through lateral surface)   + (heat diffusing through front)     - ( heat diffusing through back)
=	h(2W + 2T)t c(usur - u)   + TW t Kux (x + h)     - TW t Kux (x).	(1)

For the tip of the fin with x + h = L, we use Ku_x(L) = c(u_sur - u(L)) and

0	(heat through lateral surface of tip)   + (heat diffusing through front)     - ( heat diffusing through back)
=	(h/2)(2W + 2T)t c(usur - u(L - h/2))   + TW t c(usur - u(L))     - TW t Kux (L - h/2).	(2)

Note, the volume element near the tip of the fin has one half of the volume of the interior elements.

These are only approximations because the temperature changes continuously with space. In order to make these approximations in (1) and (2) more accurate, we divide by h t TW and let h go to zero:

0 = (2W + 2T)/(TW) c(usur - u) + (Kux )x .

(3)

Let C = (2W + 2T)/(TW) c and F = C u_sur . The continuum model is given by the following differential equation and two boundary conditions.

-(Kux )x + Cu = F,	(4)
u(0) = given, and	(5)
Kux (L) = c(usur - u(L)).	(6)

The boundary condition in (6) is often called a derivative or flux or Robin boundary condition.

3.2.3.   Model.

The above derivation is useful because (2) and (3) suggest a way to discretize the continuum model. Let u_i be an approximation of u(ih) where h = L/n. And approximate the derivative u_x(ih + h) by (u_i+1 - u_i )/h. Then equations (2) and (3) yield the finite difference approximation, a discrete model, of the continuum model (4), (5) and (6):

    For 1 i n - 1:

-[K(ui+1 - ui )/h - K(ui - ui-1 )/h] + h C ui = h F(ih) with u0 given,

(7)

    For i = n:

-[c(usur - un)/h - K(un - un-1 )/h] + (h/2) C un = (h/2) F(nh).

(8)

The discrete system (7) and (8) may be written in matrix form. For ease of notation we let n = 4, multiply (7) by h and (8) by 2h, B = 2K + h²C so that there are 4 equations and 4 unknowns.

      Bu1 - Ku2                   = h2f1 + Ku0

     -Ku1 + Bu2 - Ku3             = h2f2

           -Ku2 + Bu3 -       Ku4 = h2f3

                -2Ku3 + (B+2hc)u4 = h2f4 + 2chusur

The matrix form of this is

AU = F where

  A is, in general, n×n matrix and

  U and F are n×1 column vectors.

Here with n = 4 we have

3.2.4.   Method.

The solution can be obtained by either using the tridiagonal (or Thomas) algorithm, or using a solver that is provided with your computer software. Let us consider the tridiagonal system Ax = d where A is an n×n matrix and x and d are n×1 column vectors. We assume the matrix A has components as indicated:

.

In previous section we used the Gaussian elimination algorithm on tridiagonal matrices, and we noted the matrix could be factored into two matrices A = LU where L has nonzero components only in its diagonal and subdiagonal, and U has nonzero components only in its diagonal and superdiagonal. For the above 4×4 matrix this is

.

The plan of action is

(i). solve for _i and _i in terms of a_i , b_i and c_i by matching components of the above matrix equation,

(ii). solve Ly = d and

(iii). solve Ux = y.

• Step (i):

  For i = 1,

  a₁ = ₁*1 and c₁ = _{1 1}.

  Therefore, ₁ = a₁ and ₁ = c₁/₁.

  For 2 i n - 1,

  a_i = b_{i i-1} + _i and c_i = _{i i} .

  Therefore, _i = a_i - b_{i i-1} and _i = c_i/i .

  For i = n,

  a_n = b_{n n-1} + _n .

  Therefore, _n = a_n - b_nn-1 .

  (_n = 0 because c_n = 0.)

  These steps can be executed provided _i are not zero or too close to zero!

• Step (ii): Solve Ly = d.

  y₁ = d₁/₁ and for i = 2,...,n,   y_i = (d_i - b_i y_i-1 )/_i .

• Step (iii): Solve Ux = y.

  x_n = y_n and for i = n-1,...,1,   x_i = y_i - _i x_i+1 .

Some of the loops can be combined to form the following very important algorithm.

Tridiagonal Algorithm.

   (1) = a(1)
   (1)= c(1)/(1)
   y(1) = d(1)/(1)
   for i = 2, n
        (i) = a(i)- b(i)*(i-1)
        (i) = c(i)/(i)
        y(i) = (d(i) - b(i)*y(i-1))/(i)
   endloop
   x(n) = y(n)
   for i = n - 1,1
        x(i) = y(i) - (i)*x(i+1)
   endloop.

3.2.5.   Implementation.

In this section we use Maple and the tridiagonal algorithm to solve the finite difference equations in (7) and (8). We have done several calculations with variable c in the derivative boundary condition (6). Note the Maple code has two procedures. The procedure TRID is called from within the procedure BVPF.

Maple Code for Cooling Fin.

   Input Data:
        > with(linalg):
        > n:=20:
        > T:=.1:
        > W:=10.:
        > L:=1.:
        > csur:=.1:
        > usur:= 70.:
        > uleft:= 160:
        > K:=.001:

   Procedure BVPF is Defined (components are defined):
        > BVPF:=proc(n,K,T,W,L,csur,usur,uleft)
        > a:=vector(n,0):
        > b:=vector(n,0):
        > c:=vector(n,0):
        > d:=vector(n,0):
        > h:=L/n:
        > C:=csur*2*(W + T)/(T*W):
        > for i from 1 to n-1 do
        >     a[i]:=2*K + h*h*C:
        >     c[i]:=-K:
        >     b[i]:=-K:
        >     d[i]:=h*h*C*usur:
        > od:
        > d[1]:= d[1]+K*uleft:
        > a[n]:= 2*K + h*h*C + 2*h*csur:
	> b[n]:=-2*K:
        > d[n]:= h*h*C*usur + 2*csur*usur*h:
        > TRID(n,a,b,c,d):
        > end:

   Procedure TRID is Defined (tridiagonal algorithm):
        > TRID:=proc(n,a,b,c,d)
        > alpha:=vector(n,0):
        > Gamm:=vector(n,0):
        > y:=vector(n,0):
        > alpha[1]:=a[1]:
        > Gamm[1]:= c[1]/alpha[1]:
        > y[1]:=d[1]/alpha[1]:
        > for i from 2 to n do
        >     alpha[i]:=a[i] - b[i]*Gamm[i-1]:
        >     Gamm[i]:=c[i]/alpha[i]:
        >     y[i]:=(d[i] - b[i]*y[i-1])/alpha[i]:
        > od:
        > x[n]:=y[n]:
        > for i from n-1  by -1 to 1 do
        >     x[i]:= y[i] - Gamm[i]*x[i+1]:
        > od:
        > end:

   Output of Data:
        > with(plots):

We use the smallest c = .0001, and it gives the highest temperature.

        > BVPF(20,.001,.1,10.,1.,.0001,70.,160):
        > L:=[[0,160],seq([i*h,x[i]],i=1..20)]:

The c is increased to c = .001, and the temperature decreases.

        > BVPF(20,.001,.1,10.,1.,.001,70.,160):
        > LL:=[[0,160],seq([i*h,x[i]],i=1..20)]:

The largest c is c = .01, and this gives the lowest temperature.

        > BVPF(20,.001,.1,10.,1.,.01,70.,160):
        > LLL:=[[0,160],seq([i*h,x[i]],i=1..20)]:

        > plot({L,LL,LLL},x=0..1);

Figure: Temperatures versus Space: c = .0001, .001, .01

The Matlab version is very similar and one should consult the previous section.

3.2.6.   Assessment.

In the derivation of the model for the fin we made several assumptions. If the thickness T of the fin is too large, there will be a varying temperature with the vertical coordinate. By assuming the W parameter is large, one can avoid any end effects on the temperature of the fin. Another problem arises if the temperature varies over a large range in which case the thermal conductivity K will be temperature dependent. We will return to these problems.

Once the continuum model is agreed upon and the finite difference approximation is formed, one must be concerned about the appropriate mesh size. Here an analysis much the same as in the previous section can be given. In more complicated problems several computations with decreasing mesh sizes are done until little variation in the numerical solutions is observed.

The tridiagonal algorithm is not always applicable. Difficulties will arise if either the _i are zero or near zero. The following theorem gives conditions on the components of the tridiagonal matrix so that the tridiagonal algorithm works very well.

Existence and Stability Theorem. Consider the tridiagonal algebraic system.

If

|a₁ | > |c₁ | > 0,

|a_i | > |b_i | + |c_i |, c_i 0, b_i 0 for 1 < i < n,

|a_n | > |c_n | > 0,

then

(i). 0 < |a_i | - |b_i | < |_i | < |a_i | + |b_i | for 1 i n. (This avoids division by small numbers.)

(ii). |_i | < 1 for 1 i n - 1. (This is the stability part in the backward solve loop.)

Proof. The proof uses mathematical induction on n.

Set i =1:

b₁ = 0 and |₁ | = |a₁ | > 0 and |₁ | = |c₁ | / |a₁ | < 1.

Set i > 1 and assume it is true for i-1:

_i = a_i - b_{i i-1} and _i = c_i /_i .

So, a_i = b_{i i-1} + _i and so

|ai |	|bi | |i-1 | + |i |
<	|bi |*1 + |i |.

Therefore |_i | > |a_i | - |b_i | |c_i | > 0.

Also, |i |	= |ai - bi i-1 |
	|ai | + |bi | |i-1 |
	|ai | + |bi |*1.

Next, consider _i = c_i /_i .

|_i | = |c_i | / |_i | < |c_i |/(|a_i | - |b_i |) 1.

3.2.7.   Homework.

1. Show that the tridiagonal algorithm fails for the following problem

              x1 -  x2      = 1
             -x1 + 2x2 - x3 = 1
                   -x2 + x3 = 1
   

2. Use the above code and experiment with different T,W and c. Explain your results and evaluate the accuracy of the model.

3. In the derivation of the tridiagonal algorithm we combined some of the loops. Justify this.

4. Modify the above model and code for a tapered fin where T = .2(1 - x) + .1x.

5. Find the exact solution of the fin problem and experiment with different mesh sizes. Observe the order of convergence of the discrete solution to the continuum solution.

6. Prove the analogous convergence theorem to the theorem in the previous section.