A Poisson Solver Based on Iterations on a Sylvester System

Show more

1. Introduction

Poisson’s equation ${\nabla}^{2}u=f$ , an elliptic partial differential equation [1] , was first published in 1813 in the Bulletin de la Société Philomatique by Siméon-Denis Poisson. The equation has since found wide utility in applications such as electrostatics [2] , fluid dynamics [3] , theoretical physics [4] , and engineering [5] . Due to its expansive applicability in the natural sciences, analytic and efficient approximate solution methods have been sought for nearly two centuries. Analytic solutions to Poisson’s equation are unlikely in most scientific applications because the forcing or boundary conditions on the system cannot be explicitly represented by means of elementary functions. For this reason, numerical approximations have been developed, dating back to the Jacobi method in 1845. The linear systems arising from these numerical approximations are solved either directly, using methods like Gaussian elimination, or iteratively. To this day, there are applications solved by direct solvers and others solved by iterative solvers, depending largely on the structure and size of the matrices involved in the computation. The 1950s and 1960s saw an enormous interest in relaxation type methods, prompted by the studies on optimal relaxation and work by Young, Varga, Southwell, Frankel and others. The books by Varga [6] and Young [7] give a comprehensive guide to iterative methods used in the 1960s and 1970s, and have remained the handbooks used by academics and practitioners alike [8] .

The Problem Description

The Poisson equation on a rectangular domain is given by

${\nabla}^{2}u=f\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\Omega =\left\{x,y|0\le x\le a,0\le y\le b\right\}$ (1)

where $u=u\left(x,y\right)$ is to be solved in the 2D domain $\Omega $ , and $f\left(x,y\right)$ is the forcing function. Typical boundary conditions for this equation are either Dirichlet, where the value of $f\left(x,y\right)$ is specified on the boundary, or Neumann, where the value of the normal derivative is specified on the boundary. These are given mathematically as,

$u={g}_{D}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{or}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\partial u/\partial n\equiv \stackrel{^}{n}\cdot \nabla u={g}_{N}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\partial \Omega ,$ (2)

where $\stackrel{^}{n}$ is the outward unit normal along $\partial \Omega $ and ${g}_{D}$ and ${g}_{N}$ are the function values specified by Dirichlet or Neumann boundary conditions. It is also possible to have mixed boundary conditions along the boundary, where some edges have Dirichlet and some have Neumann, so long as the problem is well-posed. Furthermore, edges could be subject to Robin boundary conditions of the form ${c}_{1}u+{c}_{2}\partial u/\partial n={g}_{R}$ . Any numerical scheme designed to solve the Poisson equation should be robust in its ability to incorporate any form of boundary condition into the solver. A detailed discussion of boundary condition implementation is given in the Appendix. Discretizing (1) using central differences with equal grid size $\Delta x=\Delta y=h$ leads to an $M\times N$ rectangular array of unknown U, such that ${U}_{i\mathrm{,}j}\approx u\left({x}_{i}\mathrm{,}{y}_{j}\right)$ (assuming that a and b are both integer multiples of h). This discretization leads to a linear system of the form $AU+UB=F$ , the Sylvester equation, which can be solved either directly or iteratively. The direct method utilizes the Kronecker product approach [9] , given by

$Ku=f\text{\hspace{1em}}\text{\hspace{0.05em}}\text{where}\text{\hspace{0.05em}}\text{\hspace{0.17em}}K=\text{kron}\left(I,A\right)+\text{kron}\left({B}^{\text{T}},I\right)$ (3)

where
$u$ and
$f$ are appropriately ordered
$MN\times 1$ column vectors obtained from the
$M\times N$ arrays U and F, and K is a sparse
$MN\times MN$ matrix. The Kronecker product,
$\text{kron}\left(P\mathrm{,}Q\right)$ , of any two matrices P and Q is a partitioned matrix whose ij^{th} partition contains matrix Q multiplied by component
${p}_{ij}$ of P. Due to the potentially large size of the system given in (3), direct solvers are not the preferred solution approach. Specifically addressed here is an iterative approach to solving the Sylvester equation,

$A{U}^{\mathrm{int}}+{U}^{\mathrm{int}}B={F}^{\mathrm{int}},$ (4)

where ${U}^{\mathrm{int}}$ is the $m\times n$ array of interior unknowns (not including the known boundary values when Dirichlet boundary conditions are given) with $m=M-2$ and $n=N-2$ . Operator matrices A and B are given by the $O\left({h}^{2}\right)$ finite difference approximation to the second derivative,

$\frac{1}{{h}^{2}}\left[\begin{array}{ccccc}-2& 1& & & \\ 1& -2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& -2& 1\\ & & & 1& -2\end{array}\right].$ (5)

For an array of unknowns ${U}^{\mathrm{int}}\in {\mathbb{R}}^{m\times n}$ , the operator matrices are of dimension $A\in {\mathbb{R}}^{m\times m}$ , and $B\in {\mathbb{R}}^{n\times n}$ . The matrix structure of U should be modified such that the first index i corresponds to the x-direction, and the second index j corresponds to the y-direction. With this orientation, multiplying ${U}^{\mathrm{int}}$ by the matrix A on the left approximates the second derivative in the x-direction, and multiplying ${U}^{\mathrm{int}}$ by the matrix B on the right approximates the second derivative in the y-direction.

2. Sylvester Iterative Scheme

Examining (4) it might seem natural to move one term to the right-hand side of the equation to achieve an iterative scheme such as:

$AU+UB=F\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}AU=F-UB$

$\Rightarrow \text{\hspace{1em}}{U}^{k+1}={A}^{-1}F-{A}^{-1}{U}^{k}B,$ (6)

However, this scheme diverges, and an alternative approach is required to iterate on the Sylvester system. An appropriate method is to break up the iterative scheme into two “half-steps’’ as follows

