# Backward Difference Method (Fully Implicit Method)

## Transcript Of Backward Difference Method (Fully Implicit Method)

Implicit Numerical Methods

Backward Difference Method (Fully Implicit Method)

1 ∂T

m

=

1

⎡T m ⎢O

−T m−1 ⎤ O ⎥+

∆t

∂ 2T

m +"

α ∂t o α ⎢⎣ ∆t ⎥⎦ 2α ∂t 2 O

Truncation error

∈T = − ∆t ∂ 2T

m + ∆t 2 ∂ 3T

m

+

"

⎡ ⎢

2

∆

x

2

∂ 4T

m + 2∆x 4 ∂ 6T

m⎤ ⎥

2α ∂t 2 0 3!α ∂t 3 0

⎣ 4! ∂x 4 0

6! ∂x 6 0 ⎦

⎡ 2∆y 2 ∂ 4T +⎢

m + 2∆y 2 ∂ 6T

m + "⎥⎤

⎣ 4! ∂y 4 0

6! ∂x6 0

⎦

Difference scheme is first order accurate in time and the time difference is compatible with the time differential.

Stability

For a problem with no internal heat generation the difference equation is:

[ ] 1

T

m

−T

m−1

= TN m

−T m 0

+ TE m

−T m 0

+ TS m

−T m 0

Tm +W

−T m 0

α∆t 0 0

∆y 2

∆x 2

∆y 2

∆x 2

Here conditions at time step m-1 are all known and those at step m must be determined To determine stability criteria, the equation must be arranged in the form T m = ET m−1 = g

T

m−1

⎡T m = − α∆t⎢ N

+T m S

+

TE m

−T m W

⎤ ⎥

+

⎢⎡1+ 2α∆t[

1

0

⎣ ∆y2

∆x2 ⎦ ⎣

∆x 2

T m−1 = HT m − g

+

1

⎤m ]⎥ T0

∆y2 ⎦

or T m = H −1T m−1 + H −1 g

Stability requires

λSR H −1 ≤1 or λSR H ≥ 1

Look at node O

− α∆tT m N

− α∆t

TE m

− α∆t

TS m

− α∆t

Tm W

+

⎢⎡1+ 2α∆t[

1

+

1

⎤m ]⎥ TO

=

T m−1 O

∆y 2

∆x 2

∆y 2

∆x 2

⎣

∆x2 ∆y 2 ⎦

The left hand side of the above equation provides the coefficients in the H matrix. Apply Gerschgorin's Theorem.

λ − (1 + 2α∆t[ 1 + 1 ]) ≤ − α∆t + − α∆t + − α∆t + − α∆t

∆x2 ∆y 2

∆y 2

∆x2 ∆y 2

∆x 2

≤

2α∆t

⎡ ⎢

1

+

1

⎤ ⎥

⎣ ∆x 2 ∆y 2 ⎦

The interval containing the eigenvalues is:

← 2α∆t⎜⎜⎛

1

+

1

⎞ ⎟⎟

→

← 2α∆t⎜⎜⎛

1

+

1

⎞ ⎟⎟

→

⎝ ∆x2 ∆y 2 ⎠

⎝ ∆x2 ∆y 2 ⎠

λ 2

1+ 2∆tα ⎢⎡

1

+

1

⎤ ⎥

λ 1

⎣ ∆x2 ∆y2 ⎦

The upper bound of this interval is:

λ1 =1+ 2∆tα ⎢⎡

1

+

1

⎥⎤ + 2∆tα ⎢⎡

1

+

1

⎤ ⎥

⎣ ∆x2 ∆y 2 ⎦

⎣ ∆x2 ∆y 2 ⎦

⎡ =1+ 4∆tα ⎢

1

+

1

⎤ ⎥

⎣ ∆x 2 ∆y 2 ⎦

The lower bound of the interval is:

Which is always greater than 1.

λ 2

=1+

2∆tα

⎡ ⎢

1

+

1 ⎥⎤ − 2∆tα ⎢⎡ 1

+

1⎤ ⎥ =1

⎣ ∆x2 ∆y 2 ⎦

⎣ ∆x2 ∆y2 ⎦

So the eigenvalues of the H matrix associated with the interior equations are always greater

than or equal to one. This means that the eigenvalues associated with the inverse of H are

always less than or equal to one, producing unconditional stability.

Summary of the Backward Difference Method

