Yue Hao

Valeria Simoncini
email: valeria.simoncini@unibo.it

Valeria Simoncini Matrix Equation

MATRIX EQUATION SOLVING OF PDES ON REGULAR DOMAINS

Keywords: Elliptic partial differential equations, Schwarz-Christoffel mapping, Matrix equations, Iterative methods AMS subject classifications, 65F10, 65N22

whether they are published or not. The documents may come

1. Introduction. We are interested in exploring unconventional algebraic strategies for numerically solving two-dimensional partial differential equations (PDEs) of the type (1.1) -α 1 u xx -α 2 u yy + w • ∇u + βu = f, (x, y) ∈ Ω, where u = u(x, y), f = f (x, y), Ω ⊂ R 2 is the interior of a polygon, the vector w = (w 1 (x, y), w 2 (x, y)) accounts for the convection, and α 1 , α 2 , β are coefficients possibly also depending on the space variables. Dirichlet conditions are considered throughout the whole domain boundaries.

Commonly employed discretization methods such as finite elements lead to largescale sparse linear systems to be solved, typically by means of preconditioned Krylov subspace methods; see, e.g., [START_REF] Elman | Finite Elements and Fast Iterative Solvers, with applications in incompressible fluid dynamics[END_REF], [START_REF] Olshanskii | Iterative methods for linear systems, Theory and Applications[END_REF], [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Recently, the use of tensor basis discretizations such as finite differences (FD), spectral methods and Isogeometric Analysis (ISA) (see, e.g., [START_REF] Canuto | Contraction and optimality properties of an adaptive Legendre-Galerkin method: the multi-dimensional case[END_REF], [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF], [START_REF] Sangalli | Isogeometric preconditioners based on fast solvers for the Sylvester equation[END_REF]) has shown that under certain hypotheses on the domain and coefficients, it is possible to write the discretized version of (1.1) in terms of linear matrix equations; in addition to a better understanding of the structural properties of the discretized quantities, this perspective can allow significant savings not only in memory requirements, but also in computational costs. In particular, the FD approach presented in [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF] focussed on the case of rectangular domains.

Our aim is to significantly expand the matrix-oriented FD discretization approach by addressing the case where Ω is a more general domain, such as a polygon. In our context, this can be addressed by using structured grids. A large range of computational strategies to associate the physical domain with a simpler reference domain such as rectangles or disks have been proposed in the literature. Classically, bilinear projectors have been used in case of curved boundaries and are at the basis of transfinite interpolation in two and three dimensions; see, e.g., [START_REF] Knupp | The fundamentals of grid generation[END_REF], [START_REF] Farrashkhalvat | Basic Structured Grid Generation[END_REF]. Nonetheless it is recognized that, at least from a mathematical standpoint, it may be advantageous to embed the problem in the complex plane and use conformal and quasi-conformal mappings, although their used is strictly constrained to two-dimensional problems.

Other explored approaches include variational and other PDE based methods, especially elliptic grid generators, extensively employed in fluid and air dynamics in the 1980s. The literature is extremely vast, here we give the reader a couple of pointers [START_REF] Castillo | Mathematical Aspects of NUMERICAL grid generation[END_REF], [START_REF] Knupp | The fundamentals of grid generation[END_REF], [START_REF] Thompson | Numerical grid generation[END_REF]. In this paper we probe the applicability of conformal Schwarz-Christoffel (SC) mappings, which have several theoretical and computational advantages, not least the fact that quite reliable easy-to-use Matlab software is available [START_REF] Driscoll | Algorithm 756: A MATLAB Toolbox for Schwarz-Christoffel mapping[END_REF], [START_REF] Driscoll | Algorithm 843: Improvements to the Schwarz-Christoffel Toolbox for MAT-LAB[END_REF]; other available software leading to the same matrix setting could effectively be employed.

Our aim is to provide a proof of concept methodology that leads to new algebraic formulations for the differential problem in (1.1) defined in a polygonal domain. To this end, we propose a general use of the Schwarz-Christoffel mappings, that goes beyond the workhorse Laplace equation. To the best of our knowledge, this application of SC mappings is new.

We take advantage of the recently developed matrix equation framework for solving PDEs discretized using tensor bases ( [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF]), and bridge the gap between these bases and the treatment of more general, convex domains for which alternative discretization strategies are commonly employed. Our derivation discusses in detail some of the building blocks, the most crucial one being that the Jacobian determinant of the SC mapping depends on the reference domain variables, giving rise to more complex computations involving Hadamard products, if the tensor basis framework is to be exploited. We discuss how the numerical linear algebra needs to adapt to the considered operator, while pinpointing both memory and computational properties. In particular, both vector and matrix oriented equations are presented, and their pros and cons outlined. In our numerical experience, it is possible to take advantage at various stages of the properties of the involved operators, as opposed to other discretization strategies such as finite elements, where the differential operator and grid contributions cannot be distinguished.

The outline of the paper is as follows. After a brief introduction of the Schwarz-Christoffel conformal mapping in section 2, we explore the application of this mapping to the numerical discretization of PDEs on polygonal domains in section 3, which leads to a linear matrix equation framework. Then several specific examples of classes of PDEs are given in sections 4 and 5 to illustrate the derivation of the algebraic formulation. In section 6, various algebraic strategies are considered for solving the resulting matrix equation, and numerical experiments are discussed in section 7 to

show the performance of the discretization procedures and of the adopted solvers. We also analyze the selection of the grids in section 8. Finally, conclusions of this work are given in section 9.

All our numerical computations have been performed in Matlab [20].

Discretization and coordinate transformations. Discretization strate-

gies usually aim at transforming a physical region associated with a possibly complex geometry to a simple region through automatic transformations; these transformations can be constructed at the global level (such as in spectral methods) or at the local level (IGA, finite elements). These mappings allow one to then state the problem as an algebraic equation, whose solution approximates the given problem continuous solution in the original (physical) domain.

A large class of two-dimensional regions can be parameterized by means of special coordinate systems. This is the case for a circle and more generally, for a sector or an annulus, for which computational boundary-conforming grids such as polar coordinates, log-polar coordinates, parabolic cylinder coordinates, elliptic cylinder coordinates and bipolar coordinates transformations can be employed. It should be noticed, however, that if either the region or the boundary parameterization is changed, even slightly, these coordinates system may become useless.

This manuscript is for review purposes only.

As an illustrative example we consider the Poisson equation -∆u = f with zero boundary conditions, defined (in Cartesian coordinates) in the annulus sector Ω = {(x, y) : x = r cos θ, y = r sin θ, and r ∈ [r (0) , r (1) ], 0 ≤ r (0) < r (1) , θ ∈ [0, π 4 ]}. Using polar coordinates (r, θ), the Poisson equation is given by -r 2 u rr -u θθ -ru r = r 2 f (r, θ), (r, θ) ∈ [r (0) , r (1) ] × [0,

π 4 ].
The discretization of this equation in (r (0) , r (1) ) × (0, π 4 ) by FD gives the linear matrix equation (see, e.g., [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF])

Φ 2 T U + U T -ΦBU = F ,
where the matrices B, T ∈ R n×n account for the discretization coefficients of the first and second one-dimensional derivatives, respectively (see also section 3.1), while Φ = diag(r 1 , r 2 , . . . , r n ) and the r i 's are the interior discretization nodes of the interval [r (0) , r (1) ]. Here and in the following, for a vector d, we denote with diag(d) the diagonal matrix having the components of d on the diagonal.

Although a variety of special coordinate systems are available (see, e.g., [18, Chapter 1]) that can lead to the use of tensor-product based domains, this is not so for simple domains such as polygons with more than four edges. Moreover, it is important to have a computational procedure that automatically transforms a given two-dimensional domain into a rectangle. In the following we explore the use of the Schwarz-Christoffel mapping, a convenient conformal mapping whose computational realization is available via the Matlab software package sc [START_REF] Driscoll | Algorithm 756: A MATLAB Toolbox for Schwarz-Christoffel mapping[END_REF], [START_REF] Driscoll | Algorithm 843: Improvements to the Schwarz-Christoffel Toolbox for MAT-LAB[END_REF].

Other transformations could be considered to arrive at our algebraic framework, as long as the physical domain can be brought into a rectangular computational domain, on which the transformed problem can be solved; we refer the reader to, e.g., [START_REF] Knupp | The fundamentals of grid generation[END_REF] for a broad discussion of mappings associated with grid generations.

2.1. Schwarz-Christoffel mapping. Conformal mappings provide a mathematically robust tool for handling complex two-dimensional geometries, also allowing for an easy representation of the physical boundary conditions. From a numerical viewpoint, the resulting grids have the desirable property of being orthogonal, leading to reduced truncation errors when dealing with the transformed differential equations.

Moreover, the Jacobian is positive definite by construction, so that the mapping is guaranteed to be well defined (one-to-one), thus preserving the type of partial differential equation after transformation. A conformal mapping can be constructed algebraically or by solving either an integral equation [START_REF] Driscoll | Cambridge Monographs on Applied and Computational Mathematics[END_REF][START_REF] Fornberg | A numerical method for conformal mappings[END_REF][START_REF] Ives | Conformal grid generation[END_REF][START_REF] Symm | An integral equation method in conformal mapping[END_REF][START_REF] Thompson | Numerical Grid Generation: Foundations and Applications[END_REF], or (systems of) partial differential equations [START_REF] Chakravarthy | Numerical conformal mapping[END_REF][START_REF] Challis | A Numerical Method for Conformal Mapping[END_REF][START_REF] Mastin | Quasiconformal mappings and grid generation[END_REF]. The Schwarz-Christoffel mapping is one such operator, expressed as an integral [START_REF] Castillo | Mathematical Aspects of NUMERICAL grid generation[END_REF], [START_REF] Driscoll | Cambridge Monographs on Applied and Computational Mathematics[END_REF], mapping the complex upper half-plane into, e.g., the interior of a polygon. The "polygon" may in fact be a quite general region, and may have cracks or vertices at infinity.

Let the polygon vertices be denoted by z 1 , ..., z n , and let the numbers φ 1 π, ..., φ n π be the interior angles at the vertices. The pre-images of the vertices, or pre-vertices, are real and denoted by ω 1 , ..., ω n , which satisfy

ω 1 < ω 2 < • • • < ω n = ∞.
Then the Schwarz-Christoffel (SC) map g can be expressed in the form

(2.1) g(ω) = g(ω 0 ) + c ω ω0 n-1 j=1 (ζ -ω j ) φj -1 dζ,
This manuscript is for review purposes only.

known as the Schwarz-Christoffel formula.

There are two main difficulties associated with computing conformal maps from the Schwarz-Christoffel formula: i) determining the pre-vertices ω j , which is the first step for any SC mapping -the so-called SC parameter problem, and ii) integrating the right-hand side of (2.1). Except for special cases, the pre-vertices ω j cannot be determined in closed form; solving the system of equations for the pre-vertices can be a computationally expensive process. As of ii), finding an explicit expression by evaluating the integral analytically can significantly reduce the time needed to compute the map. Historically, these problems have been difficult to implement numerically, and first successful attempts have been reported in the early 1980s [START_REF] Trefethen | Numerical computation of the schwarz-christoffel transformation[END_REF], [START_REF] Ives | Conformal grid generation[END_REF]. In the past few years the Schwarz-Christoffel Toolbox for Matlab [START_REF] Driscoll | Algorithm 756: A MATLAB Toolbox for Schwarz-Christoffel mapping[END_REF][START_REF] Driscoll | Algorithm 843: Improvements to the Schwarz-Christoffel Toolbox for MAT-LAB[END_REF] has arisen as a numerical reliable tool for computationally handling these mappings; the toolbox can solve the parameter problem for various typical domains, and obtain the forward and inverse mappings. That provides us with a powerful and effective tool for solving numerically partial differential equations on a polygonal domain by the FD method, leading to matrix equations. Other software packages associated with conformal or other mappings could also be considered, whenever available [START_REF] Knupp | The fundamentals of grid generation[END_REF]. formation is an increase in the complexity of the transformed equations, but the domain becomes much simpler, the equations and the boundary conditions become easier to be approximated accurately by the FD method.

