Alternating Direction Implicit Methods (ADI)

The fully implicit method involves solving a large number, (M − 1) × (N − 1), of simultaneous equations each time-step. An important alternative is the ADI method, whose basic idea is to be implicit in only one of the x, y directions at a time, and to alternate between the two. In a complete cycle, it is only necessary to solve (M − 1) and then (N − 1) equations, which is considerably more efficient. Furthermore, the systems of equations are tri-diagonal.

During the first half time-step we are implicit in x, say, so that

\frac{U_{mn}^{j+1/2} − U_{mn}^{j}}{k/2} = \frac{1}{h^{2}} (δ_{x}^{2}U_{mn}^{j+1/2} + δ_{y}^{2}U_{mn}^{j})

or

(1 − \frac{1}{2}rδ_{x}^{2}) U_{mn}^{j+1/2} = (1 + \frac{1}{2}rδ_{y}^{2})U_{mn}^{j} .  (8.4)

Over the next half time-step we treat y implicitly, so that

\frac{U_{mn}^{j+1} − U_{mn}^{j+1/2}}{k/2} = \frac{1}{h^{2}} (δ_{x}^{2}U_{mn}^{j+1/2} + δ_{y}^{2}U_{mn}^{j+1})

or

(1 − \frac{1}{2}rδ_{y}^{2}) U_{mn}^{j+1} = (1 + \frac{1}{2}rδ_{x}^{2})U_{mn}^{j+1/2} .  (8.5)

If we want, we can eliminate the values at the half-time level by operating on (8.4) with (1 + \frac{1}{2}rδ_{x}^{2}) and (8.5) with (1 − \frac{1}{2}rδ_{x}^{2}) to obtain

(1 − \frac{1}{2}rδ_{x}^{2}) (1 − \frac{1}{2}rδ_{y}^{2}) U_{mn}^{j+1} = (1 + \frac{1}{2}rδ_{x}^{2})(1 + \frac{1}{2}rδ_{y}^{2})U_{mn}^{j} .  (8.6)

We now consider stability of the scheme to a double Fourier perturbation of the form U_{mn}^{j} = Û^{j} exp(imξ + inη). In the familiar manner, the operator δ_{x}^{2} becomes −4 sin^{2} \frac{1}{2}ξ, while δ_{y}^{2} becomes −4 sin^{2} \frac{1}{2}η. Then (8.4) and (8.5) imply

\begin{matrix}Û^{j+1/2} & = \frac{1 − 2r sin^{2} \frac{1}{2}η}{1 + 2r sin^{2} \frac{1}{2}ξ} Û^{j} = λ_{1}Û^{j} \\ Û^{j+1} & = \frac{1 − 2r sin^{2} \frac{1}{2}ξ}{1 + 2r sin^{2} \frac{1}{2}η} Û^{j+1/2} = λ_{2}Û^{j+1/2} \end{matrix} \}  say.

Then |λ_{1}| ⩽ 1 for all ξ, η only if r ⩽ 1, which is equivalent to “r ⩽ \frac{1}{2}” for a time-step \frac{1}{2}k rather than k. λ_{2} is similar. However, over the complete time-step, Û^{j+1} = λ_{1}λ_{2}Û^{j} and

|λ_{1}λ_{2}| = |\frac{1 − 2r sin^{2} \frac{1}{2}ξ}{1 + 2r sin^{2} \frac{1}{2}ξ} | |\frac{1 − 2r sin^{2} \frac{1}{2}η}{1 + 2r sin^{2} \frac{1}{2}η} | ⩽ 1  (8.7)

for every ξ and η. Thus, although each half time-step is unstable if repeated indefinitely, the alternation of one and then the other is unconditionally stable!

ADI in Three Dimensions

The obvious generalisation to three dimensions would be

