Implementation of Alternating Direction Implicit Method to Solve Two Dimensional Heat Conduction with Mathematica

Scott Beckman
August 2000

University of California at Berkeley
Material Science and Engineering Department
577 Evans Hall #1760
Berkeley, CA 94720
sbeckman@uclink.berkeley.edu

Introduction  1

Alternating Direction Method  1.1

The alternating direction implicit, ADI, method is a common algorithm used in finite difference methods.  Below it is demonstrated how to implement ADI with Mathematica.  Although Mathematica isn't the most efficient package for finite difference problems, it is an excellent tool for visualizing how ADI works.  Large scale application should inevitably involve translating the below into C or FORTRAN code/

Parabolic partial differential equation, PDE, can be expressed at implicit finite difference equations.  When such a PDE is has two spatial dimensions, the number of terms in the equation are too large to be solved efficiently.  The basic gist of ADI is to break a large implicit finite difference equation, over time step δt, into two step, each over the time step .  The first step implicitly solves for heat conduction in one direction only.  The input to this step are the current temperatures T(x,y,n) and the results from this step are an intermediate temperature distribution, called , which is sometimes incorrectly referred to as T(x,y,n+).  The second step solves implicitly for the heat conduction in the other directional only.  The input to this step is, and the resulting output is the temperature at n+1, T(x,y,n+1).

So thermal conduction in 2D

expressed as a finite difference equation

can be broken down to

The former of the above two is first solved for and the result is substituted in to the latter equation to find T(x,y,n+1).  Notice that upon arranging the first equation to find the expression is an implicit finite different equation in the x direction only.  The second equation, rearranged is an implicit equation in the y direction only.

Details of the ADI method, such as error and proof of stability, can be found in Carnahan's text.

    Carnahan, Luther,and Wilkes.  Applied Numerical Methods. New York: John Wiley & Sons, Inc., 1969.

Difference Equations and Boundary Conditions  1.2

As shown above, the first derivative is written

Similarly second derivatives are written

Boundary conditions require a special formulation of the spatial derivative.  For our problem we will be using Newton's law of cooling as our boundary condition

which can be arranged to the form

Carnahan shows that the second derivative of this boundary condition, which he calls, "third condition," is written

Solving this problem is now a matter of substituting the appropriate approximations for the derivatives and solving the resulting matrices, which will be tridiagonal.  

Conventions and Initiation  2

To reduce the number of constants that are carried throughout this problem, the following grouping conventions will be used.

The below commands turn off the error code [General::spell] which is a common warning given in this notebook because of the naming scheme selected.  For debugging purposes this can be turned back on with the command On[General::spell]

The below add-on Mathematica add-on packages are needed for this problem.

We define the number of x step Nx and number of y step Ny so (1≤x≤Nx) and (1≤y≤Ny)

Set the temperature, T to be the function Tfct.  Note that the hold operator is being used only for the convenience of demonstration.  The actual electrothermal does not need such.  

Implicit in X Direction  3

Prepare Input to TridiagonalSolve  3.1

Mathematica has a function, TridiagonalSolve which will be used.  For an exact explanation of this function, the reader is referred to the Mathematica on-line use guide.  In short the TridiagonalSolve function operates as

T=TridiagonalSolve[A,B,C,D]

where T,A,B,C, and D are vectors

For the ADI method, solving for A,B,C, and D vectors is a matter of substituting equations and rearranging.  In this example, the solution for each instance will not be given.  As it will be shown below, it is possible for the reader to use this notebook to back out the system of equations for each point in space.  In the accompanying notebook, Implementation of Brian's Method to Solve Three Dimensional Heat Conduction with Mathematica, it is shown how to use Mathematica to derive the vectors for input to the TridiagonalSolve function.  

Check the Input  3.2

We can now extract columns from the matrices to check that the solution we entered matches the solution achieved when solving by hand.

For the case middle of the region, (1<i<Nx) and (1<j<Ny), take values, say

Which, when substituted, should match the equation

In this manner, all of the systems of equations can be verified.

Assign Values to Constants  3.3

To achieve a numerical answer from the TridiagonalSolve function, values must be assigned to the constants

Loop Over Entire Region  3.4

Implicit in Y Direction  4

Prepare Input to Tridiagonal Solve  4.1

Clearing values so the correctness of step two can be verified symbolically.

Writing the equations for the second step ones notices that there is a direct parallelism between the first and second step.  Take for instance the equation which is comparable to the sample given in section 3.2

This suggests that the vector inputs to the TridiagonalSolve function can be created by manipulating the existing Xstar matrices.

For Astar, Bstar, and Cstar this is a matter of substituting ϕy for ϕx adding a rows, removing a columns, and then tranposing

Dstar, unfortunately isn't as easy to rearrange because the function's conditionals also change.  Rather than transform the Dstar matrix, we'll just rewrite the function exchanging x for y and correcting the order direction of the T values from the y direction to the x direction.  

Check the Input  4.2

We can now extract columns from the matrices to check that the solution we entered matches the solution achieved when solving by hand.

For the case middle of the region, (1<i<Nx) and (1<j<Ny), take values, say

Which, when substituted, should match the equation

In this manner,all of the systems of equations can be verified.

Assign Values to Constants  4.3

Loop Over Entire Region  4.4

Multiple Loops and Results  5

Although the particular nature of the engineering problem designates the values used for the constants and the exact method which the ADI is employed, below is a brief example of using a loop to solve for the evolution of the temperature distribution with time.