Let (x, y) ∈ Ω and (ξ, η) ∈ Π. Then identifying these two sets with corresponding regions in the complex plane, there holds

(3.1) z = g(ω) = g(ξ + iη) = x(ξ, η) + iy(ξ, η),
where x = x(ξ, η) and y = y(ξ, η) are real valued functions. Since g is a conformal mapping, the Cauchy-Riemann equations hold, that is

(3.2)
x ξ = y η and x η = -y ξ .

The functions u(x, y) and f (x, y) in Ω are transformed to functions u and f in Π using

u = u(ξ, η) = u(x(ξ, η), y(ξ, η)), f = f (ξ, η) = f (x(ξ, η), y(ξ, η)),
and similarly for α i (x, y), w i (i = 1, 2) and β(x, y). Without ambiguity, we use u, f , α i , w i and β to denote these transformed functions in the sequel.

The Jacobian matrix of the transformation g is

J = x ξ x η y ξ y η .
Thanks to (3.2) its determinant satisfies

J = J (ξ, η) = x ξ y η -x η y ξ = x 2 ξ + x 2 η > 0.
Therefore, the relations u ξ = u x x ξ + u y y ξ and u η = u x x η + u y y η yield (3.3)

u x = 1 J ( u ξ y η -u η y ξ ) = 1 J [( uy η ) ξ -( uy ξ ) η ] , u y = 1 J ( u η x ξ -u ξ x η ) = 1 J [( ux ξ ) η -( ux η ) ξ ] .
This manuscript is for review purposes only.

Combining with (3.2) gives (3.4)

u xx = 1 J   u ξ x 2 ξ + u η x ξ x η J ξ + u η x 2 η + u ξ x ξ x η J η   , u yy = 1 J   u ξ x 2 η -u η x ξ x η J ξ + u η x 2 ξ -u ξ x ξ x η J η   .
Using the formulas (3.2)-(3.4) and after simple calculations, the partial differential equation (1.1) on the physical domain Ω can be transformed into the following equation on the canonical domain Π,

(3.5) -( α 1 x 2 ξ + α 2 x 2 η ) u ξξ -( α 1 x 2 η + α 2 x 2 ξ ) u ηη -2( α 1 -α 2 )x ξ x η u ξη + J ( w 1 x ξ -w 2 x η ) - 1 J 2 ( α 1 -α 2 )(3x ξ x 2 η x ξξ + x 3 η x ξη -x 3 ξ x ξξ -3x 2 ξ x η x ξη ) u ξ + J ( w 1 x η + w 2 x ξ ) - 1 J 2 ( α 1 -α 2 )(x 3 η x ξξ + x 3 ξ x ξη -3x ξ x 2 η x ξη -3x 2 ξ x η x ξξ ) u η + J 2 β u = J 2 f , (ξ, η) ∈ Π,
or, with obvious notation,

(3.6) -γ 1 u ξξ -γ 2 u ηη + γ 3 u ξη + γ 4 u ξ + γ 5 u η + γ 6 u = f , (ξ, η) ∈ Π.
The transformed equation thus includes mixed second-order derivatives, in addition to the first and second order derivatives of the original problem. Clearly, the problem simplifies if, for instance, α 1 = α 2 . In particular, for α 1 = α 2 = 1 from (3.5) we obtain

-J u ξξ -J u ηη + J ( w 1 x ξ -w 2 x η ) u ξ + J ( w 1 x η + w 2 x ξ ) u η + J 2 β u = J 2 f ,
where a factor J can be eliminated throughout and the vector w highlighted, so that

(3.7) -u ξξ -u ηη -wJ ∇ u + J β u = J f , ∇ u = ( u ξ , u η ) T .
This simplified form stresses the role of the mapping on the various terms, elegantly and clearly discerning between the operator and domain discretization parts of the numerical procedure.