1. Set of equations are unconditionally stable.

2. Computational time per time step will be longer than that for the forward difference since the method is implicit, i.e. the set of finite difference equations must be solved simultaneously at each time step.

3. The influence of a perturbation is felt immediately throughout the complete region.

Crank-Nicolson Method

Crank-Nicolson splits the difference between Forward and Backward difference schemes. In terms of the equations used to introduce transient conduction methods, the time weighting factor f is 0.5. The equation is evaluated halfway between the old (m) and new (m+1) time levels.

1 ∂T m+1/ 2 ∂ 2T m+1/ 2 ∂ 2T m+1/ 2

α ∂t i, j

= ∂x 2 i, j

+ ∂y 2 i, j

Use Taylor’s series for time derivative

m+1

m+1/ 2 ∆t ∂T m+1/ 2 ∆t 2 ∂ 2T m+1/ 2 ∆t 3 ∂ 3T m+1/ 2

Ti, j = Ti, j

+ 2 ∂t

+ 222! ∂t 2

+ 233! ∂t 3

i, j

i, j

i, j

m

m+1/ 2 ∆t ∂T m+1/ 2 ∆t 2 ∂ 2T m+1/ 2 ∆t 3 ∂ 3T m+1/ 2

Ti, j = Ti, j

− 2 ∂t

+ 222! ∂t 2

− 233! ∂t 3

i, j

i, j

i, j

m

m+1/ 2 ∆t ∂T m+1/ 2 ∆t 2 ∂ 2T m+1/ 2 ∆t 3 ∂ 3T m+1/ 2

Ti, j = Ti, j

+ 2 ∂t i, j

+ 22 2! ∂t 2 i, j

+ 233! ∂t 3 i, j

Subtract the second from the first, and rearrange

∂T

m+1/ 2

T m+1 − T m

i, j

i, j

2∆t 2 ∂ 2T

m+1/ 2

∂t i, j =

∆t

− 22 2! ∂t 2 i, j

error 0 (∆t 2 )

m+1/ 2

Spatial derivatives (central difference about i, j )

∂ 2T m+1/ 2

∂ 2T m+1/ 2

∂ 2T m+1 = ∂ 2T m+1/ 2 + ∆t ∂( ∂x 2 )

+ ∆t 2 ∂( ∂x2 )

∂x 2

∂x 2

2 ∂t

222! ∂t 2

i, j

i, j

i, j

i, j

∂ 2T

∂ 2T

∂ 2T m ∂ 2T m+1/ 2 ∆t ∂( ∂x 2 ) m+1/ 2 ∆t 2 ∂( ∂x 2 ) m+1/ 2

∂x 2 i, j = ∂x 2 i, j

+ 2 ∂t i, j

+ 222! ∂t 2 i, j

∂ 2T m+1/ 2

∂ 2T m+1/ 2

∂ 2T m = ∂ 2T m+1/ 2 − ∆t ∂( ∂x 2 )

+ ∆t 2 ∂( ∂x 2 )

∂x 2

∂x 2

2 ∂t

222! ∂t 2

i, j

i, j

i, j

i, j

Add and divide by 2

∂ 2T

∂ 2T m+1/ 2 1 ⎡ ∂ 2T m+1 ∂ 2T m ⎤ ∆t 2 ∂( ∂x 2 ) m+1/ 2

∂x 2 i, j

=

2

⎢ ⎣

∂x

2

i, j

+ ∂x 2

i, j

⎥− 2 ⎦ 2 2!

∂t 2

i, j

So the spatial terms are also second order accurate in time. ∂ 2T m+1/ 2

The same treatment holds for ∂y 2 i, j

In compass notation the full Crank-Nicholson equation is:

1

⎡T m+1 ⎢O

−T m O

⎤ ⎥

=

1

⎢⎡TN m+1

− TS m+1

+

T m+1 E

− T m+1 W

−

2

⎡ ⎢

1

+

1

⎥⎤ TO m+1 ⎥⎤

α ⎢⎣ ∆t ⎥⎦ 2 ⎣ ∆y2

∆x 2

⎣ ∆x2 ∆x2 ⎦

⎦

+

1

⎡T m ⎢N

+ TS m

+ TE m

+T m W

⎡ − 2⎢

1

+

1

⎤ m⎤ ⎥T ⎥

2 ⎢⎣ ∆y 2

∆x 2

⎣ ∆x 2 ∆y 2 ⎦ O ⎥⎦

