A Poisson Solver Based on Iterations on a Sylvester System

We present an iterative scheme for solving Poisson’s equation in 2D. Using finite differences, we discretize the equation into a Sylvester system, AU UB F + = , involving tridiagonal matrices A and B. The iterations occur on this Sylvester system directly after introducing a deflation-type parameter that enables optimized convergence. Analytical bounds are obtained on the spectral radii of the iteration matrices. Our method is comparable to Successive Over-Relaxation (SOR) and amenable to compact programming via vector/array operations. It can also be implemented within a multigrid framework with considerable improvement in performance as shown herein.


Introduction
Poisson's equation 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 where ( ) , u u x y = is to be solved in the 2D domain Ω , and ( ) , f x y is the forcing function.Typical boundary conditions for this equation are either Dirichlet, where the value of ( ) , f x y is specified on the boundary, or Neumann, where the value of the normal derivative is specified on the boundary.These are given mathematically as, ôr on , where n is the outward unit normal along ∂Ω and D (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 ite- ratively.The direct method utilizes the Kronecker product approach [9], given by ( ) where u and f are appropriately ordered  .
For an array of unknowns int 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 int U by the matrix A on the left approximates the second derivative in the x-direction, and multiplying int U by the matrix B on the right approximates the second derivative in the y-direction.

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: 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: 2) Second half-step: where * U is some intermediate solution between steps k and 1 k + .Rearrang- ing, this leads to 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, 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 treating the x-direction implicitly and the y-direction explicitly.
The second time step then goes from time 1 2 1 k k + → + , treating the y-direction implicitly and the x-direction explicitly.The two half-steps are, which leads to This iteration procedure looks nearly identical to our Sylvester iterations given in (9) with 2 t ∆ 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 into the second step of the iterative process; this yields the single update equation for Assuming that an exact solution exact U exists that exactly satisfies the linear system (4), i.e. k k Finding an update equation for the error is done by subtracting the error at the k th step from the error at the ( ) where the matrices P and Q are given by Denoting the m eigenvalues of , . 1 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, ( ) ( ) max 1 and max 1.
The error at each consecutive iteration is decreased by the product of ( ) where 0 E is the initial error.Often in practical applications, the exact solution is not known, so the error k E 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.

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 A k λ and B k λ , 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 1 k E + is found by multip- lying 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 A λ or B λ get large in magnitude, the values of P λ or Q λ approach α β − and β α − , respectively.This implies that if α β ≠ , Applied Mathematics , 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, .
This optimal value of α β has a larger range of eigenvalues.Figure 1 shows m n > , (i.e.m =  ) so it can be seen that ( ) ( ) . 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 * α β α = = in (17) gives the following expression for the eigenvalues of P and Q: , .
Finding the optimal parameter * α is done by considering the error reduc- tion 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 * α such that the spectral radii ( ) 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 .
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 ( ) ( ) , so convergence will be limited by ( ) where absolute values are introduced as a reminder that , This value of the parameter * α most uniformly smooths all frequencies for any arbitrary 0 U 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 * α in (22) controls the intersection point, and thus creates an effective "optimal smoothing region'', denoted smooth R * .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 mg α * such that λ and B λ .The operator matrices A and B are each of tridiagonal form, Tridiagonal matrices with constant diagonals, such as A and B for Dirichlet boundary conditions, have analytical expressions for their eigenvalues given by 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, which achieves minimum and maximum values given by respectively.Using (24), (25), and (29) the analytic expressions for optimal parameters for both conservative and multigrid approaches are given by λ 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.

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 ( ) Rewriting the last expression of (19), we see that If we want to reduce our error to and we wish to know how many iterations it will take to achieve this error reduction, using (32) we set , and solving for k, we find it will take ( ) ( ) ( ) ( ) iterations to reduce the error by  .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 exact U 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 ( ) U x y on the order of 2 h so we cannot achieve better accuracy than this no matter how well we solve the linear system.Thus, it is practical to take  to be something proportional to the expected global error, e.g.
for some fixed C [12].
To calculate the order of work required asymptotically as 0 h → , (i.e.m → ∞ ) using (33) and our choice for  , we see that The expressions for ( ) where, from (31), the form of x is something like ( ) ) which clearly approach one or zero, respectively, in the limit that m → ∞ .Using these expansions, along with the fact that ( ) ( ) simplify the spectral radii, we arrive at the following when , 1 m n  .Since ( ) + , (34) combined with (36) gives the following order of work needed for convergence to within where only linear terms are used from (36), and the latter simplified expression can be deduced by using the property that ( ) ( ) Note that when m n = , the order of work for Sylvester iterations is m , 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.

Results
Problems solved by Sylvester iterations can, in general, be written shorthand as , where  is a linear operator.In the case of Poisson's equation,  is the Laplacian operator.As an error measure, the discrete 2 ⋅ norm of the residual, r F U ≡ − , can be measured at each iteration.This number provides the stopping criterion for our iterative schemes, namely the iterations are run until where ( ) 0 r is the initial residual, and tol is the tolerance.The tolerance is set to machine precision 16 tol ~10 − to illustrate the asymptotic convergence rate, however, in practice, the discretization error ( ) O h is the best accuracy that can be expected.These numerical results were run using MATLAB on a 1. 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 ( ) 0 U of the iterative scheme is a normalized random array and can be assumed to all modes of error.

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 ω at each half-step, which converges to the optimal re- laxation 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 sylvester q and sor q steadily increase, thus requiring higher numbers of ite- rations 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 16 tol 10 − = =  .

Multigrid Sylvester Iterations
In multigrid Sylvester iterations, the performance of the ( ) , V ν ν -cycle using Sylvester iterations is compared to that using the traditional Gauss-Seidel (GS)  3 ν ν ν = + ≤ , so our performance is based on the ( ) -cycle [14].In each case, the V-cycle descends to the coarsest grid having gridwidth 0 1 2 h = , and in the Sylvester implementation, the value of mg α * 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 mg sylvester q and mg gs q 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 256 × ).It can also be seen that the asymptotic convergence rates are such that mg mg sylvester gs q q < , thus convergence is met in fewer ( ) 2,1 V cycles using Sylvester smoothing versus Gauss-Seidel smoothing.

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 g and N g 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 1 2 R c u c u n g + ∂ ∂ = .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 x y h ∆ =∆ = leads to an M N × define the error between the k th iteration and the exact solution as exact .
corresponding deflated eigenvalues of the iteration matrices P and Q are

Figure 2
Figure 2. Eigenvalues contain several cosine terms which can be Taylor expanded about different values.Cosines with arguments like πx can be expanded about 1 x = or 0 x = depending on the form of x, namely 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 ( ) Seidel smoothing) and indicate significant improvement in efficiency over standard Sylvester iterations.Applied Mathematics implement this finite difference on the edge, we need to introduce a ghost layer with index 1 i = − , and pair Equation (6) to the second derivative operator AU for the following partitioned form of AU,(8) where the additional term 2 j g h − is taken to the right hand side of the Sylvester system such that int int N < < .Comparing (8) to (2), the size of int U changes from m n × to ( )1 m n + × ,and 0,1 A 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.

−
, of any two matrices P and Q is a partitioned matrix whose ij th partition contains matrix Q multiplied by component ij p 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, .Operator matrices A and B are given by the ( ) 23) can then be solved for * 5GHz Applied Mathematics