1) First half-step: ${U}^{k}\to {U}^{\mathrm{*}}$

$AU=F-UB\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}\left(A-\alpha I+\alpha I\right)U=F-UB$

$\Rightarrow \text{\hspace{1em}}\left(A-\alpha I\right){U}^{*}=F-{U}^{k}\left(B+\alpha I\right)$ (7)

2) Second half-step: ${U}^{\mathrm{*}}\to {U}^{k+1}$

$UB=F-AU\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}U\left(B-\beta I+\beta I\right)=F-AU$

$\Rightarrow \text{\hspace{1em}}{U}^{k+1}\left(B-\beta I\right)=F-\left(A+\beta I\right){U}^{*}$ (8)

where ${U}^{\mathrm{*}}$ is some intermediate solution between steps k and $k+1$ . Rearranging, this leads to

$\begin{array}{l}\left(I-\frac{1}{\alpha}A\right){U}^{*}={U}^{k}\left(I+\frac{1}{\alpha}B\right)-\frac{1}{\alpha}F,\\ {U}^{k+1}\left(I-\frac{1}{\beta}B\right)=\left(I+\frac{1}{\beta}A\right){U}^{*}-\frac{1}{\beta}F.\end{array}$ (9)

The iterative scheme (9) is similar to the Alternating Direction Implicit (ADI) formulation [10] , where Poisson’s equation is reformulated to have pseudo-time dependency,

$\frac{\text{d}U}{\text{d}t}=AU+UB-F,$ (10)

which achieves the solution to Equation (4) when it reaches steady-state. This method is separated into two half-steps, the first time step going from time $k\to k+1/2$ treating the x-direction implicitly and the y-direction explicitly. The second time step then goes from time $k+1/2\to k+1$ , treating the y-direction implicitly and the x-direction explicitly. The two half-steps are,

$\begin{array}{l}\frac{{U}^{k+1/2}-{U}^{k}}{\Delta t/2}=\frac{1}{{h}^{2}}\left[A{U}^{k+1/2}+{U}^{k}B\right],\\ \frac{{U}^{k+1}-{U}^{k+1/2}}{\Delta t/2}=\frac{1}{{h}^{2}}\left[A{U}^{k+1/2}+{U}^{k+1}B\right],\end{array}$ (11)

which leads to

$\begin{array}{l}\left(I-\frac{\Delta t}{2}A\right){U}^{k+1/2}={U}^{k}\left(I+\frac{\Delta t}{2}B\right)-\frac{\Delta t}{2}F,\\ {U}^{k+1}\left(I-\frac{\Delta t}{2}B\right)=\left(I+\frac{\Delta t}{2}A\right){U}^{k+1/2}-\frac{\Delta t}{2}F.\end{array}$ (12)

This iteration procedure looks nearly identical to our Sylvester iterations given in (9) with $\Delta t/2$ replaced by the unknown parameters 1/α and 1/β. However in our formulation, there is no pseudo-time dependency introduced. Instead, the eigenvalues of our operator matrices A and B are deflated to produce an iterative scheme that optimally converges, and finding the values of the parameters α and β becomes an optimization problem.

Convergence

After the Sylvester Equation (4) is modified into the iterative system (9), the iterative scheme can be written as a single step by substituting the expression for the intermediate solution ${U}^{\mathrm{*}}$ into the second step of the iterative process; this yields the single update equation for ${U}^{k+1}$ given by

$\begin{array}{c}{U}^{k+1}=\left(I+\frac{1}{\beta}A\right){\left(I-\frac{1}{\alpha}A\right)}^{-1}{U}^{k}\left(I+\frac{1}{\alpha}B\right){\left(I-\frac{1}{\beta}B\right)}^{-1}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left[\frac{1}{\alpha}\left(I+\frac{1}{\beta}A\right){\left(I-\frac{1}{\alpha}A\right)}^{-1}-\frac{1}{\beta}I\right]F{\left(I-\frac{1}{\beta}B\right)}^{-1}.\end{array}$ (13)

Assuming that an exact solution
${U}^{\text{exact}}$ exists that exactly satisfies the linear system (4), i.e.
$A{U}^{\text{exact}}+{U}^{\text{exact}}B=F$ , we define the error between the k^{th} iteration and the exact solution as

${E}^{k}\equiv {U}^{k}-{U}^{\text{exact}}.$ (14)

Finding an update equation for the error is done by subtracting the error at the k^{th} step from the error at the
${\left(k+1\right)}^{st}$ step, noting that the expressions involving the forcing F disappear, we arrive at

${E}^{k+1}=P{E}^{k}Q,$ (15)

where the matrices P and Q are given by

$P=\left(I+\frac{1}{\beta}A\right){\left(I-\frac{1}{\alpha}A\right)}^{-1},\text{\hspace{1em}}Q=\left(I+\frac{1}{\alpha}B\right){\left(I-\frac{1}{\beta}B\right)}^{-1}.$ (16)

Denoting the m eigenvalues of $A\in {\mathbb{R}}^{m\times m}$ by ${\lambda}_{k}^{A}$ , and the n eigenvalues of $B\in {\mathbb{R}}^{n\times n}$ by ${\lambda}_{k}^{B}$ , the corresponding deflated eigenvalues of the iteration matrices P and Q are

$\begin{array}{l}{\lambda}_{k}^{P}=\frac{1+\left({\lambda}_{k}^{A}/\beta \right)}{1-\left({\lambda}_{k}^{A}/\alpha \right)}=\frac{\alpha}{\beta}\left(\frac{\beta +{\lambda}_{k}^{A}}{\alpha -{\lambda}_{k}^{A}}\right),\text{\hspace{1em}}k=1,2,\cdots ,m,\\ {\lambda}_{k}^{Q}=\frac{1+\left({\lambda}_{k}^{B}/\alpha \right)}{1-\left({\lambda}_{k}^{B}/\beta \right)}=\frac{\beta}{\alpha}\left(\frac{\alpha +{\lambda}_{k}^{B}}{\beta -{\lambda}_{k}^{B}}\right),\text{\hspace{1em}}k=1,2,\cdots ,n.\end{array}$ (17)

