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 ⎦
α∆γ