ξ i , η j ) ∈ Π = (a 1 , a 2 )×(b 1 , b 2 ) by h 1 = a 2 -a 1 n 1 + 1 , ξ i = a 1 + ih 1 , 1 ≤ i ≤ n 1 , h 2 = b 2 -b 1 n 2 + 1 , η j = b 1 + jh 2 , 1 ≤ j ≤ n 2 .
To avoid a pedantic notation, with some abuse of notation in the following we shall not distinguish between n i and n i -2, which corresponds to either include or exclude the interval boundary points.

This manuscript is for review purposes only.

The transformation x = x(ξ, η) and y = y(ξ, η) carries the canonical-space grid to a physical-space grid (x i,j , y i,j ), where

x i,j = x(ξ i , η j ), y i,j = y(ξ i , η j ), 1 ≤ i ≤ n 1 , 1 ≤ j ≤ n 2 .
One of the convenient features of these conformal mappings is that the derivatives

x η , x ξ can be cheaply and accurately computed by using the mapping itself [9, Table III]. In matrix notation, these will be denoted by

(X η ) ij ≈ x η (ξ i , η j ) and (X ξ ) ij ≈ x ξ (ξ i , η j ).
The discretization also induces discrete values for all of the transformed variables in (3.6). For example,

u i,j = u(ξ i , η j ) = u(x(ξ i , η j ), y(ξ i , η j )) = u(x i,j , y i,j ) = u i,j .
Let U i,j be the value of the approximation to u i,j . In order to attain second-order accuracy of the discretization of the two-dimensional problems (3.6), the following standard FD approximations are adopted for i = 1, ..., n 1 , j = 1, ..., n

u ξξ (ξ i , η j ) ≈ U i+1,j -2U i,j + U i-1,j h 2 1 , u ξη (ξ i , η j ) ≈ U i+1,j+1 -U i+1,j-1 -U i-1,j+1 + U i-1,j-1 4h 1 h 2 , u ηη (ξ i , η j ) ≈ U i,j+1 -2U i,j + U i,j-1 h 2 2 , u ξ (ξ i , η j ) ≈ U i+1,j -U i-1,j 2h 1 , u η (ξ i , η j ) ≈ U i,j+1 -U i,j-1 2h 2 . 2 , (3.8)
These approximations, performed element-wise in the grid, allow one to derive linear matrix equations in the unknown matrix U i,j , directly associated with the grid points (ξ i , η j ), and thus with the corresponding nodes (x i,j , y i,j ) in the physical domain. In particular, by collecting the coefficients in the linear combinations of adjacent nodes, we see that (see, e.g., [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF])

(3.9) u ξξ (ξ i , η j ) ≈ (T 1 U ) i,j , T 1 = 1 h 2 1        2 -1 -1 2 -1 . . . . . . . . . -1 2 -1 -1 2        ∈ R n1×n1 and (3.10) u ηη (ξ i , η j ) ≈ (U T 2 ) i,j , T 2 = 1 h 2 2        2 -1 -1 2 -1 . . . . . . . . . -1 2 -1 -1 2        ∈ R n2×n2 ,
A more detailed description of the other terms will be given in later sections.

In the following, several specific examples of classes of PDEs are given to illustrate the derivation of the actual matrix equation formulation.

This manuscript is for review purposes only.

-u xx -u yy = f, (x, y) ∈ Ω,
with Dirichlet boundary conditions. The Laplace operator is the most representative of the strength of Schwarz-Christoffel transformations, as the operator in invariant under conformal mappings; see, e.g., [START_REF] Driscoll | Cambridge Monographs on Applied and Computational Mathematics[END_REF], [START_REF] Chakravarthy | Numerical conformal mapping[END_REF], [START_REF] Bazant | Conformal mapping of some non-harmonic functions in transport theory[END_REF]. Hence, according to (3.7), the transformed equation associated with (4.1) is With the FD approximations (3.9)-(3.10), the discretized matrix form of equation (4.2) is given by

(4.2) -u ξξ -u ηη = J f , (ξ, η) ∈ Π. -3 -2 -1 0 1 2 3 -3 -2 -
(4.3) T 1 U + U T 2 = F, where (4.4) F = F + b.c., and
F i,j = (J f )(ξ i , η j ), 1 ≤ i ≤ n 1 , 1 ≤ j ≤ n 2 .
The linear equation (4.3) is a Sylvester matrix equation, and can be solved by the Bartels-Stewart method or its variants [START_REF] Bartels | Solution of the matrix equation AX + XB = C[END_REF][START_REF] Golub | Matrix Computations[END_REF][START_REF] Sorensen | Direct methods for matrix Sylvester and Lyapunov equations[END_REF], or by iterative methods [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF],

depending on the matrix size. We note that the right-hand side matrix F is computed element-wise as (J f )(ξ i , η j ). Let J be the matrix accounting for J (ξ i , η j ) for all interior grid nodes. The occurrence of element-wise computations with J provides an important challenge in the solution of the matrix equation form, which will become apparent when solving more involved equations in the next sections.

Let Ω be the interior of the hexagon in Using the SC Matlab Toolbox [START_REF] Driscoll | Algorithm 756: A MATLAB Toolbox for Schwarz-Christoffel mapping[END_REF], [START_REF] Driscoll | Algorithm 843: Improvements to the Schwarz-Christoffel Toolbox for MAT-LAB[END_REF], the obtained canonical (rectangular) domain is

given by Π = [-1.60167736, 1.60167736] × [0, 2.71248244]. A sample of the graphical representation of this grid mapped by g in Ω is also shown in Figure 1, while a finer grid (with n 1 = n 2 = 30 grid points in each direction in Π) is used in our experiment.

This manuscript is for review purposes only.

We solve the problem using two different settings: we first use homogeneous (zero)

Dirichlet boundary conditions, and then we use the following boundary conditions: The approximate solution U to (4.1) for f = 1 so that F = J under zero boundary conditions is displayed in Figure 2(left), while that associated with the non-uniform boundary conditions is reported in Figure 2(right). Due to the small problem size, the computational results were obtained with the Matlab function lyap [20].

5. Beyond the Laplace operator. In (1.1), the presence of space-dependent coefficients α 1 , α 2 , w and β provides an important challenge to the discretized matrix equation setting. If Ω were a rectangle, non-constant but separable coefficients would be easy to deal with. For the sake of the discussion, let us consider the zero-order term. Assuming for instance that β(x, y) = β 1 (x)β 2 (y), then the term βu would be discretized as β β β 1 Uβ β β 2 , where β β β 1 = diag(β 1 (x 1 ), . . . , β 1 (x n1 )) and analogously for β β β 2 ; see, e.g., [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF]. In the case of a polygon, equation (3.7) shows that the term βu is replaced by J βu, where all quantities depend on the space variables, thus the product needs to be carried out element by element, even in the case that β is separable, since β and J are not separable in general.

If G denotes the matrix associated with the discretization of J β, then the term J βu leads to a matrix term in the matrix equation of the form

G • U . It should be noticed, however that if G is a low rank matrix, that is it satisfies G ≈ ZQ T with Z = [z 1 , . . . , z ], Q = [q 1 ,
. . . , q ] and min{n 1 , n 2 }, then thanks to the properties of the Hadamard product, we could write G • U = diag(z 1 )U diag(q 1 ) + . . . + diag(z )U diag(q ). Hence, in the case G can be well approximated by a low rank matrix, the Hadamard product can be substituted with a sum of matrix terms in the usual matrix product.

In the next sections we dwell with the algebraic formulations for a few classes of differential problems.

This manuscript is for review purposes only. we can write

(5.2) -u ξξ -u ηη + (J β) u = J f , (ξ, η) ∈ Π, u = 0, (ξ, η) ∈ ∂Π.
Proceeding with the FD discretization, we denote1

(5.3) G 6 (i, j) := (J β)(ξ i , η j ), with 1 ≤ i ≤ n 1 and 1 ≤ j ≤ n 2 ,
thus giving the following matrix equation, corresponding to equation (5.2),