A sufficient condition for convergence of the iterative process is achieved if the spectral radii of both iteration matrices P and Q are less than one,

$\rho \left(P\right)\equiv \mathrm{max}\left|{\lambda}_{k}^{P}\right|<1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{and}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\rho \left(Q\right)\equiv \mathrm{max}\left|{\lambda}_{k}^{Q}\right|<1.$ (18)

The error at each consecutive iteration is decreased by the product of $\rho \left(P\right)$ and $\rho \left(Q\right)$ ,

$\begin{array}{l}{E}^{k+1}=P{E}^{k}Q,\\ \Vert {E}^{k+1}\Vert =\rho \left(P\right)\rho \left(Q\right)\Vert {E}^{k}\Vert ,\\ \Vert {E}^{k+1}\Vert ={\left(\rho \left(P\right)\rho \left(Q\right)\right)}^{k}\Vert {E}^{0}\Vert ,\end{array}$ (19)

where ${E}^{0}$ is the initial error. Often in practical applications, the exact solution is not known, so the error ${E}^{k}$ cannot be computed directly. In this case, the preferred measure in iterative schemes is given by the residual, which measures the difference of the left and right hand sides of the linear system being solved. This will be further discussed in the Results section.

3. Finding Optimal Parameters α and β

Finding α and β is an optimization problem for achieving the fastest convergence rate of the Sylvester iterative scheme (9). Given the operator matrices A and B and their respective eigenvalues ${\lambda}_{k}^{A}$ and ${\lambda}_{k}^{B}$ , it seems feasible to find optimal values of α and β to minimize the spectral radii of the iteration matrices P and Q given in Equations (17). From (15) the error ${E}^{k+1}$ is found by multiplying by P on the left, and Q on the right, thus the convergence is governed by the spectral radii of both P and Q.

Figure 1 shows the eigenvalues of P and Q for arbitrary m and n, given by Equation (17), plotted vs. the eigenvalues of A and B for some parameters α and β. It can be seen that as ${\lambda}^{A}$ or ${\lambda}^{B}$ get large in magnitude, the values of ${\lambda}^{P}$ or ${\lambda}^{Q}$ approach $-\alpha /\beta $ and $-\beta /\alpha $ , respectively. This implies that if $\alpha \ne \beta $ ,

Figure 1. Eigenvalues ${\lambda}^{P}\left({\lambda}^{A};\alpha ,\beta \right)$ and ${\lambda}^{Q}\left({\lambda}^{B}\mathrm{;}\alpha \mathrm{,}\beta \right)$ vs. ${\lambda}^{Q}$ and ${\lambda}^{B}$ for given constant $\alpha $ and $\beta $ . In this figure, $\alpha \mathrm{>}\beta $ , so $\rho \left(P\right)\mathrm{=}\alpha /\beta \mathrm{>1}$ , and the scheme will diverge.

the high frequency eigenvalues of P and Q, in magnitude, will be greater than one, thus convergence condition (18) will not be satisfied. This provides the restriction for convergence that,

$\alpha =\beta .$ (20)

This optimal value of $\alpha =\beta $ will henceforth be called ${\alpha}^{\mathrm{*}}$ . It is important to note that the operator $A\in {\mathbb{R}}^{m\times m}$ or $B\in {\mathbb{R}}^{n\times n}$ with the larger dimension

$\mathcal{l}\equiv max\left(m\mathrm{,}n\right)\mathrm{,}$ (21)

has a larger range of eigenvalues. Figure 1 shows $m>n$ , (i.e. $\mathcal{l}=m$ ) so it can be seen that $\mathrm{min}\left({\lambda}^{A}\right)<\mathrm{min}\left({\lambda}^{B}\right)$ and $\mathrm{max}\left({\lambda}^{A}\right)>\mathrm{max}\left({\lambda}^{B}\right)$ . This property of the eigenvalues is important when calculating an expression for c, which will soon prove to be a highly useful parameter for an adaptive approach to smoothing. Letting $\alpha =\beta ={\alpha}^{*}$ in (17) gives the following expression for the eigenvalues of P and Q:

$\begin{array}{l}{\lambda}_{k}^{P}=\frac{{\alpha}^{*}+{\lambda}_{k}^{A}}{{\alpha}^{*}-{\lambda}_{k}^{A}},\text{\hspace{1em}}k=1,2,\cdots ,m,\\ {\lambda}_{k}^{Q}=\frac{{\alpha}^{*}+{\lambda}_{k}^{B}}{{\alpha}^{*}-{\lambda}_{k}^{B}},\text{\hspace{1em}}k=1,2,\cdots ,n.\end{array}$ (22)

Finding the optimal parameter
${\alpha}^{\mathrm{*}}$ is done by considering the error reduction of Sylvester iterations on an arbitrary initial condition U_{0}. Assume that U_{0} can be decomposed into its constituent error (Fourier) modes, ranging from low frequency (smooth) to high frequency (oscillatory) modes. Given that U_{0} contains error modes of all frequencies, the most conservative method would be to choose
${\alpha}^{\mathrm{*}}$ such that the spectral radii
$\rho \left(P\right)$ and
$\rho \left(Q\right)$ are minimized over the full range of frequencies. This ensures that all modes of error are efficiently relaxed, and convergence is governed by the product of spectral radii.

Referring to the lower curve in Figure 2, the conservative method of

