alternating-direction implicit method

The ADI method is a two step iteration process that alternately updates the column and row spaces of an approximate solution to $AX\; -\; XB\; =\; C$. One ADI iteration consists of the following steps:{{Cite journal|last=Wachspress|first=Eugene L.|date=2008|title=Trail to a Lyapunov equation solver|journal=Computers & Mathematics with Applications|volume=55|issue=8|pages=1653–1659|doi=10.1016/j.camwa.2007.04.048|issn=0898-1221|doi-access=free}}

1. Solve for $X^\{(j\; +\; 1/2)\}$, where $\backslash left(\; A\; -\; \backslash beta\_\{j\; +1\}\; I\backslash right)\; X^\{(j+1/2)\}\; =\; X^\{(j)\}\backslash left(\; B\; -\; \backslash beta\_\{j\; +\; 1\}\; I\; \backslash right)\; +\; C.$

2. Solve for $X^\{(j\; +\; 1)\}$, where $X^\{(j+1)\}\backslash left(\; B\; -\; \backslash alpha\_\{j\; +\; 1\}\; I\; \backslash right)\; =\; \backslash left(\; A\; -\; \backslash alpha\_\{j+1\}\; I\backslash right)\; X^\{(j+1/2)\}\; -\; C$.

The numbers $(\backslash alpha\_\{j+1\},\; \backslash beta\_\{j+1\})$ are called shift parameters, and convergence depends strongly on the choice of these parameters.{{Cite journal|last1=Lu|first1=An|last2=Wachspress|first2=E.L.|date=1991|title=Solution of Lyapunov equations by alternating direction implicit iteration|journal=Computers & Mathematics with Applications|volume=21|issue=9|pages=43–58|doi=10.1016/0898-1221(91)90124-m|issn=0898-1221|doi-access=}}{{Cite journal|last1=Beckermann|first1=Bernhard|last2=Townsend|first2=Alex|date=2017|title=On the Singular Values of Matrices with Displacement Structure|journal=SIAM Journal on Matrix Analysis and Applications|language=en|volume=38|issue=4|pages=1227–1248|doi=10.1137/16m1096426|issn=0895-4798|arxiv=1609.09494|s2cid=3828461}} To perform $K$ iterations of ADI, an initial guess $X^\{(0)\}$ is required, as well as $K$ shift parameters, $\backslash \{\; (\backslash alpha\_\{j\},\; \backslash beta\_\{j\})\backslash \}\_\{j\; =\; 1\}^\{K\}$.

If $A\; \backslash in\; \backslash mathbb\{C\}^\{m\; \backslash times\; m\}$ and $B\; \backslash in\; \backslash mathbb\{C\}^\{n\; \backslash times\; n\}$, then $AX\; -\; XB\; =\; C$ can be solved directly in $\backslash mathcal\{O\}(m^3\; +\; n^3)$ using the Bartels-Stewart method.{{Cite book|title=Matrix computations|last1=Golub|first1= G.|publisher=Johns Hopkins University|last2=Van Loan|first2= C|year=1989|isbn=1421407949|edition=Fourth |location=Baltimore|oclc=824733531}} It is therefore only beneficial to use ADI when matrix-vector multiplication and linear solves involving $A$ and $B$ can be applied cheaply.

The equation $AX-XB=C$ has a unique solution if and only if $\backslash sigma(A)\; \backslash cap\; \backslash sigma(B)\; =\; \backslash emptyset$, where $\backslash sigma(M)$ is the spectrum of $M$. However, the ADI method performs especially well when $\backslash sigma(A)$ and $\backslash sigma(B)$ are well-separated, and $A$ and $B$ are normal matrices. These assumptions are met, for example, by the Lyapunov equation $AX\; +\; XA^*\; =\; C$ when $A$ is positive definite. Under these assumptions, near-optimal shift parameters are known for several choices of $A$ and $B$. Additionally, a priori error bounds can be computed, thereby eliminating the need to monitor the residual error in implementation.

The ADI method can still be applied when the above assumptions are not met. The use of suboptimal shift parameters may adversely affect convergence, and convergence is also affected by the non-normality of $A$ or $B$ (sometimes advantageously).{{Cite thesis|last=Sabino|first=J|date=2007|title=Solution of large-scale Lyapunov equations via the block modified Smith method|journal=PhD Diss., Rice Univ.|hdl=1911/20641|type=Thesis}} Krylov subspace methods, such as the Rational Krylov Subspace Method,{{Cite journal|last1=Druskin|first1=V.|last2=Simoncini|first2=V.|date=2011|title=Adaptive rational Krylov subspaces for large-scale dynamical systems|journal=Systems & Control Letters|volume=60|issue=8|pages=546–560|doi=10.1016/j.sysconle.2011.04.013|issn=0167-6911}} are observed to typically converge more rapidly than ADI in this setting, and this has led to the development of hybrid ADI-projection methods.