Error is O (∆t 2 , ∆x 2 , ∆y 2 )

Stability of the One-dimensional Crank-Nicholson Method

Continue to assume uniform spacing in the spatial mesh

1

⎡T m+1 ⎢i

−T m i

⎤ ⎥=

1

⎡ Ti +m1+1 ⎢

+ Ti−1m+1

− 2T m+1 ⎤ i ⎥+

1

⎡ Ti +1m ⎢

+ Ti−1m

− 2T m i

⎤ ⎥

α ⎣ ∆t ⎦ 2 ⎣

∆x 2

⎦ 2⎣

∆x 2

⎦

Unknowns Tm+1

knowns Tm

153

Multiply equation by 2α∆t

[ ] 2 T m+1 −T m

=α∆t

⎡⎛ ⎢⎜

Ti

+1

m

+1

+ Ti−1m+1

−

2T m+1 i

⎟⎞ + ⎜⎛ Ti+1m

+ Ti−1m

−

2T m i

⎞⎤ ⎟⎥

i

i

⎢⎣⎜⎝

∆x 2

⎟⎠ ⎜⎝

∆x 2

⎟⎠⎥⎦

Arrange into the form ATm+1 = BTm + d

−

α∆t

⎡ ⎢

Ti+1m+1+

T m+1 i −1

⎤ ⎥

+

2

⎢⎡1

+

α∆t

⎤ ⎥

Ti

m+1

=

α∆t

⎡ ⎢

Ti+1

+

Tm i−1

⎤ ⎥

+

2⎢⎡1

−

α∆t

⎤ ⎥

Ti

m

⎣ ∆x2 ⎦ ⎣ ∆x2 ⎦

⎣ ∆x2 ⎦ ⎣ ∆x2 ⎦

Pick matrix H such that h = − 2α∆t and h = h = α∆t

i ,i

∆x 2

i,i+1 i,i+1 ∆x 2

then A = (2I-H) and B = (2I + H) = (4I-A)

T m+1 = A−1B T m + A−1d = A−1(4I − A) T m + A−1d = (4IA−1 −1) T m + A−1d = ( 4I − 1) T m + A−1d A

Stability requires that ⎜⎛ 4l −1⎟⎞ have a spectral radius equal or less than 1. ⎝A ⎠

4 − 1 ≤1 for any eigenvalue of A λ A

−1≤ 4 −1≤1 λ A

− λA ≤ 4 − λA ≤ λA 0 ≤ 4 ≤ 2λA

or λA ≥ 2

Now we're back to Gerschgorin's Theorem to finish the job. Summing the absolute values of the off-diagonal elements of A gives:

λ A

− 2 ⎢⎡1 +

α∆t ⎤ 2⎥

≤

2α∆t

2

⎣ ∆x ⎦ ∆x

The real interval containing these eigenvalues is:

← 2α∆t → ← 2α∆t →

∆x 2

∆x 2

λ2

2⎢⎡1+ α∆2t ⎥⎤

λ1

⎣ ∆x ⎦

so

λ = 2 + 4 α∆t and λ = 2

1

∆x 2

2

Therefore λA ≥ 2 , always and the Crank-Nicholson finite difference equations are unconditionally stable

Transformation of time

Implicit Methods lend themselves to analyzing systems, which go through some initial transient then settle into some relatively steady behavior. Smaller time steps can be used to resolve the initial behavior and very large ones used later in time. In fact, for complex, non-linear system, the steady state is obtained by simply running a transient for a sufficiently long period of time, or an iteration to obtain a steady solution may be built around use of pseudo-transient equations.

One useful method for situations known to have a steady end state involves a transformation of the time variable. Introduce new variable

γ =1− e− At

0≤t ≤∞ 0≤ γ ≤1

The time derivative transforms to:

∂T = ∂T dγ = Ae−At ∂T = A(1 − γ ) ∂T

∂t ∂γ dt

∂γ

∂γ

A conduction equation transforms to:

A(1 − γ ) ∂T = ∂ 2T + ∂ 2T α ∂γ ∂x2 ∂y 2

The backward finite difference equation in the transformed variable (at an interior mesh point) is:

T m +T m T m −T m ⎡

N

S +E

W − 2⎢

1

+

1

⎤m ⎥ TO

=

A[1 − γ

] [T m O

− T m−1 ] O

∆y 2

∆x 2