Figure 2. Eigenvalues ${\lambda}^{P}\left({\lambda}^{A}\mathrm{;}\alpha \right)$ and ${\lambda}^{Q}\left({\lambda}^{B}\mathrm{;}\alpha \right)$ illustrating the quantities involved in computing the optimal parameters ${\alpha}^{\mathrm{*}}$ and ${\alpha}_{\text{\hspace{0.05em}}mg}^{\mathrm{*}}$ , and their respective optimal smoothing regions ${R}_{\text{\hspace{0.05em}}smooth\text{\hspace{0.05em}}}^{\mathrm{*}}$ and ${R}_{smooth\text{\hspace{0.05em}}}^{mg}$ . Note that the upper curve illustrates the method of determining ${\alpha}^{\mathrm{*}}$ in the multigrid formulation, while the lower curve is for determining ${\alpha}^{\mathrm{*}}$ in conservative Sylvester iterations.

determining ${\alpha}^{\mathrm{*}}$ would be to set ${L}_{\mathrm{min}}^{*}={L}_{\mathrm{max}}^{*}$ and according to (22),

$-\left(\frac{{\alpha}^{*}+{\lambda}_{\mathrm{min}}^{A}}{{\alpha}^{*}-{\lambda}_{\mathrm{min}}^{A}}\right)=\left(\frac{{\alpha}^{*}+{\lambda}_{\mathrm{max}}^{A}}{{\alpha}^{*}-{\lambda}_{\mathrm{max}}^{A}}\right).$ (23)

Noting that eigenvalues for all dimensions collapse onto the curves shown in Figure 2, this conservative approach “locks in’’ the value of the larger operator’s spectral radius, thus providing an upper bound for convergence. Figure 2 shows $m>n$ , so $\rho \left(P\right)>\rho \left(Q\right)$ , so convergence will be limited by $\rho \left(P\right)$ . Equation (23) can then be solved for ${\alpha}^{\mathrm{*}}$ giving

${\alpha}^{*}=\sqrt{\left|\mathrm{min}\left({\lambda}_{\mathrm{min}}^{A},{\lambda}_{\mathrm{min}}^{B}\right)\right|\times \left|\mathrm{max}\left({\lambda}_{\mathrm{max}}^{A},{\lambda}_{\mathrm{max}}^{B}\right)\right|},$ (24)

where absolute values are introduced as a reminder that ${\lambda}^{A}\mathrm{,}{\lambda}^{B}$ are negative. This value of the parameter ${\alpha}^{\mathrm{*}}$ most uniformly smooths all frequencies for any arbitrary ${U}_{0}$ containing all frequency modes of error. It can be seen that the spectral radii of P and Q shown in Figure 2 occur at either endpoint, and the minimum amplitude occurs near the intersection of the curve with the axis. Varying the parameter ${\alpha}^{\mathrm{*}}$ in (22) controls the intersection point, and thus creates an effective “optimal smoothing region’’, denoted ${R}_{\text{smooth}}^{\ast}$ . Modes of error associated with this optimal smoothing region will be damped fastest, which makes Sylvester iterations highly adaptive in nature. This adaptive nature of Sylvester iterations lends itself nicely to a multigrid formulation.

The Sylvester multigrid formulation is based on the philosophy that most iterative schemes, including Sylvester iterations, relax high frequency modes fastest, leaving low frequency components relatively unchanged [11] . On all grids traversed by a multigrid V-cycle, the high frequency modes are eliminated fastest by finding the optimal parameter value ${\alpha}_{\text{mg}}^{\ast}$ such that ${L}_{\mathrm{min}}^{\text{mg}}={L}_{\text{mid}}^{\text{mg}}$ , as shown in the upper curve of Figure 2. The height ${L}_{\text{mid}}^{\text{mg}}$ is essentially the distance above the axis associated with the approximate “middle’’ eigenvalue in the range of ${\lambda}^{A}$ or ${\lambda}^{B}$ . This equality gives the following optimal parameter value for the Sylvester multigrid method

${\alpha}_{\text{mg}}^{\ast}=\sqrt{\left|min\left({\lambda}_{\mathrm{min}}^{A}\mathrm{,}{\lambda}_{\mathrm{min}}^{B}\right)\right|\times \left|min\left({\lambda}_{\text{mid}}^{A}\mathrm{,}{\lambda}_{\text{mid}}^{B}\right)\right|}\mathrm{,}$ (25)

where, if $m\ne n$ , the minimum (i.e., most negative) middle eigenvalue ${\lambda}_{\text{mid}}\equiv \left({\lambda}_{\mathrm{min}}+{\lambda}_{\mathrm{max}}\right)/2$ is chosen to shrink the optimal smoothing region $R{\text{\hspace{0.05em}}}_{\text{smooth}}^{\text{mg}}$ such that high frequencies are smoothed most effectively. This choice of optimal parameter can be observed in Figure 2 to drastically decrease the magnitude of ${\lambda}^{P}$ and ${\lambda}^{Q}$ associated with high frequencies which significantly enhance relaxation in accordance with the multigrid philosophy.

To find analytical expressions for ${\alpha}^{\mathrm{*}}$ and ${\alpha}_{\text{mg}}^{\ast}$ , it is necessary to have values for ${\lambda}^{A}$ and ${\lambda}^{B}$ . For Dirichlet boundary conditions, analytical expressions for ${\lambda}^{A}$ and ${\lambda}^{B}$ are derived below, but for Neumann boundary conditions numerical approaches are necessary to find ${\lambda}^{A}$ and ${\lambda}^{B}$ . The operator matrices A and B are each of tridiagonal form,