(5.4)

T 1 U + U T 2 + G 6 • U = F,
where • denotes the (element-wise) Hadamard product, T 1 , T 2 and U , F are defined in section 3.1.

Variable coefficient diffusion equation.

Consider the variable coefficient equation

(5.5) -α 1 u xx -α 2 u yy = f, (x, y) ∈ Ω, u = 0, (x, y) ∈ ∂Ω.
with α i = α i (x, y), i = 1, 2. According to (3.6), the transformed equation is

(5.6) -γ 1 u ξξ -γ 2 u ηη + γ 3 u ξη + γ 4 u ξ + γ 5 u η = f , (ξ, η) ∈ Π, u = 0, (ξ, η) ∈ ∂Π.
with (5.7)

γ 1 = α 1 x 2 ξ + α 2 x 2 η , γ 2 = α 1 x 2 η + α 2 x 2 ξ , γ 3 = 2( α 2 -α 1 )x ξ x η , γ 4 = 1 J ( α 2 -α 1 )(3x ξ x 2 η x ξξ + x 3 η x ξη -x 3 ξ x ξξ -3x 2 ξ x η x ξη ), γ 5 = 1 J ( α 2 -α 1 )(x 3 η x ξξ + x 3 ξ x ξη -3x ξ x 2 η x ξη -3x 2 ξ x η x ξξ ). and f = J 2 f . Given (5.8) B =        0 1 -1 0 1 . . . . . . . . . -1 0 1 -1 0        , we denote by B 1 := 1 2h1 B ∈ R n1×n1 , B 2 := 1 2h2 B T ∈ R n2×n2 the matrices associated
with the discretization of the first order derivative in ξ and η, respectively, and for = 1, . . . , 5 we set

(5.9) (G ) i,j = γ (ξ i , η j ), i = 1, . . . , n 1 , j = 1, . . . , n 2 .
With these notations, we get the following matrix equation for the problem (5.6)

(5.10) G 1 • (T 1 U ) + G 2 • (U T 2 ) + G 3 • (B 1 U B 2 ) + G 4 • (B 1 U ) + G 5 • (U B 2 ) = F,
where T 1 , T 2 and U , F are defined in section 3.1.

Convection-diffusion equation.

Consider the numerical solution of the partial differential equation

(5.11) -ε∆u + u x = f, (x, y) ∈ Ω, u = 0, (x, y) ∈ ∂Ω.
According to (3.7) the corresponding transformed equation is

(5.12) -ε( u ξξ + u ηη ) + x ξ u ξ + x η u η = J f , (ξ, η) ∈ Π, u = 0, (ξ, η) ∈ ∂Π.
Using the approximations in (3.8) gives the following matrix equation

(5.13) ε(T 1 U + U T 2 ) + X ξ • (B 1 U ) + X η • (U B 2 ) = F,
where T 1 , T 2 , U, F and B 1 , B 2 are defined in section 3.1 and after (5.8) respectively, while we recall that (X ξ ) i,j = x ξ (ξ i , η j ), (X η ) i,j = x η (ξ i , η j ), i = 1, . . . , n 1 , j = 1, . . . , n 2 .

6. Numerical solution of the matrix equation. Discretization of the transformed equation (3.5) by standard FDs leads to a medium scale linear matrix equation, This may be considered the major advantage over vector oriented discretizations. In fact, if the original PDE were to be solved by a method leading to a vector-based algebraic equation, one would be then faced with the solution of a linear system of size of the order of n 1 n 2 , for a comparable discretization accuracy. In that case too, a direct solution would be unfeasible. On the other hand, the presence of coefficient matrices involving the Hadamard product may make the computations more cumbersome.

The complexity of the solution strategy and its computational costs for solving discretized PDEs in the form (1.1) usually depend on whether the differential operator is diffusion dominated or convection dominated. In our matrix oriented context, we experienced a similar situation, therefore more than one method may be considered.

Fixed point iteration.

Let us write the left-hand side of matrix equation (5.10) as the sum of two matrix terms, so that

L(U ) -N (U ) = F,
where the choice of L depends on the actual continuous problem. Assuming L to be an invertible operator, and given a starting approximation U 0 , a sequence of approximates {U k } can be obtained by a stationary iteration

Solve L(U k+1 ) = F + N (U k ), k = 0, 1, . . .
This manuscript is for review purposes only.

or, alternatively,

U k+1 = U k +Z k , where Z k solves L(Z k ) = R k and R k = F -L(U k )+ N (U k )
is the current residual matrix. This approach is classical in the matrix equation context, and has led to ADI-type iterations already in the 1960s [START_REF] Wachspress | Iterative solution of elliptic systems[END_REF], [START_REF] Ellner | New ADI model problem applications[END_REF]; see also [START_REF] Shank | Efficient low-rank solutions of generalized Lyapunov equations[END_REF] for a more recent use in the multiterm matrix equation context.

The effectiveness of the approach depends on the cost of solving with the chosen operator L, and on the number of iterations needed to reach the required accuracy, which in turn depends on the spectral properties of the iteration operator, that is of

L -1 (N (•)).
In general, we found it hard to get convergence of this approach without any a-priori spectral knowledge of the operators involved.

6.2. Vector iterative solver with matrix-oriented preconditioning. A typically effective alternative to stationary iterations is the use of standard (vector)

iterative solvers. The discretization by FDs (or other methods) of (3.6) on the rectangle Π classically leads to a system of the form (see, e.g., [START_REF] Olshanskii | Iterative methods for linear systems, Theory and Applications[END_REF], [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF])

(6.1) Au = f ,
where the components of u collect the entries of U ij in lexicographic order, while A accounts for the discretization of all derivatives in the two space variables. In our context, we can derive the vector form by using the Kronecker operator, defined as (see, e.g.,[16, sec.1.3.6])

A ⊗ B =    A 1,1 B • • • A 1,n B . . . . . . . . . A n,1 B • • • A n,n B    .
This matrix product satisfies

(6.2) vec(AXB) = (B T ⊗ A)vec(X),
where vec(X) stacks the columns of X one below the other. We also recall that vec(G • M ) = diag(vec(G))vec(M ) for any matrices G, M of equal dimensions. The vectorization of (5.10) is thus given by (6.1) obtained as

D 1 (I ⊗ T 1 ) + D 2 (T T 2 ⊗ I) + D 3 (B T 2 ⊗ B 1 ) + D 4 (I ⊗ B 1 ) + D 5 (B T 2 ⊗ I) u = f,
where D i = diag(vec(G i )). The system to be solved has very large dimensions, with

A of size n 1 n 2 × n 1 n 2 , so that preconditioned iterative Krylov subspace methods are usually employed, such as CG, MINRES or GMRES [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. To limit memory consumptions for A nonsymmetric, a restarted version is usually employed, with a restart occurring every m iterations.

A large variety of preconditioners can be employed, such as incomplete LU, algebraic multigrid, and also operator strategies; see, e.g., [START_REF] Olshanskii | Iterative methods for linear systems, Theory and Applications[END_REF], [START_REF] Benzi | Preconditioning techniques for large linear systems: a survey[END_REF], [START_REF] Elman | Finite Elements and Fast Iterative Solvers, with applications in incompressible fluid dynamics[END_REF]. In the last class we can include operators that take into account the matrix equation structure of the problem, resulting from the performed transformation. This last approach has recently been used in algebraic equations stemming from different discretization techniques, see, e.g., [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF], [START_REF] Sangalli | Isogeometric preconditioners based on fast solvers for the Sylvester equation[END_REF]. The performance of the iterative scheme thus strongly depends on the effectiveness of the preconditioner. We shall see that operator based preconditioners can be quite robust, and may lead to mesh independent performance, in terms of number of iterations; see, e.g., [START_REF] Olshanskii | Iterative methods for linear systems, Theory and Applications[END_REF]Ch.4] for a general discussion on mesh independence of preconditioning strategies.

This manuscript is for review purposes only.

Galerkin projection method.