⎣ ∆x2 ∆y 2 ⎦

α∆γ

Backward Difference Method (Fully Implicit Method)

1 ∂T

m

=

1

⎡T m ⎢O

−T m−1 ⎤ O ⎥+

∆t

∂ 2T

m +"

α ∂t o α ⎢⎣ ∆t ⎥⎦ 2α ∂t 2 O

Truncation error

∈T = − ∆t ∂ 2T

m + ∆t 2 ∂ 3T

m

+

"

⎡ ⎢

2

∆

x

2

∂ 4T

m + 2∆x 4 ∂ 6T

m⎤ ⎥

2α ∂t 2 0 3!α ∂t 3 0

⎣ 4! ∂x 4 0

6! ∂x 6 0 ⎦

⎡ 2∆y 2 ∂ 4T +⎢

m + 2∆y 2 ∂ 6T

m + "⎥⎤

⎣ 4! ∂y 4 0

6! ∂x6 0

⎦

Difference scheme is first order accurate in time and the time difference is compatible with the time differential.

Stability

For a problem with no internal heat generation the difference equation is:

[ ] 1

T

m

−T

m−1

= TN m

−T m 0

+ TE m

−T m 0

+ TS m

−T m 0

Tm +W

−T m 0

α∆t 0 0

∆y 2

∆x 2

∆y 2

∆x 2

Here conditions at time step m-1 are all known and those at step m must be determined To determine stability criteria, the equation must be arranged in the form T m = ET m−1 = g

T

m−1

⎡T m = − α∆t⎢ N

+T m S

+

TE m

−T m W

⎤ ⎥

+

⎢⎡1+ 2α∆t[

1

0

⎣ ∆y2

∆x2 ⎦ ⎣

∆x 2

T m−1 = HT m − g

+

1

⎤m ]⎥ T0

∆y2 ⎦

or T m = H −1T m−1 + H −1 g

Stability requires

λSR H −1 ≤1 or λSR H ≥ 1

Look at node O

− α∆tT m N

− α∆t

TE m

− α∆t

TS m

− α∆t

Tm W

+

⎢⎡1+ 2α∆t[

1

+

1

⎤m ]⎥ TO

=

T m−1 O

∆y 2

∆x 2

∆y 2

∆x 2

⎣

∆x2 ∆y 2 ⎦

The left hand side of the above equation provides the coefficients in the H matrix. Apply Gerschgorin's Theorem.

λ − (1 + 2α∆t[ 1 + 1 ]) ≤ − α∆t + − α∆t + − α∆t + − α∆t

∆x2 ∆y 2

∆y 2

∆x2 ∆y 2

∆x 2

≤

2α∆t

⎡ ⎢

1

+

1

⎤ ⎥

⎣ ∆x 2 ∆y 2 ⎦

The interval containing the eigenvalues is:

← 2α∆t⎜⎜⎛

1

+

1

⎞ ⎟⎟

→

← 2α∆t⎜⎜⎛

1

+

1

⎞ ⎟⎟

→

⎝ ∆x2 ∆y 2 ⎠

⎝ ∆x2 ∆y 2 ⎠

λ 2

1+ 2∆tα ⎢⎡

1

+

1

⎤ ⎥

λ 1

⎣ ∆x2 ∆y2 ⎦

The upper bound of this interval is:

λ1 =1+ 2∆tα ⎢⎡

1

+

1

⎥⎤ + 2∆tα ⎢⎡

1

+

1

⎤ ⎥

⎣ ∆x2 ∆y 2 ⎦

⎣ ∆x2 ∆y 2 ⎦

⎡ =1+ 4∆tα ⎢

1

+

1

⎤ ⎥

⎣ ∆x 2 ∆y 2 ⎦

The lower bound of the interval is:

Which is always greater than 1.

λ 2

=1+

2∆tα

⎡ ⎢

1

+

1 ⎥⎤ − 2∆tα ⎢⎡ 1

+

1⎤ ⎥ =1

⎣ ∆x2 ∆y 2 ⎦

⎣ ∆x2 ∆y2 ⎦

So the eigenvalues of the H matrix associated with the interior equations are always greater

than or equal to one. This means that the eigenvalues associated with the inverse of H are

always less than or equal to one, producing unconditional stability.

Summary of the Backward Difference Method

1. Set of equations are unconditionally stable.