The problem of finding good shift parameters is nontrivial. This problem can be understood by examining the ADI error equation. After $K$ iterations, the error is given by

$X\; -\; X^\{(K)\}\; =\; \backslash prod\_\{j\; =\; 1\}^K\; \backslash frac\{(A\; -\; \backslash alpha\_j\; I)\}\{(A\; -\; \backslash beta\_j\; I)\}\; \backslash left\; (\; X\; -\; X^\{(0)\}\; \backslash right\; )\; \backslash prod\_\{j\; =\; 1\}^K\; \backslash frac\{(B\; -\; \backslash beta\_j\; I)\}\{(B\; -\; \backslash alpha\_j\; I)\}.$

Choosing $X^\{(0)\}\; =\; 0$ results in the following bound on the relative error:

$\backslash frac\{\backslash left\; \backslash |X\; -\; X^\{(K)\}\; \backslash right\; \backslash |\_2\}\{\backslash |X\backslash |\_2\}\; \backslash leq\; \backslash |\; r\_K(A)\; \backslash |\_2\; \backslash |\; r\_K(B)^\{-1\}\backslash |\_2,\; \backslash quad\; r\_K(M)\; =\; \backslash prod\_\{j\; =\; 1\}^K\; \backslash frac\{(M\; -\; \backslash alpha\_j\; I)\}\{(M\; -\; \backslash beta\_j\; I)\}.$

where $\backslash |\; \backslash cdot\; \backslash |\_2$ is the operator norm. The ideal set of shift parameters $\backslash \{\; (\backslash alpha\_j,\; \backslash beta\_j)\backslash \}\_\{j\; =\; 1\}^K$ defines a rational function $r\_K$ that minimizes the quantity $\backslash |\; r\_K(A)\; \backslash |\_2\; \backslash |\; r\_K(B)^\{-1\}\backslash |\_2$. If $A$ and $B$ are normal matrices and have eigendecompositions $A\; =\; V\_A\backslash Lambda\_AV\_A^*$ and $B\; =\; V\_B\backslash Lambda\_BV\_B^*$, then

$\backslash |\; r\_K(A)\; \backslash |\_2\; \backslash |\; r\_K(B)^\{-1\}\backslash |\_2\; =\; \backslash |\; r\_K(\backslash Lambda\_A)\; \backslash |\_2\; \backslash |\; r\_K(\backslash Lambda\_B)^\{-1\}\backslash |\_2$.

Near-optimal shift parameters are known in certain cases, such as when $\backslash Lambda\_A\; \backslash subset\; [a,\; b]$ and $\backslash Lambda\_B\; \backslash subset\; [c,\; d]$, where $[a,\; b]$ and $[c,\; d]$ are disjoint intervals on the real line. The Lyapunov equation $AX\; +\; XA^*\; =\; C$, for example, satisfies these assumptions when $A$ is positive definite. In this case, the shift parameters can be expressed in closed form using elliptic integrals, and can easily be computed numerically.

More generally, if closed, disjoint sets $E$ and $F$, where $\backslash Lambda\_A\; \backslash subset\; E$ and $\backslash Lambda\_B\; \backslash subset\; F$, are known, the optimal shift parameter selection problem is approximately solved by finding an extremal rational function that attains the value

Z_K(E, F) : = \inf_{r} \frac{ \sup_{z \in E} |r(z)| }{ \inf_{z \in F} |r(z)| },

where the infimum is taken over all rational functions of degree $(K,\; K)$. This approximation problem is related to several results in potential theory,{{Cite book|title=Logarithmic potentials with external fields|last1=Saff|first1=E.B.|last2=Totik|first2=V.|publisher= Springer Science & Business Media |isbn=9783662033296|location=Berlin|oclc=883382758|date = 2013-11-11}}{{Cite journal|last=Gonchar|first=A.A.|date=1969|title=Zolotarev problems connected with rational functions|journal= Mathematics of the USSR-Sbornik|volume=7|issue=4|pages=623–635|bibcode=1969SbMat...7..623G|doi=10.1070/SM1969v007n04ABEH001107}} and was solved by Zolotarev in 1877 for $E$ = [a, b] and $F=-E.${{Cite journal|last=Zolotarev|first=D.I.|date=1877|title=Application of elliptic functions to questions of functions deviating least and most from zero|journal=Zap. Imp. Akad. Nauk. St. Petersburg|volume=30|pages=1–59}} The solution is also known when $E$ and $F$ are disjoint disks in the complex plane.{{Cite journal|last=Starke|first=Gerhard|date=July 1992|title=Near-circularity for the rational Zolotarev problem in the complex plane|journal=Journal of Approximation Theory|volume=70|issue=1|pages=115–130|doi=10.1016/0021-9045(92)90059-w|issn=0021-9045|doi-access=free}}