An approach that has recently shown its effectiveness in solving linear matrix equations with several terms consists of reducing the problem dimension by projecting the equation onto a smaller space.

To derive the method in its generality, we write the linear matrix equation as

G 1 • (T 1 U ) + G 2 • (U T 2 ) + G 3 • (B 1 U B 2 ) + G 4 • (B 1 U ) + G 5 • (U B 2 ) + G 6 • U = F, or, in short, S(U ) = F .
Let V and W be two subspaces of R n (here for simplicity we assume

n 1 = n 2 = n),
and let the k columns of V k (W k ) be orthonormal bases for V (W) with k n. Then we look for an approximation (6.3)

U (k) = V k Y k W T k ≈ U,
and let (6.4)

R k := F -S(U (k) )
be the associated residual. To determine Y k , an orthogonality (Galerkin) condition is imposed to the residual matrix. With the matrix inner product, this corresponds to imposing (6.5)

V T k R k W k = 0, that is V T k S(V k Y k W T k )W k = V T k F W k .
As an example, consider the equation (5.4), where G 1 = G 2 = 1 (the matrix of all ones) so that G 1 • (T 1 U ) = T 1 U and analogously for G 2 , while all other G i s are zero except G 6 . We thus obtain (6.6)

V T k T 1 V k Y k + Y k W T k T 2 W k + V T k (G 6 • (V k Y k W T k ))W k = V T k F W k .
This reduced dimension equation is solved by resorting to its vectorization by means of the Kronecker product in (6.2). Therefore, we can write u (k) := vec(U (k) ) = (W k ⊗ V k )vec(Y k ). For the example above, we thus have

(6.7) A k u (k) := (I ⊗ (V T k T 1 V k ))y + ((W T k T T 2 W k ) ⊗ I)y + (W k ⊗ V k ) T D 6 (W k ⊗ V k )y = f k ,
where D 6 = diag(vec(G 6 )) and f k = (W k ⊗ V k ) T vec(F ). This system can be solved by either direct or iterative methods. To make the whole procedure efficient, the cost of forming the third term should be decreased, as it involves matrices with dimension n vectors. The following results show how this computational cost can be significantly decreased.

Proposition 6.1. Let M, N ∈ R n×k1 , Z, S ∈ R n×k2 , and let D ∈ R n×n be diagonal. Define the k 1 k 2 × k 1 k 2 matrix (6.8) H := (M ⊗ Z) T diag(vec(D))(N ⊗ S) =       H 11 H 12 • • • H 1k1 H 21 H 22 H 2k1 . . . . . . . . . . . . . . . H k1k1       , H ij ∈ R k2×k2 .
Then the (i, j)-block of H is (6.9)

H ij = Z T DS, where D = diag (D(M :,i • N :,j )) ,
where M :,i (N :,j ) denotes the ith column of M (the jth column of N ).

This manuscript is for review purposes only.

Proof. Let D i = diag(D :,i ) ∈ R n×n . Then for the (i, j)-block of H, we have (6.10)

H ij = n l=1 M l,i Z T D l N l,j S = Z T n l=1 M l,i D l N l,j S = Z T diag D (M :,i • N :,j ) S. .
The next result describes how to update the matrix H when the building blocks increase their size.

Proposition 6.2. Assume the notation of Proposition 6.1 holds, and that the columns of the matrices M, N and Z, S are increased by one, i.e., their dimensions are n × (k 1 + 1) and n × (k 2 + 1), respectively.

(i) If i ≤ k 1 and j ≤ k 1 , then H ij = Z T k2 DS k2 Z T k2 DS :,k2+1 Z T :,k2+1 DS k2 Z T :,k2+1 DS :,k2+1
, where Z k2 collects the first k 2 columns of Z. In particular, the first block Z T k2 DS k2 is the matrix H ij in the previous step;

(ii)

If i = k 1 + 1 or j = k 1 + 1, then H ij = Z T DS.
The results in Propositions 6.1 and 6.2 give a simplified way to compute the quantity

V T k (G 6 • (V k Y k W T k ))W k and other similar quantities that arise in explicitly computing V T k S(V k Y k W T k )W k . In particular, using D 6 = diag(vec(G 6 )) we can write the (i, j)-block of (W k ⊗ V k ) T D 6 (W k ⊗ V k ) as [(W k ⊗ V k ) T D 6 (W k ⊗ V k )] ij = V T k diag (G 6 (W :,i • W :,j )) V k .
Here and in the following we denote by W :,i the ith column of the matrix W k , omitting the subscript k to lighten the notation.

The following result provides the application of the expressions above.

Proposition 6.3. Partition the matrix A k ∈ R k 2 ×k 2 in the equation (6.7) into k × k blocks with the size of each block being k × k, then the (i, j)-block of matrix

A k is A k ij = V T k GV k ,
where

G := diag (W :,i • W :,j ) T 1 + diag (W :,i • (T 2 W k ) :,j ) + diag (G 6 (W :,i • W :,j )) .
Note that the expression for A k ij can be given in a simplified formula as follows:

A k ij = δ ij V T k T 1 V k + (W T k T 2 W k ) :,j I + V T k diag (G 6 (W :,i • W :,j )) V k ,
where δ ij is the Kronecker delta function. Furthermore, when the dimension of the subspaces V and W is increased by one, the (i, j)-block of matrix A k+1 is

         A k+1 ij = A k ij V T k GV :,k+1 V T :,k+1 GV k V T :,k+1 GV :,k+1
, for i ≤ k and j ≤ k;

A k+1 ij = V T k+1 GV k+1 , for i = k + 1 or j = k + 1.
It was shown in [START_REF] Palitta | Optimality properties of Galerkin and Petrov-Galerkin methods for linear matrix equations[END_REF] that if the operator associated with the problem is symmetric and positive definite, then the Galerkin projection is optimal, in the sense that it minimizes the error in the corresponding norm. More precisely, the following result holds.

This manuscript is for review purposes only.

Proposition 6.4. ( [START_REF] Palitta | Optimality properties of Galerkin and Petrov-Galerkin methods for linear matrix equations[END_REF]) Let S(X) = F with S : X → j=1 A j XB j be a symmetric and positive definite operator 2 , with A j ∈ R n1×n1 , B j ∈ R n2×n2 . Let U be the exact solution to the problem S(U ) = F , and let range(V k ), range(W k ) be the constructed approximation spaces, so that

U k = V k Y k W T k is the Galerkin approximate solution. Then U -U k S = min Z=V k Y W T k Y ∈R k×k U -Z S ,
where the norm is defined as X 2 S = trace j=1 X T A j XB j .

In our setting the situation may be different whenever Hadamard product terms occur in the operator. Nonetheless, this is not the case for the constant coefficient

Poisson equation in (4.1), where indeed error minimization is achieved. This property will be used in the discussion of section 8.1.

In terms of memory requirements, the advantage of the matrix projection approach over vector methods is that only vectors of length n 1 and n 2 are stored, as opposed to vectors of length n 1 n 2 . Their number depends on the maximum dimensions of the spaces that need to be generated, and this number is hard to determine a-priori. Nonetheless, some intuition can be obtained by a preliminary experiment with a coarse grid, since we have experimental evidence that the space dimensions grow only sublinearly with the number of nodes in each direction.

6.4.

Choosing the approximation spaces. The choice of the approximation space is important for the overall effectiveness of the projection. To this end, we consider rational Krylov subspaces, which have been shown to be applicable to the case of multiterm linear matrix equations [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF]. In its general form, a rational Krylov space is defined as

RK k (A, B, V 0 ) = range([V 0 , (A + s 2 B) -1 V 0 , . . . , (A + s k B) -(k-1) V 0 ]),
where V 0 is a tall matrix with ˆ linearly independent columns, and A and B are matrices of equal size, compatible with those of V 0 ; The space dimension will thus be not greater than ˆ k. The parameters s j can be computed a-priori or dynamically during the iteration.