$\left[\begin{array}{cccccc}{d}_{0}& {d}_{1}& & & & \\ {d}_{1}& {d}_{0}& {d}_{1}& & & \\ & \ddots & \ddots & \ddots & & \\ & & {d}_{1}& {d}_{0}& {d}_{1}& \\ & & & {d}_{1}& {d}_{0}& \end{array}\right]\in {\mathbb{R}}^{p\times p}\mathrm{.}$ (26)

Tridiagonal matrices with constant diagonals, such as A and B for Dirichlet boundary conditions, have analytical expressions for their eigenvalues given by

${\lambda}_{k}={d}_{0}+2{d}_{1}\mathrm{cos}\left(\frac{k\text{\pi}}{p+1}\right),\text{\hspace{1em}}k=1,2,\cdots ,p$ (27)

where p is the arbitrary dimension of the matrix [12] . Neumann boundary conditions alter the upper and lower diagonals of A or B, thus there is no analytical form of eigenvalues for Neumann boundary conditions. Using (5) and (27) gives the following analytic form of the eigenvalues of the tridiagonal matrices A and B,

${\lambda}_{k}=\frac{2}{{h}^{2}}\left(-1+\mathrm{cos}\left(\frac{k\text{\pi}}{p+1}\right)\right),\text{\hspace{1em}}k=1,2,\cdots ,p,$ (28)

which achieves minimum and maximum values given by

${\lambda}_{\mathrm{min}}=\frac{2}{{h}^{2}}\left(-1+\mathrm{cos}\left(\frac{p\text{\pi}}{p+1}\right)\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{and}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\lambda}_{\mathrm{max}}=\frac{2}{{h}^{2}}\left(-1+\mathrm{cos}\left(\frac{\text{\pi}}{p+1}\right)\right),$ (29)

respectively. Using (24), (25), and (29) the analytic expressions for optimal parameters for both conservative and multigrid approaches are given by

$\begin{array}{l}{\alpha}^{*}=\frac{2}{{h}^{2}}\sqrt{\left(-1+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)\left(-1+\mathrm{cos}\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)\right)},\\ {\alpha}_{\text{mg}}^{\ast}=\frac{\sqrt{2}}{{h}^{2}}\sqrt{\left(-1+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)\left(-2+\mathrm{cos}\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)},\end{array}$ (30)

where again $\mathcal{l}\equiv max\left(m\mathrm{,}n\right)$ . Having expressions for ${\alpha}^{\mathrm{*}}$ and ${\alpha}_{\text{mg}}^{\ast}$ allows ${\lambda}^{P}$ and ${\lambda}^{Q}$ to be found analytically using (22) which subsequently allows the spectral radii of the iteration matrices P and Q to be calculated. Knowing the spectral radii of the iteration matrices P and Q is highly advantageous, as it allows for an analysis of the Sylvester iterative scheme.

4. Analysis

The analysis of standard Sylvester iterations can be performed and describes the error reduction with each consecutive iteration using (19). Having the optimal parameters given by (30) and eigenvalues of P and Q in (22), the spectral radii can be calculated to be

$\begin{array}{l}\rho \left(P\right)=\frac{-1+\mathrm{cos}\left(\frac{m\text{\pi}}{m+1}\right)+\sqrt{\left(-1+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)\left(-1+\mathrm{cos}\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)\right)}}{-1+\mathrm{cos}\left(\frac{m\text{\pi}}{m+1}\right)-\sqrt{\left(-1+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)\left(-1+\mathrm{cos}\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)\right)}},\\ \rho \left(Q\right)=\frac{-1+\mathrm{cos}\left(\frac{n\text{\pi}}{n+1}\right)+\sqrt{\left(-1+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)\left(-1+\mathrm{cos}\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)\right)}}{-1+\mathrm{cos}\left(\frac{n\text{\pi}}{n+1}\right)-\sqrt{\left(-1+\mathrm{cos}\left(\frac{\mathcal{l}\text{\pi}}{\mathcal{l}+1}\right)\right)\left(-1+\mathrm{cos}\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)\right)}}.\end{array}$ (31)

Rewriting the last expression of (19), we see that

$\frac{\Vert {E}^{k}\Vert}{\Vert {E}^{0}\Vert}~{\left(\rho \left(P\right)\rho \left(Q\right)\right)}^{k}\mathrm{.}$ (32)

If we want to reduce our error to $\Vert {E}^{k}\Vert ~\u03f5\Vert {E}^{0}\Vert $ and we wish to know how many iterations it will take to achieve this error reduction, using (32) we set ${\left(\rho \left(P\right)\rho \left(Q\right)\right)}^{k}~\u03f5$ , and solving for k, we find it will take

$k~\frac{log\left(\u03f5\right)}{log\left(\rho \left(P\right)\rho \left(Q\right)\right)}$ (33)

iterations to reduce the error by $\u03f5$ . Here log can be with respect to any base, as long as the same one is used in both the numerator and denominator; e.g., the natural log can be used. Recall that the exact solution ${U}^{\text{exact}}$ of (4) is only an approximate solution of the differential Equation (1) we are actually solving. Due to this, we can only expect accuracy of the truncation error of the approximation. With an $O\left({h}^{2}\right)$ method, ${U}_{i\mathrm{,}j}^{\text{exact}}$ differs from $U\left({x}_{i}\mathrm{,}{y}_{j}\right)$ on the order of ${h}^{2}$ so we cannot achieve better accuracy than this no matter how well we solve the linear system. Thus, it is practical to take $\u03f5$ to be something proportional to the expected global error, e.g. $\u03f5=C{h}^{2}$ for some fixed C [12] .

To calculate the order of work required asymptotically as $h\to 0$ , (i.e. $m\to \infty $ ) using (33) and our choice for $\u03f5$ , we see that

$k~\frac{log\left(C\right)+2log\left(h\right)}{log\left(\rho \left(P\right)\rho \left(Q\right)\right)}\mathrm{.}$ (34)