When less is known about $\backslash sigma(A)$ and $\backslash sigma(B)$, or when $A$ or $B$ are non-normal matrices, it may not be possible to find near-optimal shift parameters. In this setting, a variety of strategies for generating good shift parameters can be used. These include strategies based on asymptotic results in potential theory,{{Cite journal|last=Starke|first=Gerhard|date=June 1993|title=Fejér-Walsh points for rational functions and their use in the ADI iterative method|journal=Journal of Computational and Applied Mathematics|volume=46|issue=1–2|pages=129–141|doi=10.1016/0377-0427(93)90291-i|issn=0377-0427|doi-access=free}} using the Ritz values of the matrices $A$, $A^\{-1\}$, $B$, and $B^\{-1\}$ to formulate a greedy approach,{{Cite journal|last=Penzl|first=Thilo|date=January 1999|title=A Cyclic Low-Rank Smith Method for Large Sparse Lyapunov Equations|journal=SIAM Journal on Scientific Computing|language=en|volume=21|issue=4|pages=1401–1418|doi=10.1137/s1064827598347666|bibcode=1999SJSC...21.1401P |issn=1064-8275}} and cyclic methods, where the same small collection of shift parameters are reused until a convergence tolerance is met. When the same shift parameter is used at every iteration, ADI is equivalent to an algorithm called Smith's method.{{Cite journal|last=Smith|first=R. A.|date=January 1968|title=Matrix Equation XA + BX = C|journal=SIAM Journal on Applied Mathematics|language=en|volume=16|issue=1|pages=198–201|doi=10.1137/0116017|issn=0036-1399}}

In many applications, $A$ and $B$ are very large, sparse matrices, and $C$ can be factored as $C\; =\; C\_1C\_2^*$, where $C\_1\; \backslash in\; \backslash mathbb\{C\}^\{m\; \backslash times\; r\},\; C\_2\; \backslash in\; \backslash mathbb\{C\}^\{n\; \backslash times\; r\}$, with $r\; =\; 1,\; 2$. In such a setting, it may not be feasible to store the potentially dense matrix $X$ explicitly. A variant of ADI, called factored ADI, can be used to compute $ZY^*$, where $X\; \backslash approx\; ZY^*$. The effectiveness of factored ADI depends on whether $X$ is well-approximated by a low rank matrix. This is known to be true under various assumptions about $A$ and $B$.

Historically, the ADI method was developed to solve the 2D diffusion equation on a square domain using finite differences. Unlike ADI for matrix equations, ADI for parabolic equations does not require the selection of shift parameters, since the shift appearing in each iteration is determined by parameters such as the timestep, diffusion coefficient, and grid spacing. The connection to ADI on matrix equations can be observed when one considers the action of the ADI iteration on the system at steady state.

File:ADI-stencil.svg

The traditional method for solving the heat conduction equation numerically is the Crank–Nicolson method. This method results in a very complicated set of equations in multiple dimensions, which are costly to solve. The advantage of the ADI method is that the equations that have to be solved in each step have a simpler structure and can be solved efficiently with the tridiagonal matrix algorithm.

Consider the linear diffusion equation in two dimensions,

: $\{\backslash partial\; u\backslash over\; \backslash partial\; t\}\; =$

\left({\partial^2 u\over \partial x^2 } +

{\partial^2 u\over \partial y^2 }

\right)

= ( u_{xx} + u_{yy} )

The implicit Crank–Nicolson method produces the following finite difference equation:

: $\{u\_\{ij\}^\{n+1\}-u\_\{ij\}^n\backslash over\; \backslash Delta\; t\}\; =$

{1 \over 2(\Delta x)^2}\left(\delta_x^2+\delta_y^2\right)

\left(u_{ij}^{n+1}+u_{ij}^n\right)

where:

: $\backslash Delta\; x\; =\; \backslash Delta\; y$

and $\backslash delta\_p^2$ is the central second difference operator for the p-th coordinate