For equations with several terms, there is no obvious recipe on which matrices A, B should be used to build the left and right approximation spaces. The principle is that the generated space should well approximate the spectral properties of all coefficient matrices. In our case, except for the Hadamard product, the coefficient matrices stem from operators of different orders, of which the second order ones may be dominant.

We tested different options, and found that the choice B = I and A = L -1 j B j L -T j was the most appropriate, where T j = L j L T j is the Cholesky factorization of T j , j = 1, 2.

The presence of the Hadamard product makes any theoretical analysis very hard, since this product does not easily preserve spectral properties. A deeper study would be of great interest, but it is beyond the goals of this paper.

For selecting the shifts s j we followed the greedy procedure first proposed in [START_REF] Druskin | Adaptive rational Krylov subspaces for large-scale dynamical systems[END_REF] in the nonsymmetric case, and used the projected matrices

V T k (T 1 + B 1 )V k and W T k (T 2 + B 2 )
W k for that purpose. We refer the reader to [START_REF] Druskin | Adaptive rational Krylov subspaces for large-scale dynamical systems[END_REF], [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF] for additional details and a more general picture associated with the actual procedure. 2 The operator S is symmetric and positive definite if and only if its Kronecker form matrix

j=1 B T j ⊗ A j is.
This manuscript is for review purposes only.

In our context the selection of V 0 for each space was also to be made. Since the right-hand side F is not low rank, the truncated singular value decomposition (SVD) of F was performed, namely F ≈ P Σ Z T , where P , Z have columns each, while Σ = diag(σ 1 , . . . , σ ) and the σ j s are the first singular values of F , in decreasing order. Hence, the first ˆ ≤ columns of P (of Z ) were selected to start building V k (W k ), where ˆ is the first integer such that (6.11)

ˆ j=1 σ j j=1 σ j > 1 -tol,
where tol is the tolerance used in the stopping criterion of the projection method.

The method did not seem to be overly sensitive to the choice of ˆ : while changing the tolerance in (6.11) changed ˆ , the final subspace dimensions to reach convergence did not significantly differ.

7. Numerical illustration. In this section, we discuss some of our numerical experiments to illustrate the performance of the discretization procedures and of the adopted solvers. In the following we report the CPU time required to reach a desired accuracy, which is measured in terms of the relative residual

(7.1) Res := ||R k || F ||F || F ,
where || • || F is the Frobenius norm. For both theoretical purposes and approximation quality evaluation, the relative error norm is also of interest,

(7.2) Err := ||U k -U || F ||U || F ,
where U is the exact algebraic solution, obtained by solving the discretized problem by a direct solver, while U k is the current approximation.

For simplicity of exposition, in all reported experiments we take n 1 = n 2 = n and f (x, y) ≡ 1, and except where explicitly stated, we consider two different physical domains: a convex polygon and an L-shaped domain (cf. Figure 3). The reported numerical results were obtained with the choice of discretization grid that gave us consistently the best results. A detailed discussion on the grid selection is postponed to section 8.

This manuscript is for review purposes only.

A few comments on the implementation details of the considered algorithms are in order.

The stationary iteration requires a Lyapunov solver at each iteration. Whenever the coefficient matrices are the same, these are eigendecomposed once for all, and the Lyapunov solution is explicitly written down (see, e.g., [27, section 4]). Otherwise, the call to the Matlab function lyap is performed. The residual matrix is computed explicitly. All these computations are affordable for n = O(10 d ) with d moderate. This is not restrictive for a two-dimensional problem, and a good accuracy of the sought after solution is obtained already for n ≈ 1000.

Preconditioned restarted GMRES(m) is restarted every m = 20 iterations, as commonly done [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Different preconditioning strategies could be employed. For simplicity, here we consider two alternatives: i) incomplete LU preconditioning 3 , after permutation of the coefficient matrix rows and columns with symmetric minimum degree ordering, ii) operator preconditioning, which consists of solving a Lyapunov equation with a portion of the coefficient matrix, typically the terms associated with the second order operator [START_REF] Palitta | Matrix-equation-based strategies for convection-diffusion equations[END_REF].

direct
-∆u + u = f, (x, y) ∈ Ω,
with f identically one. The discretized problem takes the form in (5.4), with G 6 = J and F = J • 1 = J, where we recall that the (i, j)th entry of J corresponds to

J (ξ i , η j ).
Table 1 shows the results of solving (5.4) on the L-shaped domain, for an SC grid taking as reference vertices the nodes [1 2 4 6] (ordered counterclock-wise), where the first node is (3.2010, -2.7532). All iterations also checked the true error, and the method stopped with an error Err below 10 -5 . In general, the relative residual norm in (7.1) was smaller or significantly smaller than the relative error norm in (7.2).

As expected, the direct method (Matlab backslash [20]) becomes more and more expensive as n grows, having to solve an n 2 × n 2 system. Here and later on, we stress that comparing with the direct method is unfair, as the direct method employs a builtin Matlab function at its highest efficiency, whereas all other methods mainly rely on interpreted commands in Matlab. Comparisons between the stationary iteration and restarted preconditioned GMRES is more appropriate. The number of GMRES iterations in the last restart is reported in parenthesis. For this simple operator, taking 3 ILU factorization with threshold and pivoting was employed, with droptol tolerance 10 -3 .

This manuscript is for review purposes only.

the second-order operator as preconditioner is clearly an advantage over a classical incomplete LU preconditioning, and this is also shown by the good performance of the stationary iteration, which uses the same operator. We also observe that the number of iterations for these two methods is independent of the number of grid nodes. This is related to the fact that the iteration operator (or preconditioned operator) is given (in Kronecker form) by L -1 (L + G) = I + L -1 G where L accounts for the Laplace operator while the diagonal matrix G collects the Jacobian determinant entries. For a dominant Laplace operator, the eigenvalues of I + L -1 G tend to tightly cluster, so that convergence of Krylov subspace methods can be faster [START_REF] Olshanskii | Iterative methods for linear systems, Theory and Applications[END_REF]Ch.4].

We ran the methods for the convex domain as well, with the direct and GMRES methods behaving the same as for the L-shaped domain. On the other hand, the stationary iteration failed to converge. This appears to be due to the different deter- 3. As before, timings for the direct methods should not be compared with those of the other non-built-in algorithms. Note that GMRES preconditioned by the Lyapunov operator in this case is not effective, and its takes too long for the largest grid; it was thus stopped prematurely. This could be predicted as in this case the non-constant operator coefficients are not well approximated by the unit constant value. The projection method behaves quite well, using significantly lower memory than GMRES with a quite dense preconditioner.

A convection-diffusion problem.

We consider the PDE

(7.4) -u xx -u yy + w • ∇u = f, (x, y) ∈ Ω,
with w = (ω, 0), on the same polygonal domain considered earlier, using t = [1 2 3 5] as vertices. Discretization leads to the matrix equation in (5.13).

Here we analyze how the performance of the iterative algebraic solvers varies as ω increases, so as to make the problem more convection dominated. on number of iterations, rather than on CPU times. The results are reported in Table 4. We clearly see that GMRES (20) with the operator preconditioner is insensitive to the meshsize variation, and it is only marginally affected by the increase of ω. ILUT preconditioning does not share similar properties in terms of meshsize, while similar results can be observed as ω varies. These results are similar to those observed when using operator preconditioners with other discretization approaches. It is interesting, however, that here the conformal mapping affects the convection terms in a different way than in other discretizations, as the partial derivatives X ξ and X η explicitly arise.

The projection method seems to be the most sensitive method to both meshsize and convection. On the other hand, it should be noticed that the projection method requires O(n • (spacedim)) memory allocations, which in the worst case amounts to O(1024 • 87), whereas GMRES (20) requires O(n 2 • 20) allocations, which goes up to O(1024 2 • 20) allocations for the largest problem considered, without taking into account the memory required for the factors in the ILUT preconditioning. Hence, if memory is a constraint, GMRES methods on the vector problem can be largely penalized.

