{'fraction_non_alphanumeric': 0.0586890243902439, 'fraction_numerical': 0.02269854253420583, 'mean_word_length': 4.092587333143993, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The numerical solution of fractional integral equations via orthogonal polynomials in fractional powers

We present a spectral method for one-sided linear fractional integral equations on a closed interval that achieves exponentially fast convergence for a variety of equations, including ones with irrational order, multiple fractional orders, non-trivial variable coefficients, and initial-boundary conditions. The method uses an orthogonal basis that we refer to as Jacobi fractional polynomials, which are obtained from an appropriate change of variable in weighted classical Jacobi polynomials. New algorithms for building the matrices used to represent fractional integration operators are presented and compared. Even though these algorithms are unstable and require the use of high-precision computations, the spectral method nonetheless yields well-conditioned linear systems and is therefore stable and efficient. For time-fractional heat and wave equations, we show that our method (which is not sparse but uses an orthogonal basis) outperforms a sparse spectral method (which uses a basis that is not orthogonal) due to its superior stability.


Introduction
Fractional differential equations (FDEs) and fractional integral equations (FIEs) have received a great deal of attention in the literature, not only because they appear in many scientific fields 1 , but also because their solutions are difficult to compute with traditional approaches. Numerous standard numerical methods 2 have been applied to FDEs and FIEs, however they only achieve algebraic convergence. This is because the fractional integral operators can introduce algebraic singularities that are not well approximated by polynomials.
The fractional integral of order µ > 0 is a singular integral defined as with a < x < b and where I µ a + and I µ b − are referred to as, respectively, the left-sided and right-sided fractional integrals. Consider the left-sided fractional integral of the unit constant function, I µ −1 + [1](x), which is (1 + x) µ /Γ(1 + µ) (see (6)). Hence, I µ −1 + [1](x) has an algebraic singularity at x = −1 as illustrated by Fig. 1a and the Chebyshev polynomial approximation of I µ −1 + [1](x) on [−1, 1] converges at the slow algebraic rate of O(n −2µ ) as indicated by Fig. 1b, where n is the polynomial degree. By contrast, with the Jacobi fractional polynomial (JFP) basis that we shall introduce, algebraic singularities are incorporated into the basis, which allow us to compute fractional integrals of weighted polynomials in fractional powers exactly. Fractional derivatives have different definitions, among which the Riemann-Liouville and Caputo types are commonly used. We denote them by, respectively, D ν RL and D ν C , where ν > 0 and they are defined as compositions of fractional integral and standard derivative operators (but in different orders) as follows [33, (2.22) Here D m denotes the m-th order (standard) differential operator: D m ϕ(x) = d m dx m ϕ(x), and I µ can be either a left-sided or right-sided fractional integral. In this paper we focus on problems involving fractional integral operators. We do consider FDEs, however we reformulate them as FIEs, similar to the integral reformulation method discussed in [40,Section 5] in which ODEs are reformulated as ordinary integral equations. In future work we shall construct differentiation matrices for JFP bases, which will make our method, which we refer to as the JFP method, directly applicable to FDEs, without the need for integral reformulation.
We now give a brief overview of existing methods that can achieve exponential convergence for FDEs and FIEs and compare them to the method we shall present. All of these approaches incorporate algebraic singularities into the basis functions.
The fractional collocation methods in [31,48] and the fractional Galerkin method in [36] achieve spectral convergence when the solution is smooth. However, in addition to giving rise to dense linear systems, these methods revert back to algebraic convergence if the solution has a singularity at the left endpoint, which is expected in practice (e.g., Example 2). We shall overcome this problem in the JFP method by including an extra singularity into the basis such that it captures the singularity of the solution. This idea is also applicable to the methods in [31,48,36] and should restore spectral convergence to the examples where they achieved algebraic convergence.
Our work is strongly influenced by the method proposed by Hale and Olver [25], in which direct sums of appropriately weighted Jacobi polynomials are used as basis functions. We shall refer to this method as the sum space method and we note that the basis functions are related to the "generalized Jacobi functions" of [9] and the "polyfractonomials" of [48]. Similar to the ultraspherical spectral method [38] for ODEs, in the sum space method, the domain and range of operators in FDEs and FIEs are represented in different bases to ensure that the matrix representations of operators are banded. This enables a fast algorithm with linear complexity and exponential convergence for a wide range of problems. However, we shall find in Section 6 that the sum space method leads to tremendously large expansion coefficients (larger than 10 100 ) for the solution of an important family of FIEs that arises in the solution to the time-fractional heat/wave equation (D µ t u = D 2 x u). These large sum space coefficients require very high precision for the accurate computation of solutions and is thus very expensive. By contrast, this family of FIEs poses no difficulties to the JFP method whose solution coefficients are bounded below 1. We shall clarify this difference in the performance of the two methods by using bounds on Jacobi expansion coefficients for analytic functions to bound the JFP coefficients and by illustrating that the largest expansion coefficient in the sum space basis grows at the same (super-exponential) rate as the largest coefficient of the power series expansion of the solution.
We consider our method to be the successor of that of Bhrawy and Zaky [7]. They applied a change of variables to classical Jacobi polynomials such that the algebraic singularities of the resulting basis, the JFP basis (which is called thus for reasons we explain in Section 3), conform to those of the solution 3 . The JFP basis inherits many desirable properties of classical Jacobi polynomials, including orthogonality, see [7]. Therefore, it can potentially be used to solve FIEs and FDEs just as classical orthogonal polynomials are used to solve ODEs. We believe this approach has not been sufficiently analyzed nor developed to its full potential. In this paper, we develop new algorithms for computing the entries of matrices that represent the action of fractional integration operators on the JFP basis. We shall refer to these matrices as fractional integration matrices. We also illustrate that the algorithms for fractional integration matrices (including the one used by Bhrawy and Zaky) are unstable but can be "pseudo-stabilized" by integrating high-precision computing [5] automatically such that the algorithms output accurate results. We also investigate empirically the performance of the parameter-dependent JFP method for a wider range of parameters than those considered in [7].
Our pseudo-stabilization technique discussed in Appendix A is essential to the scalability of the JFP method and thus also for its application to practical computational problems. We emphasize that highprecision computations are required only for the computation of fractional integration matrices and not for the solution of the resulting linear systems. In addition, we do not require any of the inputs to the FIE or FDE (e.g., variable coefficients, the function on the right-hand side of the equation, boundary conditions, etc.) to be computable to high-precision accuracy.
The following is an outline of the paper: We introduce the basic constituents of the JFP method in Sections 2 and 3 (matrix representations of operators on quasimatrices, high-precision floating-point numbers, the JFP basis, etc.). We then focus on the properties (Section 4) and computation (Section 5) of fractional integration operators and matrices acting on the JFP basis. The JFP method is used in Section 6 to solve a variety of FIEs 4 , including an FDE and fractional PDE reformulated as FIEs, and comparisons are made with the sum space method. We conclude the paper with a summary and a discussion of topics for future work. Appendix A is devoted to the above-mentioned pseudo-stabilization of an algorithm for computing fractional integration matrices.

Preliminaries
Henceforth, we consider left-handed fractional integrals and derivatives on compact intervals, which we can take to be defined on [−1, 1], without loss of generality. To simplify the notation, we let and let I denote I 1 , the standard integral operator. To avoid ambiguity, we shall denote (infinite) identity matrices by 1 . We shall make use of the fact that fractional integral operators form a semigroup [33, (2.21)]:

Quasimatrix notation and fractional monomial bases
As in [25,40], we shall adopt quasimatrix notation because of its convenience. A quasimatrix can be thought of as a matrix with a function in every column rather than a vector. For instance, we let where x ∈ [−1, 1]; M0(x) (or simply M0) is an example of a quasimatrix that has the (shifted) monomial basis functions in its columns, ordered by degrees, hence it has a countable infinity of columns and we write M0 ∈ R [−1,1]×N 0 (similar to the way in which we might state A ∈ R m×n for a matrix A). Let γ ∈ R, then we define The columns of M γ 0 can be viewed as elements of a fractional monomial basis if γ has a positive non-integer value. We shall also make use of a normalized version of M γ 0 , obtained by replacing We denote the quasimatrix M γ 0 weighted by (1 + x) δ , i.e., We can express the action of the fractional integral operator on M γ δ using the fact that [33, (2.44)]: where β > −1, to ensure the integral exists. Throughout, we shall use as a dummy variable. The following result is an immediate corollary of (6): Proposition 1. Let µ, γ > 0 and µ = kγ, where k is a positive integer (i.e., k ∈ N+), then I µ maps M γ δ to itself via an infinite matrix whose k-th subdiagonal is nonzero, i.e., We only consider the case µ = kγ because, as we shall find in Section 6, this is a necessary condition for the solutions to the FIEs to have an expansion in the JFP basis.
We say that a matrix A has bandwidths (λ, ω) (or bandwidths(A) = (λ, ω)) if the entries of A satisfy Ai,j = 0 for i−j > λ and j −i > ω. For example, if A is upper triangular, then we say that bandwidths(A) = (0, ∞) if A is infinite and bandwidths(A) = (0, N ) if A is N × N . As another example, regardless of whether the matrix (7) is infinite or N × N , it has bandwidths (k, −k) .
For quasimatrices P and Q and a weight function w, we let P, Q w denote a matrix of L 2 (w) inner products between the functions forming the columns of P and Q. As an example we shall use later, let P = Q = M0 and let w be the Jacobi weight function, w = w α,β , where Also, let and define e k , k ≥ 0 to be the basis vector then M0e k is the function in column k of M0 (i.e., (1 + x) k ) and the (i, j)-entry (with i, j ≥ 0) of the infinite matrix M0, M0 w α,β , which we denote by M (α,β) , is where B is the beta function. Since M0 ∈ R [−1,1]×N 0 , in (10) we interpret M 0 M0 ∈ R N 0 ×N 0 as the outer product of two vectors 5 whose entries are functions on [−1, 1]. In the general case, the (i, j) entry of P, Q w is defined as Pei, Qej w (cf. (11)). If P = Q, then the (infinite and symmetric) matrix P, P w is known as a Gram matrix [30, p. 441], hence M (α,β) is a Gram matrix.

Jacobi polynomials
We let the columns of the quasimatrix P (α,β) consist of Jacobi polynomials, ordered by degrees, where the P (α,β) n are orthogonal with respect to the inner product (9) with 6 α, β > −1. The P (α,β) n are related to the (shifted) monomial basis functions according to 5 Not to be confused with the notation in [40], where P wQ is used to denote P, Q w . 6 In [43, 4.22] and [3,9] Jacobi polynomials are considered for α ≤ −1 or β ≤ −1.

5
which follows from the identities [13, (18.5.7)] and P (β,α) n (x) = (−1) n P (β,α) n (−x) (see [13, where C (α,β) is an upper triangular matrix with entries We can use the orthogonality of the Jacobi polynomials to obtain the inverse of C (α,β) . Define P (α,β) 2 = P (α,β) , P (α,β) w α,β , then P (α,β) 2 is a diagonal matrix with the squared norms of the Jacobi polynomials on the main diagonal, which are given in [13, (10) and (12), it follows that and thus Since C (α,β) is upper triangular, so is C (α,β) −1 and since the entries of M (α,β) are known (see (10)), can be computed as a finite sum of i + 1 terms 7 . We shall use the matrices acting on Jacobi polynomials given in Table 1 as building blocks to construct matrix representations of the operators we need. For example, the integration matrix (which is a representation of the standard integration operator) is diagonal and follows from the weighted differentiation matrix in Table 1: We shall also use the following matrix representation of the multiplication operator (which we shall refer to as a multiplication matrix): Proposition 2. We have that where the multiplication matrix L (α,β) has bandwidths (1, 1).
Proof. By the actions of the matrices in Table 1, has bandwidths (1, 1) since it is a product of matrices with bandwidths (1, 0) and (0, 1).

High-precision computation
We shall use high-precision computation [5] in the Julia programming language 8 due to the instability of our algorithms for computing fractional integration matrices, which are discussed in Section 5. As in [17], "q-bit precision" means the machine epsilon (or relative accuracy of the floating-point numbers we compute with) is 2 1−q . For example, the widely-used double-precision arithmetic has 53-bit precision and a machine epsilon of 2.22e-16. Some of the commonly used precisions are given in Table 2.

Jacobi fractional polynomials
The JFP basis functions that we shall use are weighted polynomials in fractional powers defined as follows, where the choices of b ∈ R and p > 0 will be problem (and algorithm)-dependent, as we shall find in subsequent sections. The corresponding quasimatrix is denoted by Q (α,β,b,p) . Note that the relation between x and y defines a bijection from [−1, 1] to [−1, 1], see Fig. 2a for an example. The first few Legendre fractional polynomials 9 with b = 0, p = 2 are illustrated in Figs. 2b and 2c. We use the term Jacobi fractional polynomials, or JFPs, for our basis functions (17) because they define polynomials in fractional powers (or in the mapped variable y) and to distinguish them from fractional Jacobi polynomials [22], which are defined by fractional differential equations. The JFPs are related to the basis functions used by Bhrawy and Zaky by an affine transformation, hence we refer to [7] for a discussion of some properties of JFPs.
The following result will allow us to relate the JFP basis to a weighted fractional monomial basis in x. where Proof. It follows from the definition (5) and the relation between x and y in (17) that the right-hand side of (18) simplifies to Corollary 4. The JFP basis is related to a fractional monomial basis as follows Proof. It follows from (12) and Lemma 3 that

4 Properties of fractional integrals applied to the JFP basis
The results derived in this section will motivate the two algorithms for computing fractional integration matrices that are discussed in the next section. The following result gives the action of I µ on the JFP basis and is the foundation of the first algorithm to be presented in the following section. The subsequent results in this section will inform the second algorithm.
Proof. We derive the result from (20) (where D b,p is defined in (19)), Proposition 1 with δ = b/p, γ = 1/p and the relation between x and y given in (17): where we used the fact that (7)) and hence the result (21) follows. The lower banded structure of (21) is a consequence of the fact that C (α,β) and C (α,β) −1 are upper triangular and the structure of (7).
We know from [24, Lemma 2.2] that fractional integrals satisfy recurrence relations in ultraspherical bases. We now derive a general recurrence relation satisfied by fractional integrals, which we shall apply to Jacobi fractional bases.
Proof. By definition of the fractional integral, In order to apply Proposition 6 to the JFP basis, we need to construct matrix representations (which are given in the following lemmata) of multiplication and integration operators.

Lemma 7 (Multiplication matrix). Suppose p is a positive integer and let
where 1 denotes the identity matrix, and X (α,β) b,p has bandwidths (p, p).
, and denote the matrix representing the integral operator applied to the JFP basis by I and has bandwidths (p, p).
Proof. Using (16), the matrices in Table 1, the definition of the JFP basis in (17) and making the change of variables 1+t 2 = 1+s 2 p , it follows that hence we reach (22).
According to Table 1, the bandwidths of L has bandwidths (p, p).
, then the fractional integration matrix defined in Theorem 5 satisfies the following banded Sylvester equations where X Proof. The bandwidths of the matrices follow from Lemmas 7 and 8. The first equation (23) follows immediately from the commutativity of fractional (and integer-order) integration matrices stated in (3).

Algorithms for computing fractional integration matrices
In this section we discuss the following algorithms for computing the entries of fractional integral matrices: b,p,µ by solving either of the banded Sylvester equations (23) or (24). This requires the pre-computation of the first p columns of I (α,β) b,p,µ , which can be obtained from Algorithm 1.

Algorithm 1
In Algorithm 1, each column of I (α,β) b,p,µ is computed independently (hence making Algorithm 1 amenable to parallelisation) by solving an upper triangular linear system (recall that C (α,β) is upper triangular). If I (α,β) b,p,µ has bandwidths (k, ∞), k ≥ 0, then obtaining the j-th column (with j ≥ 0) requires the solution of a (j + 1 + k) × (j + 1 + k) system. Hence, computing N columns of I (α,β) b,p,µ has a complexity of O(N 3 ). However, the condition numbers of N × N sections of C (α,β) (which we denote by C (α,β) 0:N −1,0:N −1 ) grow exponentially [21] with N and thus high-precision arithmetic is required to maintain accuracy. In q-bit precision, each multiplication costs O(q log q log log q) operations with the Schönhage-Strassen algorithm 11 , hence Algorithm 1 has a complexity of O(qN 3 log q log log q). As shown in Appendix A, q = O(N ) is required to maintain accuracy due to the exponentially growing condition numbers of C (α,β) , there a pseudo-stabilized version of Algorithm 1 has a complexity of O(N 4 log N log log N ). For Algorithm 1, we may choose any p > 0 (rational or irrational), such that µp = k ∈ N+.

Algorithm 2
The Sylvester equations (23) and (24)  where an entry is considered to be zero if an index is negative, and thus Am,n+pBn+p,n = Cm,m−p:m+pAm−p:m+p,n − Am,n−p:n+p−1Bn−p:n+p−1,n.
This shows that if the first p columns of A are known (which, for A = I (α,β) b,p,µ , can be obtained from Algorithm 1), then the subsequent columns can be obtained from the recurrence relation defined by (28).
We have found that Algorithms 1 and 2 are roughly equally unstable. In Algorithm 2, the numerical errors propagate due to column-by-column recurrence, and thus the accuracy is limited by the precision we begin with. By contrast, in Algorithm 1, each column can be computed independently to any desired precision. Another disadvantage of Algorithm 2 is that we can only use it for rational µ since we require µp = k ∈ N+ and p must be a positive integer (see Theorem 9). However, the complexity of Algorithm 2 is lower than that of Algorithm 1 by one order. Specifically, the first N columns of I (α,β) b,p,µ involve the computation of O(N 2 ) entries, each of which requiring up to 4p + 3 multiplications, leading to a total cost of O(pN 2 q log q log log q) in q-bit precision. As shown in Appendix A, q = O(N ) and we conclude that Algorithm 2 has an overall complexity of O(pN 3 log N log log N ). Fig. 4 illustrates the growth of errors for a high-precision computation of I   (24) and (23) in 256-bit precision. The difference in accuracy between the two solutions is negligible.
In Appendix A we describe how we pseudostabilize Algorithm 2 by automatically incorporating highprecision arithmetic in such a manner that the first N columns of the matrix I (α,β) b,p,µ are computed to a prescribed accuracy while minimizing computational cost.
Algorithm 2 can be stabilised with Algorithm 1 at the expense of higher computational complexity as follows. We can use Algorithm 1 to compute the first p columns of I (α,β) b,p,µ and also its N -th column for some N > p with high-precision arithmetic, which costs O(N 3 log N log log N ) operations. Then the entries in the intervening columns (columns p + 1 to N − 1) can be assembled into a single vector and approximated as the least squares solution to an overdetermined system. We have found that this system is well-conditioned and can thus be computed accurately in a lower precision (for example, we have achieved roughly 10 −13 accuracy in double precision). However, we found that this procedure is slower than the pseudo-stabilized Algorithm 2.
Another strategy to stabilize Algorithm 2, which we shall pursue in future work, is to derive asymptotic approximations to the entries of I (α,β) b,p,µ in column N as N → ∞, which could be used (as described above) to set up an overdetermined system and stably compute the entries of I (α,β) b,p,µ . This approach would have the optimal O(N 2 ) complexity.

Spectral methods using Jacobi fractional polynomials
We first test our methods by applying them to problems that have solutions expressible in terms of Mittag-Leffler functions in Section 6.1 before tackling more challenging problems in Section 6.2. In Example 3 of Section 6.1 we shall consider a striking example in which the performance of the JFP method and the sum space method of [25] are radically different.

FIEs, FDEs and a fractional PDE with Mittag-Leffler function solutions
The general Mittag-Leffler functions are defined as and they are the general solutions to the following fractional integral equations (which can be verified using (6)): Theorem 10 ([27, Theorem 13.4]). If Re ν, Re µ > 0, then the solution to the fractional integral equation is given by 12 For the problems we consider, α and β are positive real numbers, in which case E α,β (z) is an entire function. Various methods for computing Mittag-Leffler functions have been implemented, see [19,23,29,41]. We shall use [19] to benchmark the accuracy of our methods.
We can express the solution (31) in the normalised monomial basis (4) as u( If we define u α,β,b,p to be the coefficients of the solution in the JFP basis, i.e., u(x) = Q (α,β,b,p) (x)u α,β,b,p , then using (20), it follows that where For I µ u to be defined in the JFP basis, a necessary condition is b > −p, to ensure that the basis Q (α,β,b,p) is integrable (because it contains monomials of the form (1 + x) b/p ). Equations (33) and (30) therefore hold only if there are integers k * ≥ 1 and n ≥ 0 such that If we let u( ) denote the -th entry of the vector u for ≥ 0, then u b,p ( ) is the coefficient of the monomial basis function (1 + x) (b+ )/p . Likewise, uν,µ( ) is the coefficient of (1 + x) ν−1+ µ . Hence, (33) and (35) imply that u b,p (n + mk * ) = uν,µ(m), m ≥ 0, and otherwise, u b,p ( ) = 0, = n + mk * . Specifying u b,p thus, we can obtain u α,β,b,p by solving (34), which is an ill-conditioned system whose condition numbers grow exponentially and relies on knowing the exact coefficients of the solution in the monomial basis. Instead, we shall obtain the JFP coefficients u α,β,b,p by solving the integral equation (30) in the JFP basis which, as we shall see, result in systems whose condition numbers are bounded by a quantity that grows linearly in λ.
If the conditions (35) are satisfied and we set u(x) = Q (α,β,b,p) (x)u α,β,b,p , then (30) becomes 12 In [27,Theorem 13.4] there is a factor of Γ(ν) in the formula for u(x), which we believe is due to a miscalculation.
where (1 + x) ν−1 = Q (α,β,b,p) (x)f , i.e., f contains the coefficients of (1 + x) ν−1 in the JFP basis. Using (17) and the fact that (b + n)/p = ν − 1, it follows that (1 + x) ν−1 = 2 (b+n)(1/p−1) (1 + y) b+n , hence f has n + 1 nonzero entries. In the simplest case, n = 0 and f = f0e0, where f (0) = f0. By Theorem 5, the matrix I (α,β) b,p,µ has bandwidths (k * , ∞) and we obtain u α,β,b,p by truncating and solving the system If (35) is satisfied and we use the relation between x and y in (17), then the exact solution (31) can be expressed as which we represent in the JFP basis as Comparing these two equations, we conclude that by truncating and solving (36), we are approximating an entire function of y by a truncated expansion in Jacobi polynomials. Therefore, we can use bounds on Jacobi coefficients of Jacobi expansions of analytic functions, which are given in [49], to obtain bounds 13 on the entries of u α,β,b,p . For example, if we use Chebyshev polynomials, then holds for any ρ > 1 and thus the coefficients decay super-exponentially [46,Chapter 8]. Bounds similar to (38) hold for any Jacobi parameters α, β > −1, see [49].
Remark 1. In (38), ρ refers to the sum of the semi-axes of a Bernstein ellipse Eρ (in the complex y-plane) on and within which the function we are approximating (the Mittag-Leffler function) is analytic and M is the maximum modulus of the function on Eρ.
We shall find it instructive to compare the bounds on the JFP coefficients in (38) to an estimate of the largest coefficient in the (fractional) normalized monomial basis, i.e., max n≥0 | uν,µ(n)|, where (see (32)), The large-λ regime is of particular interest since it will arise in the solution to a fractional heat/wave equation. We let 0 < µ < 1 and consider the case λ 1. As n increases from 0, the magnitude of uν,µ(n) increases and reaches a maximum at which | uν,µ(n + 1)| ≈ | uν,µ(n)|, ⇒ Γ(ν + (n + 1)µ) then, as n increases further, | uν,µ(n)| decays super-exponentially as n → ∞. If λ is large, then (40) is satisfied for large n and we can use a leading order approximation of the ratio of gamma functions in (40) from [47,Eq. (1)] as n → ∞, Using the approximation (nµ) µ ≈ 2 µ λ, Γ(ν + nµ) ∼ (nµ) ν Γ(nµ), n → ∞ (which again follows from the result in [47]) and Stirling's approximation in (39), it follows that max n≥0 | uν,µ(n)| ≈ Since 0 < µ < 1, this shows that the largest monomial coefficient grows super-exponentially in λ 1. We shall find evidence of super-exponential growth (as a function of λ) in the maximum coefficient of the sum space solution in Example 3, which renders the method catastrophically unstable for large λ. By contrast, we shall find in this same example that the JFP coefficients (with bounds given in (38)) are bounded below 1 for all the values of λ that we consider. Example 1. We first consider a simple instance of (30), viz., in Theorem 9 are satisfied, which enables us to use Algorithm 2 (see Section 5). Fig. 5 illustrates that truncating and solving (36) converges super-exponentially for a wide range of fractional orders (including irrational orders), as expected. We let α = β = b = 0, vary µ and p and fix the value of k * = 1 (i.e., µp = 1) in Figs. 5a and 5b and in Figs. 5c and 5d, we let k * ≥ 1. Figs. 5c and 5d suggest that the choice k * = 1 yields faster convergence than choosing k * > 1. However, for some equations (e.g., (44)), larger k * can have (sometimes significantly) better performance. We shall pursue the question of determining optimal choices of parameter values (which might involve optimising the bounds (38)) in future work.
Example 2. We modify the previous example by letting the right-hand side have a singularity,  Fig. 6b.
Example 3. We consider the problem with the exact solution u(x) = E 1/2,1 (−λ 2 √ 1 + x) (see Theorem 10), which can also be expressed as [27] u = exp(λ 4 (1 + x))erfc( where erfc denotes the complementary error function. We shall compare the effect of λ on the rate of convergence of the sum space method of [25] and the JFP method. We let the constant λ grow quadratically in (44) since this is also the case in the time-fractional fractional heat/wave equation that we shall consider in Example 4. In the sum space method for (44), the solution is expanded as a sum of Legendre and weighted secondkind Chebyshev polynomials, or, expressed in terms of Jacobi polynomials, µ = e/π, p = π/e µ = π/2, p = 2/π µ = 1/π 2 , p = π 2 b p = 1/µ for different irrational µ where S is the 'interleaved' sum space basis, i.e., and c = (a0, b0, a1, b1, · · · ) are the interleaved coefficients 14 . As shown in [25], the half-order integral operator I 1/2 maps S to itself via a tridiagonal matrix, say I 1/2 . Therefore, in the sum space basis, (44) becomes S 1 + λ 2 I 1/2 c = Se0, and the coefficients c can be obtained by solving the tridiagonal system By contrast, with the JFP method, the system (36) (with λ → λ 2 ) has bandwidths (k * , N ), where N is the size of the truncated system and k * = µp = p/2. Fig. 7 gives comparisons of the numerical results for solving (44) with the sum space method and the JFP method (with α = β = 0, b = 0 and p = 2 in (36)). 14 We assume the coefficients c decay sufficiently fast so that the implicit change in the order of summation in (46)  When using double precision in the sum space method, Fig. 8a shows that for λ = 2, 3, numerical errors in the solution on the order of 10 −2 (and larger) are visible. This is because some coefficients are as large as 10 15 (see Fig. 7a), which can lead to large cancellation errors. Using high precision, see Fig. 7b, shows how astronomically large the sum space coefficients become as λ increases, which requires higher precision to compute and larger truncation sizes to ensure the solution achieves a specified accuracy. More precisely, solving an N × N truncated version of the tridiagonal system (48) in q-bit precision has a complexity of O(N q log q log log q). Computing an accurate solution to (44) with the sum space method requires q-bit precision, with q = O(λ 4 ), and the required truncation size N also grows as O(λ 4 ), hence the overall complexity of the sum space method is O(λ 8 log λ log log λ).
In the JFP method for (44), O(N 2 q log q log log q) operations are expended in computing the first N columns of the matrix I (0,0) 0,2,1/2 in q-bit precision using the pseudo-stabilized Algorithm 2 (see Appendix A). Solving the resulting N × N truncated version of the system (36) for (44) in q-bit precision costs O(N 2 q log q log log q) operations. As shown in Appendix A, q = O(N ) and for (44), N = O(λ) and q is independent of λ with 15 q < q, hence the JFP method has an overall complexity of O(λ 3 log λ log log λ), which is several orders lower than that of the sum space method (O(λ 8 log λ log log λ) for (44)). To reiterate, notwithstanding the fact that the JFP method requires the solution of a lower banded system (which, in addition requires high precision to compute), whereas the sum space method only requires the solution of an explicitly known tridiagonal system, the JFP method still has much lower complexity than the sum space method for (44).
The numerical results in Fig. 7b indicate that the magnitude of the largest sum space coefficient is similar to that of the largest power series coefficient in the normalized monomial basis, i.e., on the order of O exp(2λ 4 ) (see (41) but with λ replaced by λ 2 and µ = 1/2). The normalized monomial and sum space coefficients of the solution to (44) can be related using the following: we have that k,n , see (4), (12) and (13). Then since u = E 1/2,1 (−λ 2 √ 1 + x), and using the sum   Table 2). The largest coefficient in b for each λ is approximately exp(2λ 4 ), which, for λ = 2, . . . , 5, is on the order of 10 14 , 10 70 , 10 222 and 10 543 , respectively. The JFP coefficients in c are all bounded below 1, even for much larger values of λ. d shows the error of the JFP solution to (44), evaluated in double precision. For the JFP method, the truncation size (or number of coefficients) required to achieve 10 −14 accuracy grows linearly in λ, whereas for the sum space method, the truncation size required for a fixed accuracy grows as O(λ 4 ). For the sum space method e, the (2-norm) condition numbers of (48) grow roughly as exp(O(n 2 )), where n is the truncation size, and plateaus at a maximum value of approximately exp(2λ 4 ). For the JFP method f, the condition numbers of (36) asymptote to a value of approximately 1.8483λ 2 .  space expansion (46), Hence, a = C (α,β) −1 ce and b = −λ 2 C (α,β) −1 co, where a and b are the sum space coefficients and the entries of ce and co are ce(n) = 2λ 4 n /Γ(1 + n) and co(n) = 2λ 4 n Γ(3/2 + n) . Since M0 = P (α,β) C (α,β) −1 , the n-th column of C (α,β) −1 represents the Jacobi expansion coefficients of ((1 + x)/2) n , which are bounded below one and strictly positive 16 . Since the entries of ce and co are also strictly positive and their largest entries have a magnitude on the order of O exp(2λ 4 ) , we expect the largest sum space coefficients to have a similar magnitude, which is indeed the case in Fig. 7b. The numerical results in Fig. 7b indicate that the truncation size required for the sum space method to achieve a prescribed accuracy grows as O(λ 4 ). By the same reasoning used above, we also expect this to be the case if the solution is expressed in the normalized monomial basis. Indeed, if we require | uν,µ(n)| = for large λ, with uν,µ(n) defined in (39), then using Stirling's approximation in the limit n → ∞, we deduce that n = O(λ 1/µ ), λ → ∞. Replacing again λ with λ 2 and setting µ = 1/2, the truncation size n grows as O(λ 4 ), as it does for the sum space method in Fig. 7b.
Remark 2. We have considered the problem (44) with order 1/2, however similar conclusions hold if the problem has order 0 < µ < 1: the maximum sum space coefficient grows as O exp(2λ 2/µ ) , and the required precision and truncation size scale as O(λ 2/µ ), which leads to a complexity of O(λ 4/µ log λ log log λ). 16 Numerically, we found that C For the JFP method, the magnitude of the coefficients in Fig. 7c and the O(λ) growth of the truncation size in Figs. 7c and 7d can be derived from the bound (38) on the coefficients 17 . Using (45) and the relation 1 + x = (1 + y) 2 /2 from (17) with p = 2, it follows that u attains its maximum modulus on the Bernstein ellipse Eρ in the complex y-plane at y = −(ρ + ρ −1 )/2 (i.e., where Eρ intersects the negative real axis), hence in the bound (38), Setting ρ = 1 + , M becomes Hence M → 1 as → 0 (or ρ → 1) and we conclude from (38) that the JFP coefficients are bounded above by 2 for any λ, which agrees with the numerical results in Fig. 7c. Next we estimate the growth of the truncation size as a function of λ for the JFP method using the bounds (38) and (50). First, we set 2M ρ −n = c, with c 1 and for fixed λ, n 1, estimate the value of ρ that minimizes 2M ρ −n . We find that the requirement ∂ ∂ρ M ρ −n = 0 leads to the following leading order estimate for the optimal value of ρ, Substituting this ρ into 2M ρ −n = c, we find that hence the truncation size grows linearly in λ.
To make another comparison between the JFP and sum space bases, we take the perspective of orthogonal polynomials on algebraic curves [39,16,15]. In view of (44), in which µ = 1/2 and with p = 2 in (17), we consider orthogonal polynomials on the algebraic curve γ := (x, y) ∈ R 2 : 1 which is shown in Fig. 2a. We define the following inner product on γ, where dσ defines the arc length measure on the curve γ, hence we can set dσ = 1 + (x ) 2 dy, where x = dx dy = y + 1. Setting x = (1 + y) 2 /2 − 1, defining f (y) = f (x, y), g(y) = g(x, y) and choosing w(x, y) 17 The bound (38) is applicable to Chebyshev polynomials, however the JFP results in Fig. 7 were computed with Legendre polynomials (because the condition β − b ∈ N 0 arising from Algorithm 2 rules out the use of Chebyshev polynomials; with Algorithm 1, however, Chebyshev polynomials are permissible). We use the bound (38) since it is simpler than the bounds on Legendre coefficients given in [49]. Furthermore, the first n Chebyshev and Legendre coefficients of the solution to (44) can be related as follows, u 0,0,0,2 (1 : n) = Ru −1/2,−1/2,0,2 (1 : n), where R is an upper triangular matrix with [20] Hence, the Legendre coefficients are at worst larger than the Chebyshev coefficients by a factor of √ n, which is inconsequential since the Chebyshev coefficients decay superexponentially.
such that w(x, y)dσ = w α,β (y)dy, where the Jacobi weight is defined in (8), the inner product on the curve γ becomes the standard Jacobi inner product in y, i.e., Hence, the JFP basis is a weighted orthogonal polynomial basis on the curve γ (for (44), the JFP basis is not weighted because we set b = 0). The sum space basis (47), however, is not orthogonal with respect to the inner product (52) for any choice of α, β > −1. For example, setting α = β = 0 and 1 + x = (1 + y) 2 /2 in (52), we find that m γ,w = 0. The (non-orthogonal) sum space basis is an example of a frame [1, Example 3], which is a generalization of the notion of a basis. Finally, we note that the JFP and sum space bases for the problem (44) are related as follows S(x) = P (α,β) (y)R, where 1 + x = (1 + y) 2 /2, hence R = P (α,β) 2 −1 P (α,β) , S w α,β , which is an upper triangular matrix.
Remark 3. It is worth noting that if k * , and thus p, are large, e.g., if µ = 1/2 and k * = 5 (and thus p = 10), then compared to choosing k * = 1, the rate of convergence is slower for small λ but significantly faster for large λ. The effect of k * and λ on the rate of convergence of the JFP method can likely be understood by using (37) and analyzing the bound (38) as a function of k * and λ.

Example 4.
We consider a periodic one-dimensional time-fractional heat/wave equation of Caputo type, where f (x) has a Fourier series f (x) = +∞ n=−∞ fne inx . For 0 < µ < 1, the equation is known as a fractional heat (or diffusion) equation. For 1 < µ < 2, (53) is referred to as a fractional wave equation and, in addition, ∂tu needs to be specified at t = 0. For simplicity, we shall set ∂tu(x, 0) = 0.
Recalling from (1) that D µ C = I 1−µ D for 0 < µ ≤ 1 and D µ C = I 2−µ D 2 for 1 < µ ≤ 2, we apply I µ to the FDEs (54) to convert them to FIEs: where we have used the conditions un(0) = 1, u n (0) = 0 (for 1 < µ ≤ 2), the latter of which originates from setting ∂tu(x, 0) = 0. Adapting Theorem 10 to the interval [0, T ], it follows that un(t) = Eµ,1(−n 2 t µ ). In practice, we consider (55) for 0 ≤ n ≤ N , where N is the truncation size of the Fourier series that is required to approximate f to a given accuracy. Using a change of variables in (55), or using the exact solution, it can be verified that un(t) = uN This shows that in order to compute solutions to the fractional PDE (53), we need to compute solutions to the FIE considered in the previous example, (44), in the large-λ regime.
As an example, in Figs. 9 and 10 we let u(x, 0) = f (x) = e − cos 2x+1/2 sin x − 2 sin sin x, for which the truncation size of its Fourier series is N = 29 in double precision, and for all values of µ, we let p = 5.

More general FIEs and FDEs
In this subsection, we consider problems whose exact solutions are not known to us and include multiple orders (Example 5), variable coefficients (Example 6) and non-trivial boundary conditions (Example 7).
where it is assumed that µ1 > 0 and µ k , 2 ≤ k ≤ N are positive rational multiples of µ1, hence µ k = m k n k µ1, where the m k and n k are positive (and relatively prime) integers. We can set µ1 = k * /p, where k * ∈ N+ and p > 0 and we further assume that there exists a positive integer n such that the given functions a0(x), f (x), a k (x), b k (x), 1 ≤ k ≤ N have convergent expansions in non-negative powers of (1 + x) 1/(np) , where n ≥ max{n2, . . . , nN }. We allow f (x) to have an expansion in the basis 18 M 1/(np) ν−1 , ν > 0, in which case we divide by (1 + x) ν−1 in (64) and we relabel (1 + x) 1−ν u(x) as u(x). Hence, we assume without loss of generality that f has an expansion in non-negative powers of (1 + x) 1/(np) . Then the solution to (64) has an expansion in the JFP basis Q (α,β,0,np) .
Remark 5. We have demonstrated that the sum space method for the FIE (44) performs poorly compared to the JFP method as λ increases. However, as illustrated in [25], the sum space method converges exponentially fast in linear complexity (because it yields banded or almost-banded systems) and in double precision for a wide range of problems (examples of which include (56) and (60)). For these problems the sum space method is preferable to the JFP method because it converges at a similar rate but with lower complexity (because, by comparison, the JFP method always leads to a lower banded system for fractional order problems) and without the need to compute integration matrices in high precision. The sum space method has linear complexity if the variable coefficients are analytic functions of x but for an FIE such as (58), the method also leads to a lower banded system.
The sum space method is only applicable to equations of rational order µ = p/q and requires the direct sum of 2q different weighted orthogonal polynomial bases 19 . For example, for (43), q = 6 because the solution has an expansion in powers of (1 + x) 1/6 , and thus 12 weighted bases are required while the JFP method uses a single basis (but with different parameter values p and b) for all fractional order (including irrational order) problems.
We have not attempted to make a rigorous classification of the types of problems for which the JFP method is superior to the sum space method and vice versa. However, Example 3 suggests that the sum space method performs poorly for problems in which the largest monomial coefficient of the solution becomes large (e.g., on the order of 10 10 ). Otherwise, we conjecture, the sum space method performs well.
Large condition numbers of the systems that arise in the sum space method are not good indicators of the accuracy of the resulting solutions. For example, in [25,Example 11], the following fractional Airy equation is considered, which is a singularly perturbed problem with an increasingly oscillatory solution as → 0. For this problem the sum space method achieves high accuracy (around 10 −10 for = 10 −4 ) despite large condition numbers (on the order of 10 15 and larger). This is reminiscent of the high accuracy achieved by the ultraspherical spectral method in [38] for a singularly perturbed standard Airy equation (i.e., with i 3/2 D 3/2 RL in (65) replaced by D 2 ) despite the large condition numbers that arise. By our conjecture, we expect the sum space method to perform well for the problem (65), large condition numbers notwithstanding, because the monomial coefficients of the solution to (65) do not become large.

Conclusion
We have illustrated the application of the JFP method to a variety of FIEs (including FDEs and a fractional PDE reformulated as FIEs) in which exponentially fast convergence to the solution is achieved. The JFP method converges much faster and with a lower overall complexity than the sparse sum space method in [25] for solutions whose power series (or shifted monomial) coefficients become large. For such problems, a relatively large number of coefficients are needed to resolve the solution in which case the use of highprecision arithmetic is essential to scale up the JFP method while retaining sufficient accuracy. Pseudostabilization of the unstable algorithms we have introduced for fractional integration matrices incorporates high-precision automatically (see Appendix A) and increases the complexity of the JFP method from O(N 2 ) (without high-precision, but unstable) to O(N 3 log N log log N ).

Future work and open problems
We shall construct differentiation matrices for the JFP basis, which will be banded and could be used to solve FDEs and fractional PDEs without having to reformulate them as FIEs, as we have done in this paper. Adding Newton iteration in function space [8], the JFP method would also be applicable to nonlinear FIEs and FDEs.
The superior performance of the JFP method in Example 3 of Section 6.1 suggests it could be an effective method for computing Mittag-Leffler functions. Just as the ultraspherical spectral method was shown in [10] to be an effective method for the global computation of a special function (the Gauss hypergeometric function), the JFP method could be adapted in a similar way to efficiently evaluate Mittag-Leffler functions on intervals and regions in the complex plane.
The convergence and stability analysis of the JFP method for general linear FIEs (64) is a topic for future research. This analysis could prove the bounds on condition numbers that we found in Fig. 7f and reveal the dependence of the rate of convergence on the parameter p in the JFP basis (see remark 3).
Another topic for future research is stable algorithms for fractional integration matrices in the JFP basis with optimal (i.e., O(N 2 )) complexity. As mentioned in Section 5.2, one possibility is to use asymptotic approximations to the entries of the fractional integration matrices to stabilize the Sylvester equations in Algorithm 2.

A.1 Introduction
Pseudo-stabilization is a technique that uses high-precision computations to ensure that the output of an unstable algorithm has a specified accuracy. Pseudo-stabilization is called thus because it makes an unstable algorithm appear stable by outputting accurate results. This technique is algorithm-specific and here we illustrate it for Algorithm 2 (see Section 5.2), which computes the entries of fractional integration matrices via the recurrence relation (28) (with A = I (α,β) b,p,µ ) in which errors grow exponentially. Since high-precision computations are more expensive than standard double precision, our aim is to use as low a precision as possible while ensuring the integration matrices are accurate to the specified tolerance. To choose this optimal precision, we need to know how fast errors grow in (28), examples of which are shown in (14). Since the exact values of the entries of I (α,β) b,p,µ are not known, throughout this appendix, errors in a given precision are computed with reference to results computed in a much higher precision.

A.2 Growth of errors
The error in the first p columns of I (α,β) b,p,µ is determined by the precision used in Algorithm 1, while the growth rate of the error in subsequent columns depends on the recurrence coefficients (i.e., the entries of B and C in (28)) and hence also on the precision with which they are computed, as confirmed by Fig. 14. For recurrence coefficients computed in double precision, Fig. 14b shows that the errors grow roughly at a constant rate for the first circa 150 columns, after which the growth rate starts shooting up. The errors in higher (256-bit) precision grow at roughly a constant rate in Fig. 14b, however we shall find in Fig. 16 that if a sufficiently large number of columns are computed, the growth rate of the error also increases (for any precision). Due to the complicated behavior of the growth rate of the error, we shall adopt an empirical approach to model and predict the growth of the error.
Let q be the precision and the corresponding machine epsilon, then (q) = 2 1−q . We let e0 and E(q, n; e0) denote, respectively, the maximum errors in the initial p columns and in the n-th column of I  (28)) were computed with 256-bit precision (see Table 2). In Fig. 14b, the growth rate is the maximum error in a column divided by that in the previous column.
with n ≥ 0. Fig. 14a shows examples of E(q, n; e0) for q = 53 and q = 256. We let E(q, n) = E(q, n; e0)/e0, hence, by definition, E(q, p − 1) = 1 and therefore we simulate the normalized maximum error E(q, n) by E(q, n), which is defined as the maximum value in the n-th column of A when one or more of the initial values in the recurrence (28) are set to 1 and the rest are set to zero. We have found that E(q, n; e0) e0 = E(q, n) ≈ E(q, n) holds, see Fig. 15 for an example. This implies that we can model the growth of errors with E(q, n).
Next we consider the estimated growth rate of the error, E(q, n; e0) E(q, n − 1; e0) ≈ E(q, n) E(q, n − 1) := r(q, n). Fig. 14b shows r(q, n) for q = 53, 256 and Fig. 16 shows r(q, n) for much larger n (up to 10 4 ). In Fig. 14b, r(256, n) is approximately constant, however in Fig. 16 it is seen that r(256, n) also increases in a step-like manner for large enough n. Fig. 16b shows that r can be very noisy, however the denoised or averaged version of r, which we denote by r0(q, n), takes the form of a superposition of step functions satisfying the scaling property (see Figs. 16a and 16c) r0(q, n) ≈ r0(λq, λn).
which means that we can use a simulation on lower precision and fewer iterations to predict error propagation on higher precisions through a larger number of iterations. The approximation (70) is exact for the simplest case in which r(q, n) and r(λq, λn) are constant and equal (e.g., for double precision with q = 53, n ≤ 150, see Fig. 14b, and λ > 1). If we let e0(q) = 2 −q , then for constant r, E(n, q; e0) = r n 2 −q and E(λn, λq; e0) = r λn 2 −λq = E λ (n, q; e0). Remark 6. In Fig. 16, the height of the steps of r0 for different q are the same but the lengths are proportional to q (hence the property (68)). Except for the first step, which is the longest for every q, every step has the same length. We are surprised by the step-like behavior of the growth rate. It seems to be caused by the errors in computing the integration matrix I (α,β) b,p . We conjecture that if we could compute the integration matrix exactly, then the first step of r(q, n) would be infinitely long, i.e., the growth rate would be constant, which is what we would expect for a linear recurrence problem.
Remark 7. We have found that even though the approximations (68) and (70) are rather crude, especially when r is noisy, there is enough 'noise cancellation' in the products of growth rates in (69) (with r0 replaced by the actual growth rates r defined in (67)) for the approximation (70) to yield useful estimates.

A.3 Implementation
We now use the empirical estimate (70) for the growth of errors to estimate the lowest precision required to meet a given tolerance. These estimates were used to compute the fractional integration matrices that appear in the problems in Section 6.
Suppose we aim to achieve a maximum error of at most δ in the computation of the first N + 1 columns of A = I (α,β) b,p,µ using the recurrence (28). Setting e0 = 2 −q and using (67), we seek a q such that E(q, N ; e0) = 2 −q E(q, N ) < δ.
Now the idea is to simulate E with E on a chosen lower precision q0 and a smaller number of iterations m < N (which is to be determined) and use (70) to estimate the smallest q satisfying (71). Using (70) with λ = q/q0, we have E(q, N ) ≈ E(q0, N q0/q) q/q 0 , hence we require E(q0, N q0/q) < 2 q 0 δ q 0 /q .
Setting m = N q0/q and since we aim to minimize q, we seek the largest m (say m * ) such that Typical parameter choices are q0 = 53 and δ = 10 −16 which implies results in double precision, while N , the polynomial degree in the JFP basis required to resolve the solution to a given accuracy (which is typically the same as δ), is problem-dependent. To estimate m * in practice, we simulate E(q0, m) with E(q0, m) and increase m until the first value of m (say m1) is found for which (72) is not valid, then we set m * = m1 − 1 and set the precision to q = N q0/m * .

A.4 Computational complexity
Now we analyse the cost of the pseudo-stabilized Algorithm 2 with inputs p, N and δ 1, where p ∈ N+ satisfies µp = k * ∈ N+. Using (72) and (70), it follows that log 2 E(1, m/q0) < 1, so m = O(q0) and since q = q 0 m N , q = O(N ). The first N columns of I (α,β) b,p,µ involve the computation of O(N 2 ) entries, each of which requiring up to 4p + 3 multiplications, leading to a total cost of O(pN 2 q log q log log q) in q-bit precision 20 . Recalling that q = O(N ), we conclude that Algorithm 2 has an overall complexity of O(pN 3 log N log log N ).
Remark 8. Since the accuracy drops during the recurrence, the precision can be lowered gradually to reduce the complexity by a constant factor. However, that does not change the conclusions with big O notations.
Remark 9. There are other computational costs that are negligible compared to the main recurrence process: • There are divisions involving entries of the banded integration matrices which means that we can compute the inverse first with a cost of O(pN 2 log N (log log N ) 2 ) using Newton-Raphson division.
• The computation of E(q0, m) has O(pq 3 0 log q0 log log q0) complexity. Remark 10. It is also possible to pseudo-stabilize Algorithm 1 (see Section 5) since the growth of errors follow similar patterns. However, the complexity of Algorithm 1 is higher, making the overall complexity of the pseudo-stabilized version O(N 4 ). We therefore only use Algorithm 1 if µ is irrational (recall that Algorithm 2 is only applicable for rational µ).