2. Computational time per time step will be longer than that for the forward difference since the method is implicit, i.e. the set of finite difference equations must be solved simultaneously at each time step.

3. The influence of a perturbation is felt immediately throughout the complete region.

Crank-Nicolson Method

Crank-Nicolson splits the difference between Forward and Backward difference schemes. In terms of the equations used to introduce transient conduction methods, the time weighting factor f is 0.5. The equation is evaluated halfway between the old (m) and new (m+1) time levels.

1 ∂T m+1/ 2 ∂ 2T m+1/ 2 ∂ 2T m+1/ 2

α ∂t i, j

= ∂x 2 i, j

+ ∂y 2 i, j

Use Taylor’s series for time derivative

m+1

m+1/ 2 ∆t ∂T m+1/ 2 ∆t 2 ∂ 2T m+1/ 2 ∆t 3 ∂ 3T m+1/ 2

Ti, j = Ti, j

+ 2 ∂t

+ 222! ∂t 2

+ 233! ∂t 3

i, j

i, j

i, j

m

m+1/ 2 ∆t ∂T m+1/ 2 ∆t 2 ∂ 2T m+1/ 2 ∆t 3 ∂ 3T m+1/ 2

Ti, j = Ti, j

− 2 ∂t

+ 222! ∂t 2

− 233! ∂t 3

i, j

i, j

i, j

m

m+1/ 2 ∆t ∂T m+1/ 2 ∆t 2 ∂ 2T m+1/ 2 ∆t 3 ∂ 3T m+1/ 2

Ti, j = Ti, j

+ 2 ∂t i, j

+ 22 2! ∂t 2 i, j

+ 233! ∂t 3 i, j

Subtract the second from the first, and rearrange

∂T

m+1/ 2

T m+1 − T m

i, j

i, j

2∆t 2 ∂ 2T

m+1/ 2

∂t i, j =

∆t

− 22 2! ∂t 2 i, j

error 0 (∆t 2 )

m+1/ 2

Spatial derivatives (central difference about i, j )

∂ 2T m+1/ 2

∂ 2T m+1/ 2

∂ 2T m+1 = ∂ 2T m+1/ 2 + ∆t ∂( ∂x 2 )

+ ∆t 2 ∂( ∂x2 )

∂x 2

∂x 2

2 ∂t

222! ∂t 2

i, j

i, j

i, j

i, j

∂ 2T

∂ 2T

∂ 2T m ∂ 2T m+1/ 2 ∆t ∂( ∂x 2 ) m+1/ 2 ∆t 2 ∂( ∂x 2 ) m+1/ 2

∂x 2 i, j = ∂x 2 i, j

+ 2 ∂t i, j

+ 222! ∂t 2 i, j

∂ 2T m+1/ 2

∂ 2T m+1/ 2

∂ 2T m = ∂ 2T m+1/ 2 − ∆t ∂( ∂x 2 )

+ ∆t 2 ∂( ∂x 2 )

∂x 2

∂x 2

2 ∂t

222! ∂t 2

i, j

i, j

i, j

i, j

Add and divide by 2

∂ 2T

∂ 2T m+1/ 2 1 ⎡ ∂ 2T m+1 ∂ 2T m ⎤ ∆t 2 ∂( ∂x 2 ) m+1/ 2

∂x 2 i, j

=

2

⎢ ⎣

∂x

2

i, j

+ ∂x 2

i, j

⎥− 2 ⎦ 2 2!

∂t 2

i, j

So the spatial terms are also second order accurate in time. ∂ 2T m+1/ 2

The same treatment holds for ∂y 2 i, j

In compass notation the full Crank-Nicholson equation is:

1

⎡T m+1 ⎢O

−T m O

⎤ ⎥

=

1

⎢⎡TN m+1

− TS m+1

+

T m+1 E

− T m+1 W

−

2

⎡ ⎢

1

+

1

⎥⎤ TO m+1 ⎥⎤

α ⎢⎣ ∆t ⎥⎦ 2 ⎣ ∆y2

∆x 2

⎣ ∆x2 ∆x2 ⎦

⎦

+

1

⎡T m ⎢N

+ TS m

+ TE m

+T m W

⎡ − 2⎢

1

+

1

⎤ m⎤ ⎥T ⎥

2 ⎢⎣ ∆y 2

∆x 2

⎣ ∆x 2 ∆y 2 ⎦ O ⎥⎦