To complete this experimental analysis, we consider a different non-convex domain, displayed in Figure 4, for the vertex choice t=[1 4 6 8] (the first vertex is the This manuscript is for review purposes only.

Table 5 Performance of iterative solvers on the convection-diffusion problem in (7.4) for ω = 1, with the non-convex domain in Figure 4.

most right-south one). Schwarz-Christoffel mappings are particularly appropriate to handle this type of complex domains. Table 5 reports the results of our experiments for problem (7.4) and ω = 1, to be compared with those in Table 4 for the same choice of ω. The performance of all methods is completely analogous to that for the polygonal convex domain, emphasizing that once a good vertex choice is made, the performance of this discretization procedure is algebraically robust. Indeed, for the projection method we found the selection of the vertices quite crucial, and this is further discussed in the next section.

8. Analysis of the grid selection. In our numerical experiments with the Galerkin method we found that the method was overly sensitive to the grid choice.

In particular, the number of iterations, and thus the space dimension of V k and W k , varied quite significantly for different vertex selections. We quote a paragraph from Trefethen and Driscoll's book [10, page 107ff] as a warning for such an occurrence:

An inspection of many of the figures in this book makes it clear that the size of equally spaced cells in the computational domain can vary greatly and rapidly in the physical domain [..]. This variation makes a certain amount of physical sense, at least in some contexts, but is undesirable in many applications. A particular problem [...] in SC mapping is the presence of singularities in f due to the corners. As a mapping, the SC formula handles corners elegantly and completely, but here they remain a significant challenge in the computational domain.

Before we continue we wish to stress that we did not encounter any performance This manuscript is for review purposes only. sensitivity in using the vector oriented GMRES. Hence, we conjecture that the computational sensitivity arises when the matrix properties of J play a role, whereas the vectorized version of J only influences the computation with its positive entries.

We thus proceed with a deeper analysis of the role of J in the algorithmic performance of the Galerkin approach for the convex polygon. A corresponding analysis was performed for the other considered domains with similar results, hence these are not reported. We focus our analysis on the Poisson equation in (4.2), where J only enters in the right-hand side and it is easier to analyze. Moreover, we distinguish between a good grid, for which we experienced fast convergence, and a bad grid, which led to a significant delay as the discretization is refined.

With the good grid t=[1 2 3 5], for instance, the Galerkin method applied with n = 2 10 selects an approximation to the right-hand side of rank 17, and with a space of dimension 45 reaches the accuracy required by (7.1). On the contrary, with the bad grid t=[1 2 3 6], in spite of a slightly smaller rank for the right-hand side approximation (rank 11), the method does not reach completion within a space of dimension 100. We stress that the singular values of J in the two cases show similar patterns and they did not appear to play any role in the different behavior. Figure 5 displays the residual norm history as the space dimension grows, for the two cases, illustrating the significantly different slopes in the convergence curve. The plot also includes the error F-norm history (see (7.2)) for the two cases. We recall that the Galerkin method minimizes the A-norm of the error for the given approximation space. Thus, the error behavior displayed in the plot, though in a different norm, indicates that the space generated by starting with vectors from the bad grid does not provide as good information as for the good grid. In passing, we also observe that the error is in both cases much smaller than the residual norm would predict. This appears to be unrelated to the low rank truncation of the right-hand side used to build the approximation space, but rather to classical results on the discrepancy between error and residual norms. For this specific problem, this is given by

R k F ≤ I ⊗ T 1 + T 2 ⊗ I U k -U F .
Since T 1 = T 2 and they are tridiagonal, symmetric and Toeplitz matrices of size n as defined in section 3.1, we get

I ⊗ T 1 + T 2 ⊗ I = λ max (I ⊗ T 1 + T 2 ⊗ I) = 2λ max (T 1 ) = 2 h 2 1 2 + 2 cos( π (n -1)
) .

This manuscript is for review purposes only. For n = 2 10 , it thus holds that I ⊗ T 1 + T 2 ⊗ I ≈ 8 • 10 6 , which largely accommodates the gap between error and residual norms observed in the plot.

8.1. Action of good and bad grids. Table 6 reports some grid information for different selections of the polygon vertices that are used by the mapping g to construct the four rectangle pre-vertices. More precisely, for each vertex selection, the aspect ratio is the ratio between the largest and smallest rectangle edges; max df is the maximum absolute value of the Jacobian determinant square root, max ξ,η x 2 ξ + x 2 η for the adopted grid; maxmin ratio is the ratio between the longest and shortest grid edges on the polygon boundary for the chosen grid. The first three grid selections in the table lead to good performance of the Galerkin method, the second three choices lead to a delayed convergence. The only parameter that seems to be in favor of the good grids is the maxmin ratio, which tends to be smaller than for the bad grids.

Figure 6 shows the grid mapped by the discretization of the corresponding triangle via the SC Toolbox, for two choices of vertices (a good grid on the left, a bad grid on the right). The Schwarz-Christoffel mapping associated with these two grids leads to the Jacobian determinant J = J (ξ i , η j ), i, j = 1, . . . , n in Figure 7. The peaks refer to singularity approximations corresponding to the disregarded polygon vertices. For the good choice t=[1 2 3 5] the missing vertices [4 6] are not consecutive, therefore each corresponds to a peak on a different axis direction. Two close peaks on a single edge occur in the bad choice t=[1 2 3 6], corresponding to the disregarded vertices [4 5]. We believe this is a key fact: alternating disregarded vertices lead to a better This manuscript is for review purposes only. approximation of the solution for this particular problem. We argue that this different behavior may affect the choice of the starting blocks generating the right and left approximation spaces, V k , W k , which in turn may lead to the generation of poor spaces. We start by recalling that if T 1 = QΛQ T is the eigenvalue decomposition of T 1 with Λ = diag(λ 1 , . . . , λ n ), then (8.1)

U = Q (Q T JQ) (Λ ⊗ I + I ⊗ Λ) Q T ,
where denotes elementwise division. For this problem the solution nonsymmetry is only due to J. If J were symmetric, then in the approximation V k Y k W T k ≈ U , the spaces V k , W k could be taken to be equal. Moreover, due to the expression in (8.1), a good approximation to U can be obtained only if the eigencomponents corresponding to the smallest eigenvalues of T 1 are well approximated in both spaces V k , W k . Since these eigencomponents are discretizations of smooth functions, it would be desirable that this type of vectors could be readily detected in the spaces, hence the two spaces should behave somewhat similarly. Let us see how these properties transfer into

This manuscript is for review purposes only.

spectral requirements on J. Let J = P ΣZ T be the singular value decomposition of J. Some sort of symmetry in J is reflected in the left and right singular vectors of J behaving similarly. This implies that taking the first left (resp. right) singular vectors to build V k (resp. W k ) then the generated spaces would be similar. This behavior is displayed in Figure 8 for J stemming from the good grid (left) and from the bad grid (right). If J were symmetric, then P T Z = I. This is not the case for either grid, however the plots show that for the good grid the neighboring vectors in P are close to the corresponding vectors in Z; this is not the case for the bad grid, as vectors in P have large components onto many column vectors in Z. By taking the first columns of P and Z to build the two approximation spaces, the resulting spaces will significantly differ for the bad grid.

The presence of peaks in J results into a corresponding structure in the vectors of P and Z. This is expected since the first singular triplets of J provide the best approximation to J of rank ([16, Theorem 2.4.8]). Hence, if the first columns of P and Z have a non-smooth behavior (to match the occurrence of several peaks), these will not be able to quickly approximate the smooth eigencomponents of the operator T 1 , leading to slow convergence. In other words, the smoothness of the operator T 1 should be well captured by the starting vectors generating V k , W k , for a fast approximation to take place. The bad grid (right plot of Figure 7) has several peaks on one side of the grid, and these are only taken care of by the left singular vectors. This seems to imply that V k and W k will be very different, and that the space W k will only slowly approximate the leading eigencomponents in the solution U , causing a delay.

