Fractional partial differential equations (PDEs) have recently received a tremendous amount of attention, which can be attributed to the flexibility of fractional operators in capturing long-range effects, due to their nonlocal nature. In addition, they have fewer regularity requirements than their classical counterparts. In particular, the fractional Laplacian has been successfully used as a regularizer in imaging science in place of the total variation regularization [3,4]. Moreover, the fractional Helmholtz equation was derived in [32], using first principle arguments combined with a constitutive relation, to model geophysical electromagnetism. Other applications include: Quasigeostrophic flow [13], phase field models [2,3], porous media [14], and quantum mechanics [20]. Motivated by these applications, we introduce a new approach to solve the fractional PDE (-∆) s u = f in Ω u = 0 on ∂Ω,
where Ω ⊂ R n is a bounded open domain with boundary ∂Ω, and s ∈ (0, 1) is the fractional exponent. The operator (-∆) s denotes the spectral fractional Laplacian, whose rigorous definition will be provided in Section 2. The nonlocality of the fractional Laplacian makes it challenging to realize in practice [11,29]. However, several approaches exist. For example, the authors in [28] use a spectral discretization in space and discuss computing the spectrum of the Laplacian to realize the fractional Laplacian. Unfornately, computing the spectrum of an operator in general domains can each piece independently. As a result, this solver is easily parallelized. The convergence analysis for ultraspherical spectral methods can be found in [26]. In practice, we use the polynomial coefficients of the solution to determine if the solver has converged. Specifically, when a coefficient falls below a given threshold (e.g., we choose 10 -10 in this paper), we terminate the spectral method and use the result as the discrete approximation of the solution to the fractional PDE. In all numerical examples and for all values of s, we observe exponential convergence with just a few degrees of freedom. In this sense, the potential benefits of our method are clear. Additionally, the proposed approach has several other advantages over the existing spectral methods. We can combine our spectral methods on rectangles with the ultraspherical spectral element method (UltraSEM) [17] to develop solvers on polygonal domains with unstructured quadrilateral or triangular meshes. To be specific, we can partition the general domain Ω into quadrilateral subdomains, convert them into rectangles with a change of variables, and apply our solver. As a result, the total number of degrees of freedom is the sum of those used on each subdomain. Finally, when n = 2, we emphasize that our method can be generalized to solve fractional PDEs of the form:
where the operator L is the general elliptic operator Lu = -∇ • (A∇u) + cu. Here, A ∈ R 2×2 is matrix function that is symmetric and positive definite and 0 ≤ c ∈ R is a function. This is feasible due to [17]. However, [17] only considers the case with n = 2, motivating our restriction to n = 2 in (1.2). The remainder of the paper is organized as follows. In Section 2, we review preliminary results on fractional PDEs as well as the polynomial and tensor notation that we use in the subsequent sections. In Section 3, we introduce the direct solver for the disk domain. In Section 4, we solve the extended PDE on rectangles both directly and with a parallel domain-decomposition solver. Finally, in Section 5, we demonstrate our method on two applications: solving the fractional elliptic PDE on the cube, and solving a fractional PDE-constrained optimal control problem.
Let Ω ⊂ R n be a bounded open set with Lipschitz boundary ∂Ω. Let -∆ be the L 2 (Ω) realization of the standard Laplace operator with zero Dirichlet boundary conditions. It then follows that -∆ has compact resolvent and its eigenvalues can be arranged as 0
(Ω) the eigenfunctions corresponding to λ k . These eigenfunctions form an orthonormal basis of L 2 (Ω).
For s ≥ 0, we let
For a relation between H s (Ω) and the classical fractional-order Sobolev space H s (Ω), we refer to [5,15], and the references therein. We shall denote the dual of H s (Ω) by H -s (Ω). Specifically in this manuscript, we focus on 0 < s < 1. Now, to define the fractional Laplacian: (-∆) s is the mapping
To solve (1.1), we use the Caffarelli-Silvestre extension [11,29]. This requires introducing an extension variable ζ, and the change of variable z = ζ 2s 2s
. Then, the resulting problem aims to find U : Ω × [0, ∞) → R that satisfies the following equation
where α = 2 -1 s , and d s = s 2s-1 Γ(1-s) Γ(s) . Here, ∆ x denotes the Laplacian with respect to the original domain Ω and U zz denotes the second derivative with respect to the extended dimension z. After solving for U , we can recover the solution to (1.1) as u(x) = U (x, 0).
In general, approximating functions by polynomials on an unbounded domain is a challenging problem. Motivated by the fact that the solution U in the z-direction decays exponentially [24] (also confirmed by our numerical experiments), we consider the following truncated problem:
where R > 0 is the truncation parameter. In our numerical experiments, the choice of R is motivated by [24]. In particular, we set R = O(log(DoF Ω )) for the rectangular domains, where DoF Ω is the total degrees of freedom used for Ω. Experimentally, we notice that for the disc domain it is more appropriate to choose R = O(DoF
We emphasize that the additional variable z introduced by the extension requires that we solve a problem of one dimension higher. In particular, although Ω is chosen to be rectangles or disks in this paper, we must solve (2.2) in hexahedron or cylinders. In the subsequent sections, we review some basic polynomial bases for discretization and tensor operations for the discrete operators.
2.1. Ultraspherical Polynomial Basis and Spectral Methods. Ultraspherical (or Gegenbauer) polynomials are a special family of polynomials that are usually denoted by C (λ) n (x). Here, λ > 0 is a coefficient and n is the polynomial degree of x [25,Table 18.3.1]. They are orthogonal on the interval (-1, 1) with respect to the weight function w(x) = (1 -x 2 ) λ-1/2 and satisfy the three-term recurrence [25, Table 18.9.1]
(2.3)
For notational convenience, we use C(λ) n (x) to denote the L 2 (-1, 1)-normalized ultraspherical polynomials with respect to the weight w(x). Two well-known classes of polynomials, Chebyshev polynomials of the second kind and Legendre polynomials, are both special cases of ultraspherical polynomials, with coefficient λ = 1 and λ = 1/2, respectively. 2.2. Tensor Notations. We follow the notation for tensors found in [19], which we briefly review now for the reader's convenience. A tensor is a multidimensional array. We focus on threedimensional tensors, denoted by calligraphic upper case letters, such as X ∈ C n 1 ×n 2 ×n 3 . Comparatively, matrices are represented by upper case letters, such as X ∈ C n 1 ×n 2 . We emphasize that the following discussion directly extends to four-dimensional tensors, see Section 5.
To represent submatrices, and fibers and slices of tensors, we employ MATLAB notation. In particular, "m : n" means that we take all integers between m and n, including m and n, and a simple ":" means that we take all available indices. For example, A 1:4,: is a submatrix of A containing the first four rows and all the columns. We also use the keyword "end" to indicate the last index. For example, A end,: represents the last row of A. Fibers are higher-order analogues of matrix rows and columns, formed by fixing all but one index. A three-dimensional tensor X has column fibers X :,j,k , row fibers X i,:,k , and tube fibers X i,j,: . Slices are two-dimensional sections of tensors made by fixing all but two indices. A three-dimensional tensor X has horizontal slices X i,:,: , lateral slices X :,j,: , and frontal slices X :,:,k . The k-fold (or k-mode) product of a tensor X ∈ C n 1 ×n 2 ×n 3 with a matrix A ∈ C n k ×n k is denoted by X × k A, and defined elementwise as
(2.4)
The k-fold product corresponds to each mode-k fiber of X being multiplied by the matrix A.
In this section, we solve (2.2) on a disk domain. Without loss of generality, we let Ω be the unit circle. Otherwise, we can easily convert to this problem by scaling with the radius. We use polar coordinates to rewrite (2.2) as
To avoid the singularity at ρ = 0, we use the DFS method [33] to extend to ρ ∈ (-1, 1) by setting
Notice that both Ũ and f are now continuous at ρ = 0, which leads to simpler spectral discretization with smaller polynomial degrees. From this point, we work directly with these "doubled" functions so that the singularity at ρ = 0 does not require additional consideration. The disk domain allows us to assume that both the solution Ũ and the function f have Fourier expansions:
Ũk (ρ, z)e ikθ and f (ρ, θ)
In this way, we can decouple (3.1) into differential equations for each Fourier mode:
Following [18, § 4.1.2], we assume that the ansatz for Ũk is given by
where the term (1 -ρ 2 ) incorporates the boundary condition Ũk (±1, z) = 0. Then, we solve for Ṽk . The choice of the above ansatz only imposes partial regularity on Ũk , and we refer to [18, § 4.1.2] for a detailed discussion on why it is challenging to impose full regularity.
Since the fractional exponent s can be any number between 0 and 1, the function z 1/s makes the development of solvers for (3.2) a challenging task for our spectral method. To overcome this difficulty, we develop different solvers for varied values of s.
Polynomial Approximation of z 1/s . We first consider the case in which s ∈ (0, 1) is such that the function z → z (1/s) can be approximated accurately by a polynomial of low degree. For the discretized problem, the multiplication by z (1/s) is transformed into a matrix-matrix product. Consequently, the polynomial degree is directly related to the discretization size, and the definition of "low degree" in the above statement is related to the size of the discretized system one is capable of solving. For example, one may interpret "low degree" to mean "degree less than 50". In this case, any value of s such that 1/s is an integer between 1 and 50 falls into this class.
To make the spectral discretization easier, we make the change-of-variables w = 2 R z -1 ∈ (-1, 1). In this way, the PDE in (3.2) becomes:
and we solve for
Conventionally, one spectrally discretizes Ṽk with Chebyshev polynomials for both ρ and w, but the matrices used to represent differentiation and function multiplication in the discretized equations are dense and hard to manipulate. Instead, we set
where T j (w) is the Chebyshev polynomial of the first kind of degree j, X (k) is the matrix of coefficients of Ṽk in the C(3/2) basis, and F (k) is the vector of coefficients of fk in the Chebyshev basis. Although C(3/2) is an uncommon polynomial basis, it is easy and efficient to transform coefficients in the C(3/2) basis to Chebyshev coefficients [18]. Therefore, users of the solver do not need to know about the special ultrashperical polynomial basis.
In this way, (3.5) indicates that we have three cases:
Following [18] on operations related to C(3/2) and [26] on operations related to Chebyshev polynomials, (3.7) can be discretized to the following matrix equations:
where D is a diagonal matrix representing the second derivatives of (1
and H (k) is a matrix with two columns:
One can find a visualization of these operational matrices in Appendix A.
Using the same matrix notation, the discretized matrix equation is:
(3.9)
• If k = 0, then Ũk = (1 -ρ 2 ) Ṽk , and the discretized version of (3.2) is:
We can merge the linear equation with the linear constraint into one matrix equation in all three cases, see for instance [31] for a similar approach. It is desirable to keep the structure and the sparsity of the matrices related to ultraspherical discretization while solving the equation. However, due to the variable s, we do not know the spectrum of the matrices associated with the Chebyshev basis, which is essential if we want to use fast iterative solvers such as an alternating direction implicit (ADI) method [9]. Therefore, we treat all matrices as general dense matrices and solve all matrix equations with the Bartels-Stewart algorithm [8]. Nevertheless, QZ decompositions on sparse penta-diagonal matrices are cheaper and more stable than general dense matrices. Thus, we gain from using ultraspherical polynomials. The solutions X (k) can be transformed to the Chebyshev coefficient matrix for both variables, and then they form the frontal slices of the solution Z to the discretized DFS version of (3.1). Finally, the coefficient matrix that represents the solution of (1.1) can be calculated via
We summarize this solver in Algorithm 1.
Algorithm 1 Fractional PDE solver on the unit disk -Polynomial approximation of z 1/s 1: Input: The coefficient matrix F of f after DFS extension in Chebyshev and Fourier bases. Increase n 2
Solve (3.8), (3.9), and (3.10) for X (k) .
Stack all X (k) in the tube direction to form Z. 8: end while. 9: Calculate W = Z × 2 T 0 (-1) T 1 (-1) . . . . 10: Convert the columns of W into coefficients in Chebyshev basis. Remark 3.1 (Numerical convergence). As a test, consider s = 1/2 so that the extended PDE is a Laplace equation on a cylinder. When f = J 0 (s 01 ρ), where J 0 is the first Bessel function of the first kind, and s 01 is the first nonzero root of J 0 , the solution is
Figure 1 shows the reduction of approximation error obtained by our adaptive solver.
Piecewise Polynomial Approximation of z 1/s . The more challenging case is when z 1/s is not well-approximated by a polynomial. For these values of s, we approximate z 1/s with piecewise polynomials, i.e., where p i is a polynomial for 0 ≤ i ≤ -1, z 0 = 0, and z = R. Then, on the interval (z i , z i+1 ), we have that
with bottom and top boundary conditions
where φi and ψi are implicitly defined. We use these implicit boundary conditions to ensure that Ũi can form an overall continuous solution. On intersection surfaces, solutions on the two sides have identical Dirichlet conditions and opposite Neumann conditions. On each interval, we perform a change of variables to produce similar PDEs as (3.4). These PDEs can be discretized into m matrix equations. We solve on all intervals simultaneously by combining the matrix equations corresponding to the same Fourier mode into one equation. We use the matrix notation from Section 3.1 and get the following three matrix equations:
• If |k| ≥ 2, then we have:
• If |k| = 1, then we have:
• If k = 0, then we have:
In the matrix equations above, we set
where
is the kth Fourier mode solution on (z j , z j+1 ). The first column of H (k) is -d s R fk (ρ)/2 and all other columns are zero. E 1 and E 2 are block diagonal matrices with diagonal blocks S 2,0 B i and K i D 2 , respectively, where B i represents the multiplication of p i and K i represents the multiplication of the Chebyshev basis by (w + z i+1 +z i z i+1 -z i ) 2 . In addition, we set
Then, we can form X 0 from Y (k) and obtain the coefficient matrix of the solution u of (1.1). By construction, the first two columns of D 2 are 0. As a result, E 2 has 2 scattered zero columns. This property is undesirable when using the solver in [31]. Instead, we permute the columns of E 2 so that the first 2 columns are 0, and permute E 1 , Y (k) , B z and H (k) accordingly. After solving the permuted matrix equations, we can easily obtain the original solution. We note that the two solvers are, in fact, equivalent for different values of s. The solver in this subsection can be thought of as a generalized version of the solver in Section 3.1, where the solution only has one piece and no implicit boundary conditions are needed. This generalized solver is described in Algorithm 2.
Remark 3.2. In our numerics, we use Chebfun [16] to automatically partition the extended direction and to approximate the map z → z 1/s with Chebyshev polynomials on each domain with desired accuracy level, i.e., to determine p i and z i in (3.11). Chebfun only divides the domain when a high degree polynomial cannot achieve the resolution, so the piecewise polynomial approximation tends to have as few pieces as possible. From the above discretized systems, one notices that more pieces Algorithm 2 Fractional PDE solver on the unit disk -Piecewise approximation of z 1/s 1: Input: The coefficient matrix F of f after DFS extension in Chebyshev and Fourier basis.
> 10 -10 for any k or j do 5:
Increase n (j) 2 .
Construct the matrices in (3.13), (3.14), and (3.15), and permute E 1 , E 2 , Y (k) , B z and H (k) such that the first 2 columns of E 2 are zero columns while the equations still hold.
Solve for Y (k) .
Stack X (k) 0 to form X 0 . 9: end while. 10: Calculate W = X 0 × 2 T 0 (-1) T 1 (-1) . . . . 11: Convert the columns of W into coefficients in Chebyshev basis. lead to more complicated matrix equations to solve. Consequently, it is beneficial to use Chebfun for approximations. Figure 2 shows the number of polynomial segments needed to approximate the map z → z 1/s on (0, 10) for varying s ∈ (0, 1) with an accuracy of 10 -12 using Chebfun. If we strive for machine precision, the numbers need to be larger. However, we found in practice that 10 -12 gives sufficiently accurate PDE solutions. Comparatively, it is also possible to partition the z-direction into more segments so that z 1/s can be approximated by a low degree polynomial on each interval. For this way of partitioning, we point the readers to [6] for more details.
In this section, we solve (1.1) in the unit square (-1, 1) 2 by spectrally discretizing the truncated, extended PDE (2.2). General rectangles are straightforward to scale to unit squares by a changeof-variables and are thus easy to solve with the same method. It is also worth noting that the solver can be easily generalized to higher-dimensional domains. To this end, we demonstrate the solver application on a cube in Section 5. We first introduce a direct solver that is similar to the solvers in Section 3. The direct solver, can encounter efficiency issues as the equations become large when z 1/s needs to be approximated by a piecewise polynomial with many pieces. In those scenarios, we design a parallelizable solver using domain decomposition in Section 4.2.
4.1. Direct Solver. We first approximate z 1/s as an -piece piecewise polynomial (3.11), with the simplest case being = 1. On each interval (z i , z i+1 ), we have the PDE:
with bottom and top boundary conditions
where φ i and ψ i are implicitly defined to ensure continuity of the entire solution. We employ the change the variables w = 2 z i+1 -z i (z -z i ) -1 on (z i , z i+1 ), and assume the ansatz:
where C(3/2) p (x) is the pth normalized ultraspherical polynomial with coefficient 3/2, T r (w) is the rth Chebyshev polynomial, X (i) is a 3D tensor, and F is a matrix. As before, X (i) and F contain the coefficients of U (i) and f in C(3/2) and Chebyshev basis, respectively. We can then write the discretized problem as a tensor equation by stacking X (i) in the tube direction to form Y:
where A = D -1 M , the first frontal slice of G is -d s RF/2 and the remaining slices are zero, and all other matrices are defined in Section 3. In order to use a tensor analogue of the solver in [31], we first perform a column permutation on E 1 , E 2 and B z , and a frontal slice permutation on Y and G so that the first 2 columns of E 2 are zero and (4.3) still holds.
The
3) can be rewritten as Y × 3 L = H, where L = SB z = I L such that the leftmost 2 × 2 submatrix of L is the identity matrix, and
Then, (4.3) can be combined into a single equation:
where (E 1 ) 1:2 represents the first 2 columns of E 1 . Since the first 2 columns of E 2 are 0 and the leftmost 2 × 2 submatrix of L is the identity, we know that the first 2 frontal slices of Y do not influence the solution. Therefore, we can solve for Y 2 , which contains the rest of the frontal slices, by solving a smaller Sylvester equation:
where R 1 is the first i n
-2 rows of (E 1 ) 2 +1: i n (i)
. After computing Y 2 , it is then straightforward to use the linear constraint to calculate the first 2 frontal slices by
Since we do not know the behavior of the spectrum of R 1 and R 2 , we use a tensor analogue of the Bartels-Stewart algorithm [8] to solve for Y 2 . To be specific, we take real Schur decompositions of the pairs of matrices using the QZ decomposition [8], and then solve for each column fiber of Y 2 . Finally, we recover X (0) as the first n (0) 3 frontal slices of Y, convert it to the Chebyshev coefficient tensor Z and get the coefficient matrix of the solution of (1.1) by Z × 3 T 0 (-1) T 1 (-1) . . . . This direct solver on a unit square is summarized in Algorithm 3. Increase n (i) 3 .
Solve (4.5) for Y 2 and compute Y 1 = H -Y 2 × 3 L.
Stack Y 1 and Y 2 in the tube direction, and form X (0) to be the first n
3 frontal slices of Y. 8: end while. 9: Calculate W = X (0) × 2 T 0 (-1) T 1 (-1) • • • . 10: Convert both columns and rows of W into coefficients in Chebyshev basis.
Remark 4.1 (Numerical convergence). As a numerical example, we consider the case that u = sin(πx) sin(πy) + sin(2πx) sin(2πy), then f = (2π 2 ) s sin(πx) sin(πy) + (8π 2 ) s sin(2πx) sin(2πy) by the spectral definition. Figure 3 (Left) shows the coefficient decay along the extended direction of the discretized tensor solution for different values of s. This plot demonstrates that when the algorithm terminates, the coefficient of the polynomial terms of the discretized solution is small enough. Figure 3 (Right) shows the accuracy improvements of our adaptive algorithm. With the increase of polynomial degree to approximate the extended domain, we obtain a reduction of error between the numerical and analytic solution. 2πy) and different values of s on Ω = (-1, 1) 2 with our spectral solver. Left: the largest coefficient in magnitude on each slice along the extended direction in the discretized solution, i.e., the largest X pqr when r is fixed in (4.2). When s = 2/5 and s = 4/7, z (1/s) is approximated by piecewise polynomials with two pieces. The coefficient patterns show the decay of both pieces. Right: the accuracy achieved with different degrees of freedom in the extended direction. For s = 1/4, we achieve sufficient coefficient decay within five iterations. For s = 4/7, we require two iterations and for s = 2/5, the coefficients of a 28 × 28 × 121 discretization decay below 10 -13 in the first iteration, resulting in a single point on the plot.
Remark 4.2 (Numerical convergence for non-compatible datum). As another example, we consider a non-compatible case with f = 2 and s < 1/2, leading to low regularity of the solution to the fractional PDE [24,Sec. 6.3]. Figure 4 shows the coefficient decay along the extended direction of the discretized tensor solution for different values of s. Again, when the algorithm terminates, the coefficient of the polynomial terms of the discretized solution is small enough, and the number of coefficients to achieve this decay is smaller than that in Remark 4.1.
In the previous section, we solved (1.1) on the unit square by jointly solving (4.1) for all segments. However, the solution of (1.1) only corresponds to X (0) , the solution on the first domain. Therefore, we design a domain decomposition solver, inspired by the hierarchical Poincaré-Steklov method [22], to solve only for X (0) . Suppose we know the Robin boundary conditions for X (0) and Dirichlet boundary conditions for all other X (i) . It is then straightforward to write a tensor equation for each X (i) in the form of (4.4):
where on Ω = (-1, 1) 2 with our spectral solver. The plot shows the largest coefficient in magnitude on each slice along the extended direction in the discretized solution, i.e., the largest X pqr when r is fixed in (4.2). When s = 2/5, z (1/s) is approximated by a piecewise polynomial with two pieces. The coefficient patterns show the decay of both pieces.
conditions of X (i) that are assumed to be known. This means that, once we have G (i) , we are able to solve for X (i) . Equation (4.6) allows us to construct a solution map
such that vec X (i) = S (i) vec G (i) , where vec X (i) reshapes all elements in X (i) to a vector.
In addition, we define the Dirichlet-to-Neumann (DtN) map
The map K (i) converts the Dirichlet boundary conditions G (i) into the Neumann values B (i) on the boundaries. We solve for each column of S (i) by solving a tensor Sylvester equation:
where
3 tensor, and
3 tensor. This suggests that the computation of S (i) j is a parallelizable process. We can calculate the operator K (i) by
where
Our goal is to use the solution operator S (0) to solve for X (0) . Thus, we must obtain the implicit Dirichlet boundary condition. Next, we show that we can carefully merge the solution maps S (i) and DtN maps K (i) . As a result, we obtain solution maps and DtN maps that work on several domains simultaneously, and these merged maps can help us find the desired boundary condition.
The boundaries for each domain consist of one upper surface and one lower surface so that we can separate the boundary conditions into two parts and the DtN operator into four parts, i.e., vec
Then we can follow [22] to merge the solution and DtN operators:
In other words, given the Dirichlet boundary conditions on the lower surface of the ith domain and the upper surface of the (i + 1)st domain, K (i,i+1) enables one to compute the Neumann conditions on those surfaces, and S (i,i+1) allows one to calculate the Dirichlet condition of the overlapping surface.
We need G (0) u or G
(1) v to compute X (0) . To achieve this, we can merge the operators for X (1) , . . . , X ( -1) to get S (1, -1) and K (1, -1) , and then merge them with S (0) and K (0) in the final step. In particular, there are two ways of merging the operators:
(a) Starting from the top piece, we form S ( -1) and K ( -1) . We then iterate downwards from i = -1 to i = 0, form new operators for each piece, and merge them with the operators from the previous iteration. This is merging in a sequential way. (b) Since the maps for each piece are independent, we can form the operators on all domains in parallel. Then, we merge in a hierarchical manner, merging two of them simultaneously, and these merging operations are parallelizable.
In summary, this domain decomposition solver is Algorithm 4.
Remark 4.3. We can also use domain decomposition to construct a parallel solver for the disk domain. The solution operator S (i) ∈ C n 1 n (i) 2 m×2n 1 m on the ith domain takes in the Dirichlet boundary coefficients on the top and the bottom surfaces, and returns the coefficients of Ũ (i) . Although we cannot construct S (i) in one setting due to partial regularity, we discover that row (k -1)n 1 n
of S (i) corresponds to the (k -1 -m/2)th Fourier mode of the solution. In this way, we can construct these rows with the linear system converted from the generalized Sylvester equation. It is then straightforward to construct the DtN map K (i) from S (i) , and we can use the hierarchical Poincaré-Steklov method described for the rectangle domain. Increase n (i) 3 .
Form solution and DtN maps on each domain.
Merge the maps to get S (1, -1) and K (1, -1) .
Merge with S (0) and K (0) to get solution tensor X (0) on the first domain. 9: end while. 10: Calculate W = X (0) × 2 T 0 (-1) T 1 (-1) • • • . 11: Convert both columns and rows of W into coefficients in Chebyshev basis.
In this section, we present two examples. In the first example, we extend our 2D solver to 3D to solve the fractional elliptic PDE on Ω = (0, 1) 3 . The second example is an optimal control problem with a fractional PDE constraint.
Fractional PDE on the Cube. We can extend our solver from Section 4 to solve fractional PDEs on the unit cube. Using ultraspherical polynomials for the cube dimensions and Chebyshev polynomials for the extended direction, the discretized problem is the following:
where Y is formed by stacking 4D tensor solutions of the discretized problem on each interval along the fourth dimension, F is a 3D tensor representing the polynomial coefficients of the initial condition, and all other matrices have been defined in Section 4. We solve (5.1) by merging the linear constraint into the tensor equation and using a 4D Bartels-Stewart algorithm analogue to obtain the solution directly. For a numerical example, we consider the simple case that f = (3π 2 ) s sin(πx) sin(πy) sin(πz) + (12π 2 ) s sin(2πx) sin(2πy) sin(2πz) so that the analytic solution is u = sin(πx) sin(πy) sin(πz) + sin(2πx) sin(2πy) sin(2πz). Figure 5 (Left) shows the coefficient decay along the extended direction of the discretized tensor solution for different values of s. This plot shows that the coefficient of the discretized solution is small enough when the algorithm finishes. Figure 5 (Right) shows the accuracy improvements in our adaptive algorithm when we increase the degrees of freedom by allowing higher polynomial degrees to approximate the solution in the extended direction. Similar to the problem on the square, polynomial coefficients of the discretized solution decay exponentially along the extended direction, which allows us to find an accurate numerical solution with only a few degrees of freedom. 5.2. Optimal Control Problem. We consider the optimization problem given by min
where u d is a given function and α is the control penalty parameter. We solve this problem via a direct solver. Specifically, we express the optimality condition as two fractional PDEs, discretize both, combine them into one tensor equation and obtain the best u and z directly. In particular, we solve
where we have eliminated the so-called gradient equation. For simplicity, we consider solving (5.3) directly on the unit square. Let U and P be the coefficient tensors for the extensions of u and p, respectively. Then, U and P satisfy the following tensor equations:
where the first frontal slice of
(-1) 0 . . . 0 -U d , and U d is the coefficient matrix of u d . To distinguish between matrices used in the equations for u and p, we use superscripts (u) and (p). The matrices E 1 , E 2 and B y are defined in Section 4. We can then rearrange (5.4) so that we only have one tensor Sylvester equation and a linear constraint:
where Y is formed by stacking U and P along the tube direction, G is a zero tensor except for the last frontal slice of y . We can then solve for Y, get U and P, and calculate the coefficient tensor Q of q by Q = -1 α P. For numerical demonstration, we take Ω = (-1, 1) 2 , α = 10 -2 , and u d = (1+α(2π 2 ) 2s ) sin(πx) sin(πy). The analytic solution of (5.2) with this data is u = sin(πx) sin(πy) and q = (2π 2 ) s sin(πx) sin(πy). Figure 6 shows the performance of our adaptive direct solver. As demonstrated, our method generates accurate numerical solutions for both u and q with only a few degrees of freedom in the extended direction. ) sin(πx) sin(πy) and Ω = (-1, 1) 2 . As the degrees of freedom in the extended direction increase, the approximation of both u and q improves. When s = 1/4, we perform four iterations, which add more degrees of freedom along the extended direction, before our algorithm terminates. When s = 2/5, our first trial guarantees enough decay in the coefficients for both u and q, resulting in a single point for each on the plot.
In this paper, we present a spectral method that uses ultraspherical and Fourier polynomials to solve fractional Laplacian equations on square and disk domains via the Caffarelli-Silvestre extension. Based on the value of the fractional exponent s, we decompose the PDE along the extended domain. We show a direct method that finds solutions on all sub-domains through one tensor equation, and a parallelizable domain decomposition solver generated from the hierarchical Poincaré-Steklov method. Numerical tests suggest that coefficients of the solutions decay exponentially along the extended direction, and we can recover accurate discretized solutions with a few degrees of freedom. Our method is easily generalized to problems of higher dimensions, such as solving fractional PDEs on cubes, and it can be used to accurately compute solutions of optimal control problems. For future work, we will develop spectral solvers for fractional operators with variable coefficients (1.2) through the extension scheme, yielding an approach that can be used to solve more general PDEs and optimal control problems as in [27].
In this appendix, we provide the readers with some details of the matrices we use during discretization in Sections 3 and 4.
• D is a diagonal matrix representing second derivative of (1 -ρ 2 ) C(3/2) (ρ), with diagonal elements D j,j = -(j(j + 3) + 2). • M is a symmetric penta-diagonal matrix with 0 super and sub diagonals representing multiplication of 1 -ρ 2 in C(3/2) basis, with elements M j,j = 2(j + 1)(j + 2) (2j + 1)(2j + 5) , M j,j+1 = 0, M j,j+2 = -1 (2j + 3)(2j + 5) (j + 4)!(2j + 3) j!(2j + 7) .
• B 1 and B 2 are two banded matrices representing multiplications in different ultraspherical polynomial basis, where the bandwidth is determined by the degree of the polynomial approximation in the respective basis. In addition, B 1 can be shown to be a Toeplitz-plus-Hankel-plus-rank-1 operator [26], and both B 1 and B 2 satisfy three-term recurrence [30,Chpt.6]. • D 2 represents second derivative of Chebyshev basis, with the form
Here, elements of D 2 on the diagonal and first upper top diagonal are all 0.
• S 2,0 converts Chebyshev basis to C (2) basis, and can be calculated by S 2,0 = S 1 S 0 , where