\begin{matrix}(1 − \frac{1}{3}rδ_{x}^{2}) U_{lmn}^{j+1/3} & = (1 + \frac{1}{3}r(δ_{y}^{2} + δ_{z}^{2})) U_{lmn}^{j} \\ (1 − \frac{1}{3}rδ_{y}^{2}) U_{lmn}^{j+2/3} & = (1 + \frac{1}{3}r(δ_{x}^{2} + δ_{z}^{2})) U_{lmn}^{j+1/3} \\ (1 − \frac{1}{3}rδ_{z}^{2}) U_{lmn}^{j+1} & = (1 + \frac{1}{3}r(δ_{x}^{2} + δ_{y}^{2})) U_{lmn}^{j+2/3} \end{matrix} \}  (8.8)

However, such a scheme turns out not to be unconditionally stable, and so is little improvement on explicit schemes. If we consider U_{lmn}^{j} = Û^{j} exp(ilξ + imη + inτ), we get

Û^{j+1/3} = λ_{1}Û^{j} etc. where λ_{1} = \frac{1 − \frac{4}{3}r(sin^{2} \frac{1}{2}η + sin^{2} \frac{1}{2}τ)}{1 + \frac{4}{3}r sin^{2} \frac{1}{2}ξ} = \frac{1 − b − c}{1 + a} ,

Then

|λ_{1}λ_{2}λ_{3}| = |\frac{(1 − b − c)(1 − a − b)(1 − a − c)}{(1 + a)(1 + b)(1 + c)}| .  (8.9)

Stability occurs only if a, b, c are small enough, with the worst case being ξ = η = τ = π, when the necessary condition is r ⩽ 3/4.

A more successful generalisation is obtained by using a Crank-Nicolson idea (θ = \frac{1}{2}) in the implicit direction. Consider the 2-D scheme

\begin{matrix}U_{mn}^{∗j+1} − U_{mn}^{j} & = \frac{1}{2}rδ_{x}^{2}(U_{mn}^{∗j+1} + U_{mn}^{j}) + rδ_{y}^{2}U_{mn}^{j} \\ U_{mn}^{j+1} − U_{mn}^{j} & = \frac{1}{2}rδ_{x}^{2}(U_{mn}^{∗j+1} + U_{mn}^{j}) + \frac{1}{2}rδ_{y}^{2}(U_{mn}^{j+1} + U_{mn}^{j}) \end{matrix} \} .  (8.10)

If we eliminate U^{∗} from (8.10) we obtain the same relation (8.6) between U^{j+1} and U^{j}. However, whereas the intermediate value U^{j+1/2} in (8.4) is an approximation to the real solution at the j + \frac{1}{2} time-level, the intermediate values in (8.10) should be regarded as a first estimate of the solution at j + 1, which is subsequently improved upon.

In three dimensions, (8.10) generalises to the scheme with two estimates U^{∗} and U_{∗}^{∗},

\begin{matrix}(1 − \frac{1}{2}rδ_{x}^{2})U_{lmn}^{∗j+1} & = [1 + r(\frac{1}{2}δ_{x}^{2} + δ_{y}^{2} + δ_{z}^{2})] U _{lmn}^{j} \\ (1 − \frac{1}{2}rδ_{y}^{2})U_{∗lmn}^{∗j+1} & = [1 + r(\frac{1}{2}δ_{x}^{2} + \frac{1}{2}δ_{y}^{2} + δ_{z}^{2})] U _{lmn}^{j} + \frac{1}{2}rδ_{x}^{2}U_{lmn}^{∗j+1} \\ (1 − \frac{1}{2}rδ_{z}^{2})U _{lmn}^{j+1} & = [1 + \frac{1}{2}r(δ_{x}^{2} + δ_{y}^{2} + δ_{z}^{2})] U _{lmn}^{j} + \frac{1}{2}rδ_{x}^{2}U_{lmn}^{∗j+1} + \frac{1}{2}rδ_{y}^{2}U_{∗lmn}^{∗j+1} \end{matrix} \}

(8.11)

This scheme is unconditionally stable, as shown on example sheet 2.