The expressions for $\rho \left(P\right)$ and $\rho \left(Q\right)$ in (31) contain several cosine terms which can be Taylor expanded about different values. Cosines with arguments like $\text{\pi}x$ can be expanded about $x=1$ or $x=0$ depending on the form of x, namely

$\begin{array}{l}cos\left(\text{\pi}x\right)~-1+\frac{{\text{\pi}}^{2}}{2}{\left(x-1\right)}^{2}+\mathcal{O}\left({\left(x-1\right)}^{3}\right)\text{\hspace{1em}}\text{\hspace{0.05em}}\text{for}\text{\hspace{0.05em}}\text{\hspace{0.17em}}x\approx \mathrm{1,}\\ cos\left(\text{\pi}x\right)~1-\frac{{\text{\pi}}^{2}}{2}{x}^{2}+\mathcal{O}\left({x}^{3}\right)\text{\hspace{1em}}\text{\hspace{0.05em}}\text{for}\text{\hspace{0.05em}}\text{\hspace{0.17em}}x\approx \mathrm{0,}\end{array}$ (35)

where, from (31), the form of x is something like $m/\left(m+1\right)$ or $1/\left(m+1\right)$ , which clearly approach one or zero, respectively, in the limit that $m\to \infty $ . Using these expansions, along with the fact that $1/\left(1-x\right)~1+x+O\left({x}^{2}\right)$ for $x\ll 1$ to simplify the spectral radii, we arrive at the following

$\rho \left(P\right)~\rho \left(Q\right)~1-\frac{\text{\pi}}{\mathcal{l}+1}+\frac{1}{4}{\left(\frac{\text{\pi}}{\mathcal{l}+1}\right)}^{2}\mathrm{,}$ (36)

when $m\mathrm{,}n\gg 1$ . Since $h=1/\left(m+1\right)$ , (34) combined with (36) gives the following order of work needed for convergence to within $\u03f5~C{h}^{2}$ :

$k~\frac{-2log\left(m+1\right)}{2log\left(1-\frac{\pi}{\mathcal{l}+1}\right)}~\frac{\mathcal{l}}{\text{\pi}}log\left(m\right)\mathrm{,}$ (37)

where only linear terms are used from (36), and the latter simplified expression can be deduced by using the property that $log\left(1+x\right)~x+O\left({x}^{2}\right)$ for $x\ll 1$ . Note that when $m=n$ , the order of work for Sylvester iterations is $k~\left(m/\text{\pi}\right)log\left(m\right)$ , which is comparable to the work necessary for the Successive Over-Relaxation (SOR) algorithm to solve Poisson’s equation [12] . This will be our basis for comparison in the Results section for standard Sylvester iterations.

5. Results

Problems solved by Sylvester iterations can, in general, be written shorthand as $\mathcal{L}U=F$ , where $\mathcal{L}$ is a linear operator. In the case of Poisson’s equation, $\mathcal{L}$ is the Laplacian operator. As an error measure, the discrete ${\Vert \text{\hspace{0.05em}}\cdot \text{\hspace{0.05em}}\Vert}_{2}$ norm of the residual, $r\equiv F-\mathcal{L}U$ , can be measured at each iteration. This number provides the stopping criterion for our iterative schemes, namely the iterations are run until

$\Vert {r}^{\left(k\right)}\Vert =\Vert F-\mathcal{L}{U}^{\left(k\right)}\Vert <\text{tol}\times \Vert {r}^{\left(0\right)}\Vert ,$ (38)

where ${r}^{\left(0\right)}$ is the initial residual, and tol is the tolerance. The tolerance is set to machine precision $\text{tol}~{10}^{-16}$ to illustrate the asymptotic convergence rate,

${q}^{\left(k\right)}=\frac{\Vert {r}^{\left(k\right)}\Vert}{\Vert {r}^{\left(k-1\right)}\Vert},$ (39)

however, in practice, the discretization error $O\left({h}^{2}\right)$ is the best accuracy that can be expected. These numerical results were run using MATLAB on a 1.5 GHz Mac PowerPC G4. The model problem that is solved is given by