Error is O (∆t 2 , ∆x 2 , ∆y 2 )

Stability of the One-dimensional Crank-Nicholson Method

Continue to assume uniform spacing in the spatial mesh

1

⎡T m+1 ⎢i

−T m i

⎤ ⎥=

1

⎡ Ti +m1+1 ⎢

+ Ti−1m+1

− 2T m+1 ⎤ i ⎥+

1

⎡ Ti +1m ⎢

+ Ti−1m

− 2T m i

⎤ ⎥

α ⎣ ∆t ⎦ 2 ⎣

∆x 2

⎦ 2⎣

∆x 2

⎦

Unknowns Tm+1

knowns Tm

153

Multiply equation by 2α∆t

[ ] 2 T m+1 −T m

=α∆t

⎡⎛ ⎢⎜

Ti

+1

m

+1

+ Ti−1m+1

−

2T m+1 i

⎟⎞ + ⎜⎛ Ti+1m

+ Ti−1m

−

2T m i

⎞⎤ ⎟⎥

i

i

⎢⎣⎜⎝

∆x 2

⎟⎠ ⎜⎝

∆x 2

⎟⎠⎥⎦

Arrange into the form ATm+1 = BTm + d

−

α∆t

⎡ ⎢

Ti+1m+1+

T m+1 i −1

⎤ ⎥

+

2

⎢⎡1

+

α∆t

⎤ ⎥

Ti

m+1

=

α∆t

⎡ ⎢

Ti+1

+

Tm i−1

⎤ ⎥

+

2⎢⎡1

−

α∆t

⎤ ⎥

Ti

m

⎣ ∆x2 ⎦ ⎣ ∆x2 ⎦

⎣ ∆x2 ⎦ ⎣ ∆x2 ⎦

Pick matrix H such that h = − 2α∆t and h = h = α∆t

i ,i

∆x 2

i,i+1 i,i+1 ∆x 2

then A = (2I-H) and B = (2I + H) = (4I-A)

T m+1 = A−1B T m + A−1d = A−1(4I − A) T m + A−1d = (4IA−1 −1) T m + A−1d = ( 4I − 1) T m + A−1d A

Stability requires that ⎜⎛ 4l −1⎟⎞ have a spectral radius equal or less than 1. ⎝A ⎠

4 − 1 ≤1 for any eigenvalue of A λ A

−1≤ 4 −1≤1 λ A

− λA ≤ 4 − λA ≤ λA 0 ≤ 4 ≤ 2λA

or λA ≥ 2

Now we're back to Gerschgorin's Theorem to finish the job. Summing the absolute values of the off-diagonal elements of A gives:

λ A

− 2 ⎢⎡1 +

α∆t ⎤ 2⎥

≤

2α∆t

2

⎣ ∆x ⎦ ∆x

The real interval containing these eigenvalues is:

← 2α∆t → ← 2α∆t →

∆x 2

∆x 2

λ2

2⎢⎡1+ α∆2t ⎥⎤

λ1

⎣ ∆x ⎦

so

λ = 2 + 4 α∆t and λ = 2

1

∆x 2

2

Therefore λA ≥ 2 , always and the Crank-Nicholson finite difference equations are unconditionally stable

Transformation of time

Implicit Methods lend themselves to analyzing systems, which go through some initial transient then settle into some relatively steady behavior. Smaller time steps can be used to resolve the initial behavior and very large ones used later in time. In fact, for complex, non-linear system, the steady state is obtained by simply running a transient for a sufficiently long period of time, or an iteration to obtain a steady solution may be built around use of pseudo-transient equations.

One useful method for situations known to have a steady end state involves a transformation of the time variable. Introduce new variable

γ =1− e− At

0≤t ≤∞ 0≤ γ ≤1

The time derivative transforms to:

∂T = ∂T dγ = Ae−At ∂T = A(1 − γ ) ∂T

∂t ∂γ dt

∂γ

∂γ

A conduction equation transforms to:

A(1 − γ ) ∂T = ∂ 2T + ∂ 2T α ∂γ ∂x2 ∂y 2

The backward finite difference equation in the transformed variable (at an interior mesh point) is:

T m +T m T m −T m ⎡

N

S +E

W − 2⎢

1

+

1

⎤m ⎥ TO

=

A[1 − γ

] [T m O

− T m−1 ] O

∆y 2

∆x 2

⎣ ∆x2 ∆y 2 ⎦

α∆γ