: $\backslash delta\_p^2\; u\_\{ij\}=u\_\{ij+e\_p\}-2u\_\{ij\}+u\_\{ij-e\_p\}$

with $e\_p=(1,0)$ or $(0,1)$ for $p=x$ or $y$ respectively (and $ij$ a shorthand for lattice points $(i,j)$).

After performing a stability analysis, it can be shown that this method will be stable for any $\backslash Delta\; t$.

A disadvantage of the Crank–Nicolson method is that the matrix in the above equation is banded with a band width that is generally quite large. This makes direct solution of the system of linear equations quite costly (although efficient approximate solutions exist, for example use of the conjugate gradient method preconditioned with incomplete Cholesky factorization).

The idea behind the ADI method is to split the finite difference equations into two, one with the x-derivative taken implicitly and the next with the y-derivative taken implicitly,

: $\{u\_\{ij\}^\{n+1/2\}-u\_\{ij\}^n\backslash over\; \backslash Delta\; t/2\}\; =$

{\left(\delta_x^2 u_{ij}^{n+1/2}+\delta_y^2 u_{ij}^{n}\right)\over \Delta x^2}

: $\{u\_\{ij\}^\{n+1\}-u\_\{ij\}^\{n+1/2\}\backslash over\; \backslash Delta\; t/2\}\; =$

{\left(\delta_x^2 u_{ij}^{n+1/2}+\delta_y^2 u_{ij}^{n+1}\right)\over \Delta y^2}

The system of equations involved is symmetric and tridiagonal (banded with bandwidth 3), and is typically solved using tridiagonal matrix algorithm.

It can be shown that this method is unconditionally stable and second order in time and space.{{Citation | last1=Douglas | first1=J. Jr.| title=On the numerical integration of u_xx+ u_yy= u_t by implicit methods | mr=0071875 | year=1955 | journal=Journal of the Society for Industrial and Applied Mathematics | volume=3 | pages=42–65}}.