$\begin{array}{l}{\nabla}^{2}u=-2\left[{y}^{2}\left(1-6{x}^{2}\right)\left(1-{y}^{2}\right)+{x}^{2}\left(1-6{y}^{2}\right)\left(1-{x}^{2}\right)\right]\text{\hspace{1em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\Omega ,\\ u=0\text{\hspace{1em}}\text{\hspace{0.05em}}\text{on}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\partial \Omega ,\end{array}$ (40)

where $\Omega =\left\{x,y|0\le x\le 1,0\le y\le 1\right\}$ , and whose exact solution

${u}^{\text{exact}}\left(x\mathrm{,}y\right)=\left({x}^{2}-{x}^{4}\right)\left({y}^{4}-{y}^{2}\right)\mathrm{,}$ (41)

is known so errors can be computed [11] . This model problem is used to show performance of both standard and multigrid Sylvester iterations. In all cases, the initial guess ${U}^{\left(0\right)}$ of the iterative scheme is a normalized random array and can be assumed to contain all modes of error.

5.1. Standard Sylvester Iterations

For comparison, standard Sylvester iterations were tested against Successive Over-Relaxation (SOR) with Chebyshev acceleration (see e.g., [13] ). In SOR with Chebyshev acceleration, one uses odd-even ordering of the grid and changes the relaxation parameter $\omega $ at each half-step, which converges to the optimal relaxation parameter. The results are shown in Table 1. It is important to note that SOR iterations involve no matrix inversions, whereas Sylvester iterations do, thus the CPU time measure might not be an appropriate gauge for this particular comparison. From these results, it is clear that standard Sylvester iterations are comparable to SOR, and in most cases, converge to within tolerance in fewer iterations than SOR. An unfortunate artifact of standard iterative schemes is that as system size increases, so does the spectral radii governing the convergence. This can be observed in Table 1, as asymptotic convergence rates for each method ${q}_{\text{sylvester}}$ and ${q}_{\text{sor}}$ steadily increase, thus requiring higher numbers of iterations to solve within tolerance. The number of Sylvester iterations required to converge is consistent with the predicted number of iterations given by Equation (33) letting $\u03f5=\text{tol}={10}^{-16}$ .

5.2. Multigrid Sylvester Iterations

In multigrid Sylvester iterations, the performance of the $V\left({\nu}_{1}\mathrm{,}{\nu}_{2}\right)$ -cycle using Sylvester iterations is compared to that using the traditional Gauss-Seidel (GS)

Table 1. Sylvester iterations vs. successive over-relaxation (SOR).

Table 2. Multigrid sylvester iterations vs. multigrid gauss-seidel (GS) iterations.

iterations. The parameter ${\nu}_{1}$ represents the number of smoothing iterations done on each level of the downward branch of the V-cycle, while ${\nu}_{2}$ represents the number done on the upward branch. In practice, common choices are $\nu ={\nu}_{1}+{\nu}_{2}\le 3$ , so our performance is based on the $V\left(\mathrm{2,1}\right)$ -cycle [14] . In each case, the V-cycle descends to the coarsest grid having gridwidth ${h}_{0}=1/2$ , and in the Sylvester implementation, the value of ${\alpha}_{\text{mg}}^{\ast}$ is calculated to smooth high frequencies most effectively on each grid traversed by the cycle. The results are shown in Table 2. It can be seen that the asymptotic convergence rates ${q}_{\text{sylvester}}^{\text{mg}}$ and ${q}_{\text{gs}}^{\text{mg}}$ reach steady values independent of the gridwidth h. This is characteristic of multigrid methods, and enables the optimality of the multigrid method. It is clear when comparing the CPU times of the Sylvester multigrid formulation in Table 2 with standard Sylvester iterations in Table 1 that the multigrid framework is substantially faster (e.g., 30 times faster than standard iterations for a grid of size $256\times 256$ ). It can also be seen that the asymptotic convergence rates are such that ${q}_{\text{sylvester}}^{\text{mg}}<{q}_{\text{gs}}^{\text{mg}}$ , thus convergence is met in fewer $V\left(\mathrm{2,1}\right)$ cycles using Sylvester smoothing versus Gauss-Seidel smoothing.

6. Conclusion

Sylvester iterations provide an alternative iterative scheme to solve Poisson’s equation that is comparable to SOR in the number of iterations necessary to converge, namely converging to discretization accuracy within $k~\left(m/\text{\pi}\right)log\left(m\right)$ iterations. The true benefit of the Sylvester iterations, however, comes from its adaptive ability to smooth any range of error frequencies, thus being a perfect candidate for smoothing in a multigrid framework. Multigrid $V\left(\mathrm{2,1}\right)$ -cycles using Sylvester smoothing have an asymptotic convergence rate of ${q}_{\text{sylvester}}^{\text{mg}}=0.069$ (versus ${q}_{\text{gs}}^{\text{mg}}=0.083$ for Gauss-Seidel smoothing) and indicate significant improvement in efficiency over standard Sylvester iterations.

Appendix

1) Boundary condition implementation

Solving the Poisson equation using Sylvester iterations lends itself nicely to boundary condition implementation. Dirichlet boundary conditions of the form

$u\left(0,y\right)={u}_{1}\left(y\right),\text{\hspace{1em}}u\left(a,y\right)={u}_{2}\left(y\right),\text{\hspace{1em}}u\left(x,0\right)={u}_{3}\left(x\right),\text{\hspace{1em}}u\left(x,b\right)={u}_{4}\left(x\right),$ (1)

where ${u}_{1}\mathrm{,}{u}_{2}\mathrm{,}{u}_{3}$ and ${u}_{4}$ are functions describing the edges of U, can be implemented as follows. The unknown values in the Sylvester system given in Equation (4) are the $m\times n$ array of interior values, where $m=M-2,n=N-2$ . It is possible to incorporate Dirichlet boundary conditions directly into this interior system by examining the partitioned matrix product, for example AU, given by

(2)

with UB taking an analogous partitioned form. Multiplying through by ${h}^{2}$ associated with the operator matrices A and B, the partitioned Sylvester system for internal unknowns gives

${A}_{L}\cdot {u}_{1}^{\text{T}}+\left({A}^{\mathrm{int}}\right)\left({U}^{\mathrm{int}}\right)+{A}_{R}\cdot {u}_{2}^{\text{T}}+{u}_{3}\cdot {B}_{T}^{\text{T}}+\left({U}^{\mathrm{int}}\right)\left({B}^{\mathrm{int}}\right)+{u}_{4}\cdot {B}_{B}^{\text{T}}=\left({h}^{2}\right)\left({F}^{\mathrm{int}}\right)\mathrm{,}$ (3)

where all matrix-vector products are $m\times n$ outer products. Note that the product AU incorporates Dirichlet boundary conditions in the x-direction, and UB incorporates Dirichlet boundary conditions in the y-direction. Combining the partitioned systems incorporating both A and B matrix multiplications and boundary conditions yields

$\left({A}^{\mathrm{int}}\right)\left({U}^{\mathrm{int}}\right)=\left({h}^{2}\right)\left({F}^{\mathrm{int}}\right)-\left({A}_{L}\cdot {u}_{1}^{\text{T}}+{A}_{R}\cdot {u}_{2}^{\text{T}}\right)-\left({u}_{3}\cdot {B}_{T}^{\text{T}}+{u}_{4}\cdot {B}_{B}^{\text{T}}\right)\mathrm{,}$ (4)

which is an $m\times n$ linear system for ${U}^{\mathrm{int}}$ .

For Neumann boundary conditions, the edge at which the condition is imposed becomes part of the internal unknowns in the Sylvester system. As an example, consider a Neumann boundary condition given by

$\frac{\partial u}{\partial x}=g\left(y\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{on}\text{\hspace{0.05em}}\text{\hspace{0.17em}}x=0.$ (5)

Staying within the finite difference formulation of derivatives and letting $g\left({y}_{j}\right)\equiv {g}_{j}$ , this condition can be discretized and approximated with the $O\left({h}^{2}\right)$ central difference approximation, which yields

${\left(\frac{\partial u}{\partial x}\right)}_{i,j}\approx \frac{{U}_{i+1,j}-{U}_{i-1,j}}{2h}={g}_{j}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.05em}}\text{\hspace{0.17em}}i=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\le j\le N.$ (6)