In summary, our interpretation is that when choosing the vertices it is advisable to avoid unbalance by missing too many consecutive vertices together with keeping neighboring vertices, especially those corresponding to large values of the Jacobian. 9. Conclusions. We have analyzed the matrix computational problems associated with the discretization of linear elliptic partial differential equations of polygonal domains by means of the Schwarz-Christoffel mappings. The matrix oriented formulation leads to special linear matrix equations, involving Hadamard products, which make the use of matrix equation solvers quite challenging. On the other hand, iterative vector based methods such as GMRES preconditioned by the higher-order operator seem to behave quite well, with a number of iterations that seems to be mesh-independent, and quite insensitive to possible corner singularities of the mapping. If memory allocations are an issue, then matrix-oriented methods should be preferred: a working strategy with several implementation details is provided.

Our use of the Schwarz-Christoffel mappings was surely encouraged by the presence of an easy to use Matlab Toolbox. However, it is well known that these mappings are restricted to two-dimensional domains. Hence, the generalization of our methodology to three-dimensions would certainly need to significantly extend this setting, or employ different transformations for the generation of structured grids.

3 .

 3 Schwarz-Christoffel mappings for PDEs in polygonal domains. The Schwarz-Christoffel conformal function g maps a rectangular domain Π = (a 1 , a 2 ) × (b 1 , b 2 ) (canonical domain) onto the physical domain Ω. The cost of using such trans-

3. 1 .

 1 Discretization of the transformed PDE. To approximate the transformed equation (3.6) by FDs, we define the grid points (

4 .

 4 Poisson equation. The simplest possible case to be considered is the follow-

Fig. 1 .

 1 Fig. 1. Convex domain. Grid of the Schwarz-Christoffel mapping over the physical domain.

Figure 1

 1 having vertices (ordered counterclock-wise) z 1 = (2.06107, -2.59033), z 2 = (2.91603, 0.38168), z 3 = (1.22646, 2.58015), z 4 = (-1.90840, 2.47837), z 5 = (-2.92621, 0.58524) and z 6 = (-1.19593, -2.30534).

(4. 5 )

 5 u = (x -2.06107)(x -1.22646) + (y + 2.59033)(y -2.58015) + 1, on edges z 1 -z 2 and z 2 -z 3 u = 1 on other edges.

Fig. 2 .

 2 Fig. 2. Poisson equation solution U in the physical domain Ω with homogeneous (left) and non-constant (right) Dirichlet boundary conditions.

5. 1 .

 1 The Poisson equation with a reaction term. Consider the numerical solution of (5.1)-∆u + βu = f, (x, y) ∈ Ω, u = 0, (x, y) ∈ ∂Ω.The addition of the quantity βu, with β = β(x, y), increases the difficulty of the algebraic setting, since it adds a composite term to the matrix equation. Using(3.7)

Fig. 3 .

 3 Fig. 3. The physical domain Ω. Left: convex polygon. Right: L-shaped domain.

7. 1 .

 1 A reaction-diffusion problem. Consider the problem discussed in section 5.1 with β = 1, that is

7. 2 .

 2 A simple elliptic problem. Consider the partial differential equation (7.3) -(x 2 + 1)u xx -(y 2 + 1)u yy = 1, (x, y) ∈ Ω, u = 0, (x, y) ∈ ∂Ω, in the above two physical domains. The equation has non-constant coefficients, leading to the presence of the Jacobian determinant in the discretized matrix equation, as shown in (3.6). The results obtained with our experiments on the L-shaped domain are reported in Table

Fig. 4 .

 4 Fig. 4. General non-convex polygon. The physical domain Ω and an SC grid.

Fig. 5 .

 5 Fig. 5. Relative residual norm convergence for the good grid t=[1 2 3 5] and for the bad grid t=[1 2 3 6].

Fig. 6 .

 6 Fig. 6. Typical meshes for the polygonal domain (coarse grids are shown). Left: Good grid t=[1 2 3 5]. Right: Bad grid t=[1 2 3 6].

Fig. 7 .

 7 Fig. 7. Entry magnitudes of J on the (ξ, η) domain (log scale). Left: Good grid t=[1 2 3 5]. Right: Bad grid t=[1 2 3 6].

Fig. 8 .

 8 Fig. 8. Contour plot of entry magnitude of |P T Z| above the value 0.1, for J = P ΣZ T the SVD of J. Left: Good grid t=[1 2 3 5]. Right: Bad grid t=[1 2 3 6].

Table 1

 1 Example of section 7.1 on L-shaped domain.

		stationary	ILUT	Lyap
		iteration	GMRES(20)	GMRES(20)
	n	# iter	time # rest	time # rest	time
	64 0.007	20 0.016	1(6)	0.106	1(4) 0.033
	128 0.047	19 0.027	1(11)	0.107	1(5) 0.027
	256 0.183	19 0.121	2(1)	0.552	1(5) 0.124
	512 0.677	18 0.617	3(12)	4.741	1(5) 0.683
	1024 4.348	18 4.481	6(19) 53.686	1(5) 4.666

Table 2

 2 minant, J , whose discretized version J provides different spectral properties of the iteration matrix in this case. Example of section 7.2 on polygon, with nodes [1 23 5].

		direct	Projection	ILUT	Lyap
			method	GMRES(20)	GMRES(20)
	n		space dim	time # rest	time # rest	time
	64	0.056	29	1.320	1(9)	0.073	3(16)	0.079
	128	0.114	42	5.867	1(17)	0.166	4(17)	0.340
	256	0.519	46	10.481	2(20)	1.066	6(12)	1.924
	512	1.995	59	28.983	5(17)	13.927	8(18)	17.787
	1024 11.738	66	82.368	14(2) 119.030 11(18) 135.400

Table 3

 3 To this end, we focusThis manuscript is for review purposes only. Example of section 7.2 on L-shaped domain, with nodes [1 24 6].

		direct	Projection		ILUT	Lyap
				method		GMRES(20)	GMRES(20)
	n		space dim	time # rest	time # rest	time
	64	0.019	13	0.241	1(9)	0.083	5(17)	0.146
	128	0.097	23	0.698	1(16)	0.153 10(20)	1.291
	256	0.467	33	2.668	2(18)	0.868 25(12)	7.634
	512	1.995	50	13.700	4(13)	6.917 33(18) 77.751
	1024 10.334	70	100.430	9(2) 72.063	-(-)	-
			Projection ILUT	Lyap
				method	GMRES(20) GMRES(20)
		ω	n space dim #rest(its)	#rest(its)
			128	28	1(13)	1(9)
		1	256	36	2(13)	1(9)
			512	42	5(6)		1(9)
			1024	45	9(16)	1(9)
			128	40	1(9)		2(19)
		10	256	48	1(17)	2(19)
			512	54	7(19)	2(19)
			1024	67	14(4)	2(19)
			128	52	1(6)		6(8)
		30	256	62	1(11)	6(7)
			512	80	2(5)		6(6)
			1024	87	9(14)	6(6)

Table 4

 4 Performance of iterative solvers on the convection-diffusion problem in (7.4) for the convex domain.

Table 6

 6 Some sample grids discretizing the polygon, using n = 2 10 grid nodes on each side of the rectangle in (ξ, η).

	goodness	vertices	aspect ratio max df maxmin ratio
		[1 2 3 5]	1.10	16.33	389.4
		[1 3 5 6]	1.18	14.71	258.8
		[1 2 4 6]	1.06	12.30	106.7
		[1 2 3 6]	1.17	19.60	413.9
		[3 4 5 6]	1.14	25.21	769.1
		[1 2 3 4]	1.31	15.44	678.0

The numbering in G 6 accounts for more terms that will arise in the following, associated with the discretization of higher order terms.This manuscript is for review purposes only.

This manuscript is for review purposes only.

Acknowledgments. The second author is a member of Indam-GNCS. Its support is gratefully acknowledged. Part of this work was also supported by the Grant AlmaIdea 2017-2020 -Università di Bologna. The first author is funded by the China Scholarship Council (Contract No. 201906180033) and by the National Natural Science Foundation of China (Grant Nos. 11471150, 11401281). Her work was performed during her visit at the Università di Bologna, Italy.

This manuscript is for review purposes only.