There are more refined ADI methods such as the methods of Douglas,{{Citation | last1=Douglas | first1=Jim Jr.| title=Alternating direction methods for three space variables | doi=10.1007/BF01386295 | year=1962 | journal=Numerische Mathematik | issn=0029-599X | volume=4 | issue=1 | pages=41–63| s2cid=121455963 }}. or the f-factor method{{Citation | last1=Chang | first1=M. J. | last2=Chow | first2=L. C. | last3=Chang | first3=W. S. | title=Improved alternating-direction implicit method for solving transient three-dimensional heat diffusion problems | doi=10.1080/10407799108944957 | year=1991 | journal=Numerical Heat Transfer, Part B: Fundamentals | issn=1040-7790 | volume=19 | issue=1 | pages=69–84| bibcode=1991NHTB...19...69C | url=https://zenodo.org/record/1234469 }}. which can be used for three or more dimensions.

The usage of the ADI method as an operator splitting scheme can be generalized. That is, we may consider general evolution equations

: $\backslash dot\; u\; =\; F\_1\; u\; +\; F\_2\; u,$

where $F\_1$ and $F\_2$ are (possibly nonlinear) operators defined on a Banach space.{{cite book|last2=Verwer|first2=Jan|first1=Willem|last1=Hundsdorfer|title=Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations|date=2003|publisher=Springer Berlin Heidelberg|location=Berlin, Heidelberg|isbn=978-3-662-09017-6}}{{cite journal|last1=Lions|first1=P. L.|last2=Mercier|first2=B.|title=Splitting Algorithms for the Sum of Two Nonlinear Operators|journal=SIAM Journal on Numerical Analysis|date=December 1979|volume=16|issue=6|pages=964–979|doi=10.1137/0716071|bibcode=1979SJNA...16..964L}} In the diffusion example above we have $F\_1\; =\; \{\backslash partial^2\; \backslash over\; \backslash partial\; x^2\}$ and $F\_2\; =\; \{\backslash partial^2\; \backslash over\; \backslash partial\; y^2\}$.

It is possible to simplify the conventional ADI method into Fundamental ADI method, which only has the similar operators at the left-hand sides while being operator-free at the right-hand sides. This may be regarded as the fundamental (basic) scheme of ADI method,{{Cite journal|last=Tan|first=E. L.|date=2007|title=Efficient Algorithm for the Unconditionally Stable 3-D ADI-FDTD Method|url=https://www.ntu.edu.sg/home/eeltan/papers/2007%20Efficient%20Algorithm%20for%20the%20Unconditionally%20Stable%203-D%20ADI–FDTD%20Method.pdf|journal=IEEE Microwave and Wireless Components Letters|volume=17|issue=1|pages=7–9|doi=10.1109/LMWC.2006.887239|hdl=10356/138245 |s2cid=29025478}}{{Cite journal|last=Tan|first=E. L.|date=2008|title=Fundamental Schemes for Efficient Unconditionally Stable Implicit Finite-Difference Time-Domain Methods|url=https://www.ntu.edu.sg/home/eeltan/papers/2008%20Fundamental%20Schemes%20for%20Efficient%20Unconditionally%20Stable%20Implicit%20Finite-Difference%20Time-Domain%20Methods.pdf|journal=IEEE Transactions on Antennas and Propagation|volume=56|issue=1|pages=170–177|doi=10.1109/TAP.2007.913089|arxiv=2011.14043|bibcode=2008ITAP...56..170T|hdl=10356/138249 |s2cid=37135325}} with no more operator (to be reduced) at the right-hand sides, unlike most traditional implicit methods that usually consist of operators at both sides of equations. The FADI method leads to simpler, more concise and efficient update equations without degrading the accuracy of conventional ADI method.

Many classical implicit methods by Peaceman-Rachford, Douglas-Gunn, D'Yakonov, Beam-Warming, Crank-Nicolson, etc., may be simplified to fundamental implicit schemes with operator-free right-hand sides. In their fundamental forms, the FADI method of second-order temporal accuracy can be related closely to the fundamental locally one-dimensional (FLOD) method, which can be upgraded to second-order temporal accuracy, such as for three-dimensional Maxwell's equations {{Cite journal|last=Tan|first=E. L.|date=2007|title=Unconditionally Stable LOD-FDTD Method for 3-D Maxwell's Equations|url=https://www.ntu.edu.sg/home/eeltan/papers/2007%20Unconditionally%20Stable%20LOD-FDTD%20Method%20for%203-D%20Maxwell’s%20Equations.pdf|journal=IEEE Microwave and Wireless Components Letters|volume=17|issue=2|pages=85–87|doi=10.1109/LMWC.2006.890166|hdl=10356/138296 |s2cid=22940993}}{{Cite journal|last1=Gan|first1=T. H.|last2=Tan|first2=E. L.|date=2013|title=Unconditionally Stable Fundamental LOD-FDTD Method with Second-Order Temporal Accuracy and Complying Divergence|url=https://www.ntu.edu.sg/home/eeltan/papers/2013%20Unconditionally%20Stable%20Fundamental%20LOD-FDTD%20Method%20With%20Second-Order%20Temporal%20Accuracy%20and%20Complying%20Divergence.pdf|journal=IEEE Transactions on Antennas and Propagation|volume=61|issue=5|pages=2630–2638|doi=10.1109/TAP.2013.2242036|bibcode=2013ITAP...61.2630G|s2cid=7578037}} in computational electromagnetics. For two- and three-dimensional heat conduction and diffusion equations, both FADI and FLOD methods may be implemented in simpler, more efficient and stable manner compared to their conventional methods.{{Cite journal|last1=Tay|first1=W. C.|last2=Tan|first2=E. L.|last3=Heh|first3=D. Y.|date=2014|title=Fundamental Locally One-Dimensional Method for 3-D Thermal Simulation|journal=IEICE Transactions on Electronics|volume=E-97-C|issue=7|pages=636–644|doi=10.1587/transele.E97.C.636|bibcode=2014IEITE..97..636T|hdl=10220/20410 |hdl-access=free}}{{Cite journal|last1=Heh|first1=D. Y.|last2=Tan|first2=E. L.|last3=Tay|first3=W. C.|date=2016|title=Fast Alternating Direction Implicit Method for Efficient Transient Thermal Simulation of Integrated Circuits|journal=International Journal of Numerical Modelling: Electronic Networks, Devices and Fields|volume=29|issue=1|pages=93–108|doi=10.1002/jnm.2049|hdl=10356/137201 |s2cid=61039449 |hdl-access=free}}

• {{cite journal|accessdate=March 28, 2025

|url=https://www.researchgate.net/publication/265623760 |via=ResearchGate.net

|title=50 Years of ADI Methods: Celebrating the Contributions of Jim Douglas, Don Peaceman, and Henry Rachford

|first1=Adam |last1=Usadi |first2=Clint |last2=Dawson

|journal=SIAM News |volume=39|number=2|date=March 2006}} (Provides a review of the ADI methods and variants over the years.)

{{Reflist}}

{{Numerical PDE}}

{{DEFAULTSORT:Alternating Direction Implicit Method}}

Category:Partial differential equations

Category:Numerical differential equations

Category:Iterative methods