For a Neumann condition along the edge $x=0$ , the row vector ${u}_{1}^{\text{T}}$ described in (2) becomes a part of the internal array of unknowns ${U}^{\mathrm{int}}$ . In order to implement this finite difference on the edge, we need to introduce a ghost layer with index $i=-1$ , and pair Equation (6) to the second derivative operator AU for $i=0$ . This gives

$\begin{array}{l}\frac{{U}_{1,j}-{U}_{-1,j}}{2h}={g}_{j}\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}{U}_{-1,j}={U}_{1,j}-2h{g}_{j},\\ {\left(\frac{{\partial}^{2}u}{\partial {x}^{2}}\right)}_{0,j}\approx \frac{\left({U}_{1,j}-2h{g}_{j}\right)-2{U}_{0,j}+{U}_{1,j}}{{h}^{2}}=\frac{2{U}_{1,j}-2{U}_{0,j}}{{h}^{2}}-2\frac{{g}_{j}}{h},\end{array}$ (7)

which leads to the following partitioned form of AU,

(8)

where the additional term $-2{g}_{j}/h$ is taken to the right hand side of the Sylvester system such that ${F}_{0,j}^{\mathrm{int}}={F}_{0,j}^{\mathrm{int}}+2{g}_{j}/h$ for $0<j<N$ . Comparing (8) to (2), the size of ${U}^{\mathrm{int}}$ changes from $m\times n$ to $\left(m+1\right)\times n$ , and ${A}_{\mathrm{0,1}}$ is changed from 1 to 2 (shown boldface in (8)), and the right hand side is slightly modified along that edge. Similarly, any edge with a Neumann condition can be handled in this fashion. It is clear that both Dirichlet and Neumann boundary conditions are very simple to implement in the Sylvester iteration method, and only slightly modify the structure of the arrays involved.

References

[1] Douglass, C., Hasse, G. and Langer, U. (2003) A Tutorial on Elliptic PDE Solvers and Their Parallelization. Society for Industrial and Applied Math, Philadelphia.

https://doi.org/10.1137/1.9780898718171

[2] Feig, M., Onufriev, A., Lee, M.S., Im, W., Case, D.A. and Brooks III, C.L. (2003) Performance Comparison of Generalized Born and Poisson Methods in the Calculation of Electrostatic Solvation Energies for Protein Structures. Journal of Computational Chemistry, 25, 265-284.

https://doi.org/10.1002/jcc.10378

[3] Ravoux, J.F., Nadim, A. and Haj-Hariri, H. (2003) An Embedding Method for Bluff Body Flows: Interactions of Two Side-by-Side Cylinder Wakes. Theoretical and Computational Fluid Dynamics, 16, 433-466.

https://doi.org/10.1007/s00162-003-0090-4

[4] Trellakis, A., Galick, A.T., Pacelli, A. and Ravaioli, U. (1997) Iteration Scheme for the Solution of the Two-Dimensional Schrodinger-Poisson Equations in Quantum Structures. Journal of Applied Physics, 81, 7880-7884.

https://doi.org/10.1063/1.365396

[5] Saraniti, M., Rein, A., Zandler, G., Vogl, P. and Lugli, P. (1996) An Efficient Multigrid Poisson Solver for Device Simulations. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 15, 141-150.

https://doi.org/10.1109/43.486661

[6] Varga, R.S. (1962) Matrix Iterative Analysis. Prentice-Hall, New Jersey.

[7] Young, D.M. (1971) Iterative Solution of Large Linear Systems. Academic Press, New York.

[8] Brezinski, C. and Wuytack, L. (2001) Numerical Analysis: Historical Developments in the 20th Century. Elsevier Science B.V., Netherlands.

[9] Van Loan, C.F. (2000) The Ubiquitous Kronecker Product. Journal of Computational and Applied Mathematics, 123, 85-100.

https://doi.org/10.1016/S0377-0427(00)00393-9

[10] Peaceman, D.W. and Rachford, H.H. (1955) The Numerical Solution of Parabolic and Elliptic Differential Equations. Journal of the Society for Industrial and Applied Mathematics, 3, 28-41.

https://doi.org/10.1137/0103003

[11] Briggs, W.L., Henson, V.E. and McCormick, S.F. (2000) A Multigrid Tutorial. 2nd Edition, Society for Industrial and Applied Math, Philadelphia.

https://doi.org/10.1137/1.9780898719505

[12] LeVeque, R.J. (2007) Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. Society for Industrial and Applied Math, Philadelphia.

https://doi.org/10.1137/1.9780898717839

[13] Press, W., Vetterling, W., Teukolsky, S. and Flannery, B. (1992) Numerical Recipes in Fortran. 2nd Edition, Cambridge University Press, New York.

[14] Trottenberg, U., Oosterlee, C. and Schuller, A. (2001) Multigrid. Elsevier Academic Press, London.