A sparse spectral method for Volterra integral equations using orthogonal polynomials on the triangle

TL;DR

This paper presents a sparse spectral method leveraging bivariate orthogonal polynomials on a triangle for efficiently solving Volterra integral equations with proven exponential convergence.

Contribution

The paper introduces a novel sparse spectral approach using orthogonal polynomials on a triangle domain for Volterra equations, demonstrating high efficiency and convergence.

Findings

Achieves exponential convergence in solving Volterra equations

Demonstrates effectiveness on first and second kind equations

Proves convergence leveraging Toeplitz operator connections

Abstract

We introduce and analyse a sparse spectral method for the solution of Volterra integral equations using bivariate orthogonal polynomials on a triangle domain. The sparsity of the Volterra operator on a weighted Jacobi basis is used to achieve high efficiency and exponential convergence. The discussion is followed by a demonstration of the method on example Volterra integral equations of the first and second kind with known analytic solutions as well as an application-oriented numerical experiment. We prove convergence for both first and second kind problems, where the former builds on connections with Toeplitz operators.

Equations

$$({\mathcal{V}}_{K}u)(x):={\int }_0^{l(x)}K(x,y)u(y)\mathrm{d}y,$$

$${\mathcal{V}}_{K}u=gor(I+{\mathcal{V}}_K)u=g.$$

$$⟨P_{m,k},P_{n,j}⟩={\int }_{\Omega }{P}_{m,k}(\mathbf{x})P_{n,j}(\mathbf{x})W(\mathbf{x})\mathrm{d}A=C{\delta }_{mn}{\delta }_{jk},$$

$${\int }_{-1}^1{C}_{(\alpha ,\beta ,m,n)}{\left(1-x\right)}^{\alpha }{\left(1+x\right)}^{\beta }{P}_m^{(\alpha ,\beta )}(x)P_n^{(\alpha ,\beta )}(x)\mathrm{d}x={\delta }_{nm},$$

$$f(x)=\sum _{n=0}^{\infty }{p}_n(x)f_n=\mathbf{P}(x{)}^{\mathsf{T}}\mathbf{f},$$

$$\mathbf{P}(x):=\left(\begin{array}{c}{p}_0(x)\\ p_1(x)\\ ⋮\end{array}\right)$$

$$\mathbf{P}(x{)}^{\mathsf{T}}\mathrm{J}{\mathbf{f}}_{[0,1]}=xf(x).$$

$$\mathrm{J}=\left(\begin{array}{cccc}{a}_0& b_0& & \\ c_0& a_1& b_1& \\ & c_1& a_2& \ddots \\ & & \ddots & \ddots \end{array}\right).$$

$$xf(x)$$

$$(1-x)f(x)$$

$$T^2=\left\{(x,y):0\le x,0\le y\le 1-x\right\}.$$

$$f(x,y)=\sum _{n=0}^{\infty }\sum _{k=0}^n{p}_{n,k}(x,y)f_{n,k}.$$

$${\mathbb{P}}_n(x,y)=\left(\begin{array}{c}{p}_{n,0}(x,y)\\ p_{n,1}(x,y)\\ ⋮\\ p_{n,n}(x,y)\end{array}\right)$$

$$\mathbf{P}(x,y)=\left(\begin{array}{c}{\mathbb{P}}_0(x,y)\\ {\mathbb{P}}_1(x,y)\\ ⋮\end{array}\right).$$

$$f(x,y)=\sum _{n=0}^{\infty }\sum _{k=0}^n{p}_{n,k}(x,y)f_{n,k}=\mathbf{P}(x,y{)}^{\mathsf{T}}\mathbf{f}.$$

$$P_{k,n}^{(\alpha ,\beta ,\gamma )}(x,y)={\left(1-x\right)}^k{P}_{n-k}^{\left(2k+\beta +\gamma +1,\alpha \right)}\left(2x-1\right)P_k^{\left(\gamma ,\beta \right)}\left(\frac{2y}{1-x}-1\right).$$

$${\int }_0^1{\int }_0^{1-x}{x}^{\alpha }{y}^{\beta }(1-x-y{)}^{\gamma }{P}_{k,n}^{(\alpha ,\beta ,\gamma )}(x,y)P_{j,m}^{(\alpha ,\beta ,\gamma )}(x,y)\mathrm{d}y\mathrm{d}x=C_{(\alpha ,\beta ,\gamma )}{\delta }_{jk}{\delta }_{mn}.$$

$$P_{k,n}^{(\alpha ,\beta ,\gamma )}(x,y)=(1-x{)}^k{\tilde{P}}_{n-k}^{\left(2k+\beta +\gamma +1,\alpha \right)}\left(x\right){\tilde{P}}_k^{\left(\gamma ,\beta \right)}\left(\frac{y}{1-x}\right).$$

$$\mathbf{P}(x,y{)}^{\mathsf{T}}{\mathrm{J}}_x{\mathbf{f}}_{△}$$

$$\mathbf{P}(x,y{)}^{\mathsf{T}}{\mathrm{J}}_y{\mathbf{f}}_{△}$$

$${\mathrm{J}}_x=\left(\begin{array}{cccc}{A}_0^x& B_0^x& & \\ C_0^x& A_1^x& B_1^x& \\ & C_1^x& A_2^x& \ddots \\ & & \ddots & \ddots \end{array}\right),{\mathrm{J}}_y=\left(\begin{array}{cccc}{A}_0^y& B_0^y& & \\ C_0^y& A_1^y& B_1^y& \\ & C_1^y& A_2^y& \ddots \\ & & \ddots & \ddots \end{array}\right),$$

$${\mathcal{L}}_N(x_{\ast }){\mathbf{P}}_N^{(\alpha ,\beta )}(x_{\ast })=\left(\begin{array}{cccccc}1& & & & & \\ a_0-x_{\ast }& b_0& & & & \\ c_0& a_1-x_{\ast }& b_1& & & \\ & \ddots & \ddots & \ddots & & \\ & & c_{N-2}& a_{N-1}-x_{\ast }& b_{N-1}& \end{array}\right)\left(\begin{array}{c}{P}_1^{(\alpha ,\beta )}(x_{\ast })\\ P_2^{(\alpha ,\beta )}(x_{\ast })\\ P_3^{(\alpha ,\beta )}(x_{\ast })\\ ⋮\\ P_N^{(\alpha ,\beta )}(x_{\ast })\end{array}\right)={\mathbf{e}}_0,$$

$$f(x_{\ast })={\mathbf{P}}_N^{(\alpha ,\beta )}(x_{\ast }{)}^{\mathsf{T}}\text{a}={\mathbf{e}}_0^{\mathsf{T}}{\mathcal{L}}_N(x_{\ast }{)}^{-\mathsf{T}}\text{a},$$

$${\mathcal{L}}_N({\mathbf{x}}_{\ast }){\mathbf{P}}_N^{(\alpha ,\beta ,\gamma )}({\mathbf{x}}_{\ast })=\left(\begin{array}{ccccc}{1}_1& & & & \\ A_0^x-x_{\ast }{1}_1& B_0^x& & & \\ A_0^y-y_{\ast }{1}_1& B_0^y& & & \\ C_0^x& A_1^x-x_{\ast }{1}_2& B_1^x& & \\ C_0^y& A_1^y-y_{\ast }{1}_2& B_1^y& & \\ & \ddots & \ddots & \ddots & \end{array}\right){\mathbf{P}}_N^{(\alpha ,\beta ,\gamma )}({\mathbf{x}}_{\ast })={\mathbf{e}}_0,$$

$$\left(\begin{array}{cccc}1& & & \\ & B_0^{+}& & \\ & & B_1^{+}& \\ & & & \ddots \end{array}\right){\mathcal{L}}_N({\mathbf{x}}_{\ast })={\tilde{\mathcal{L}}}_N({\mathbf{x}}_{\ast }).$$

$$B_n=\left(\begin{array}{c}{B}_n^x\\ B_n^y\end{array}\right),$$

$$f({\mathbf{x}}_{\ast })={\mathbf{P}}_N^{(\alpha ,\beta ,\gamma )}({\mathbf{x}}_{\ast }{)}^{\mathsf{T}}\text{A}={\mathbf{e}}_0^{\mathsf{T}}{\tilde{\mathcal{L}}}_N({\mathbf{x}}_{\ast }{)}^{-\mathsf{T}}\text{A}.$$

$$\mathbf{P}(x{)}^{\mathsf{T}}{\mathrm{W}}_{\mathrm{Q}}{\mathrm{Q}}_y{\mathbf{f}}_{△}={\int }_0^{1-x}f(x,y)\mathrm{d}y,$$

$$\mathbf{P}(x{)}^{\mathsf{T}}{\mathbf{f}}_{[0,1]}=\mathbf{P}(x,y{)}^{\mathsf{T}}{\mathrm{E}}_y{\mathbf{f}}_{[0,1]}$$

$${\int }_0^{1-x}f(y)\mathrm{d}y=\mathbf{P}(x{)}^{\mathsf{T}}{\mathrm{W}}_{\mathrm{Q}}{\mathrm{Q}}_y{\mathrm{E}}_y{\mathbf{f}}_{[0,1]}$$

$${\tilde{P}}_n^{(\alpha ,\beta )}(x)=(-1{)}^n{\tilde{P}}_n^{(\beta ,\alpha )}(1-x).$$

Full text

∗Department of Mathematics, Imperial College London, UK. (t.gutleb18@imperial.ac.uk)

*‡*Department of Mathematics, Imperial College London, UK. (s.olver@imperial.ac.uk)

1. Introduction

Define the Volterra integral operator

[TABLE]

where $K(x,y)$ is called the kernel, $u(y)$ is a given function of one variable and the limits of integration are either $l(x)=x$ or $l(x)=1-x$. This paper concerns Volterra integral equations of the first and second kind, that is, to find $u$ satisfying

[TABLE]

Numerous applications and the fundamental nature of Volterra integral and integro-differential equations motivate research into efficient and accurate numerical solvers. Various forms of Volterra integral equations are analytically well-understood [12, 37, 47], have been the subject of various numerical approximation schemes [12, 11, 5, 30], and are encountered regularly in various scientific fields as well as engineering and finance applications [12, 37, 45, 47, 25, 26].

In this paper we present a method to compute Volterra integrals and solve Volterra integral equations by using orthogonal polynomials on a triangle domain [19, 36] to both resolve the kernel and to reduce the equations to banded linear systems. The method is in the same spirit as some previous contributions to the field of numerical Volterra, Fredholm, singular integral and differential equations based on operators and orthogonal polynomials such as [1, 24, 41, 23] but differs in choice of basis and domain, leading to operator bandedness properties which can be exploited for significantly increased efficiency. Notably the approach introduced in this paper can be used for a wider range of kernels than many other Volterra integral equation solvers such as the methods based on orthogonal polynomials due to Loureiro and Xu [29, 50], the recently developed ultraspherical spectral method in [23] or the Fourier extension method in [49] as it is not limited to convolution kernel cases, that is kernels of the form $K(x,y)=K(x-y)$, but works for a wider class of kernels.

The sections in this paper are organized as follows: Section 2 introduces the required aspects of univariate and bivariate polynomial function approximation on a real interval and the triangle respectively. Section 3 introduces an efficient numerical method for Volterra integrals and integral equations and discusses how to approach kernel computations using a multivariate variant of Clenshaw’s algorithm. In Section 4 we show the scheme in action in both toy and application-based examples. Proofs of convergence for well-posed problems are discussed in Section 5.

2. Function approximation with orthogonal polynomials

2.1. Jacobi polynomials on the real interval

Multivariate orthogonal polynomials are ordered sets of polynomials satisfying a particular pair-wise and weighted orthogonality condition, often of the form

[TABLE]

where $C\neq 0$ and $P_{m,k}$ are total degree $m$ polynomials. Many such sets of orthogonal polynomials are well-known and well-studied on various domains $\Omega$ such as $\mathbb{R}$, real intervals, simple $2$D and $3$D domains, as well as various higher dimensional spheres and polygons [19]. The relevant set of orthogonal polynomials for this paper are the Jacobi polynomials on the real line and on the triangle respectively. This section will thus give a quick overview of Jacobi polynomials aimed at equipping us with the tools needed to develop the Volterra integral equation solvers in later sections. We refer to [19, 20] for introductions with broader scope.

The Jacobi polynomials are orthogonal on $[-1,1]$:

[TABLE]

where $W_{(\alpha,\beta)}(x)=C_{(\alpha,\beta,m,n)}\left(1-x\right)^{\alpha}\left(1+x\right)^{\beta}$ acts as the weight function and $\delta_{nm}$ is the Kronecker delta. While the choice of $[-1,1]$ is natural, the Jacobi polynomials can be shifted to any real interval an application requires. For $\alpha=\beta=0$ the Jacobi polynomials reduce to the Legendre polynomials [19].

One of the primary applications of interest for the study of orthogonal polynomials are their applications in the expansion of non-polynomial functions:

[TABLE]

where $f_{n}$ is the function-specific coefficient of the $n$-th polynomial $p_{n}$ and we use the notation

[TABLE]

For numerical applications one uses finitely many terms in the above sum to obtain an approximation. If a distinction between different sets of polynomials and coefficient vectors on different domains is required we specify by indicating the type of polynomials using standard notation for the polynomials, such as $\mathbf{P}^{(\alpha,\beta)}(x)$ for the Jacobi polynomials on a real interval, and the domain using index notation, e.g. for the bivariate orthogonal polynomial coefficient vector of $g(x,y)$ on the triangle domain we write $\mathbf{g}_{\vartriangle}$.

To use function approximation of this type in a non-trivial numerical application one needs ways to do computations on functions represented as coefficient vectors. Basic computations such as addition and subtraction of functions have obvious implementations. Furthermore one can compute $xf(x)$ if $f(x)$ is already approximated as a coefficient vector: to do this one uses so-called Jacobi operators $\mathrm{J}$ which act as

[TABLE]

This is efficiently possible because the Jacobi polynomials satisfy a three-term recurrence relationship, making $\mathrm{J}$ a tridiagonal operator (see e.g. [19, 32, 36]):

[TABLE]

Additionally, our approach to Volterra integral equations of the second kind will require explicit constructors for raising operators $\mathrm{S}_{(\alpha,\beta)}^{(\alpha+1,\beta)},\mathrm{S}_{(\alpha,\beta)}^{(\alpha,\beta+1)}$ which are defined to increment from the Jacobi bases $\mathbf{P}^{(\alpha,\beta)}(x)$ to $\mathbf{P}^{(\alpha+1,\beta)}(x)$ and $\mathbf{P}^{(\alpha,\beta+1)}(x)$ respectively. Increments to $\alpha$ and $\beta$ can be computed using these operators but decrementing is generally only well-defined in the sense of weighted lowering operators:

[TABLE]

The explicit forms of the operators $\mathrm{J}$, $\mathrm{S}_{(\alpha,\beta)}^{(\alpha+1,\beta)}$, $\mathrm{S}_{(\alpha,\beta)}^{(\alpha,\beta+1)}$, $\mathrm{L}_{(\alpha,\beta)}^{(\alpha-1,\beta)}$ and $\mathrm{L}_{(\alpha,\beta)}^{(\alpha,\beta-1)}$ are well known in the literature, see for example [32, 36, 19] and the references therein.

2.2. Jacobi polynomials on the triangle

We now briefly discuss how function approximation using bivariate orthogonal polynomials works in general and then move on to discuss the Jacobi polynomials on the canonical unit simplex

[TABLE]

We use a basis on this triangle in the following sections to compute Volterra integrals and solve integral equations. As in the univariate case, bivariate orthogonal polynomials are said to be orthogonal with respect to an inner product akin to (2).

Analogously to how functions of a single variable may be expanded into a basis of univariate orthogonal polynomials as $f(x)=\sum_{n=0}^{\infty}p_{n}(x)f_{n}$ we can expand a function of two variables in a basis of bivariate polynomials as

[TABLE]

Writing the bivariate polynomials of total degree $n$ as

[TABLE]

allows for the following compact notation for the infinite-dimensional polynomial basis:

[TABLE]

In this notation the expansion of a function of two variables in the bivariate polynomial basis becomes

[TABLE]

For function approximation one simply uses an appropriate finite cutoff of this expansion.

On the triangle $T^{2}$ we focus on the Jacobi weights $x^{\alpha}y^{\beta}(1-x-y)^{\gamma}$. One elegant way to define the corresponding Jacobi polynomials $\mathbf{P}^{(\alpha,\beta,\gamma)}(x,y)$ on the canonical triangle $T^{2}$ is by referring to the Jacobi polynomials $\mathbf{P}^{(\alpha,\beta)}(x)$ on the real interval $[-1,1]$ (compare [19, Proposition 2.4.1]):

[TABLE]

Defined as such the triangle Jacobi polynomials are orthogonal with respect to a weighted integral over the canonical triangle domain $T^{2}$:

[TABLE]

The detailed form of the constant $C_{(\alpha,\beta,\gamma)}$ is not important here but can for example be found in [19]. We will primarily use the Jacobi polynomials shifted to the $[0,1]$ interval and denote them by $\tilde{\mathbf{P}}^{(\alpha,\beta)}$, which allows us to write the Jacobi polynomials on the triangle as:

[TABLE]

As in the 1-dimensional case we can define Jacobi operators $\mathrm{J}_{x}$ and $\mathrm{J}_{y}$, one for each variable, which respectively act as

[TABLE]

for a given bivariate polynomial basis. Unlike the 1-dimensional Jacobi polynomial case these operators are not tridiagonal but block tridiagonal operators for the triangle Jacobi polynomials [36]:

[TABLE]

where $A_{n}^{x},A_{n}^{y}\in\mathbb{R}^{(n+1)\times(n+1)}$, $B_{n}^{x},B_{n}^{y}\in\mathbb{R}^{(n+1)\times(n+2)}$ and $C_{n}^{x},C_{n}^{y}\in\mathbb{R}^{(n+2)\times(n+1)}$. Analogous operators to the raising and lowering operators discussed for the real interval case can be constructed for the Jacobi polynomials on the triangle as well, see [35, 36], but we omit their discussion as we will not make direct use of them in this paper.

To make use of Jacobi polynomials for the approximation of functions on the triangle domain in a numerical context one requires efficient algorithms to determine the coefficient vector $\mathbf{f}_{\vartriangle}$ for a given function $f(x,y)$ of two variables. This can be done using an algorithm and its implementation in a C library by Slevinsky [38, 39, 40].

2.3. Function evaluation using Clenshaw’s algorithm

Clenshaw’s algorithm provides an efficient and direct method to evaluate functions expanded into orthogonal polynomial bases at given points, i.e. to evaluate $\sum_{n=0}^{N}p_{n}(\mathbf{x})f_{n}$ at $\mathbf{x}_{*}\in\mathbb{R}^{d}$, cf. [15, 36]. The algorithm makes use of the polynomial basis’ recurrence relationships to reduce function evaluation to the solution of an upper triangular linear system using backward substitution. In this section we give an outline of how this is done for Jacobi polynomials on the real interval and the triangle, which is discussed in more detail in [36]. An operator valued variant of what is discussed in this section will be used for efficient kernel computations for Volterra integrals in section 3.2. We mention a major benefit of Clenshaw’s algorithm over building polynomials/operators via forward recurrences is that there is substantially less memory needed in the intermediary calculations.

For the case of Jacobi polynomials on a real interval, the three-term recurrence relationship seen in the Jacobi operator in (3) can be used to write

[TABLE]

where $\mathbf{e}_{0}$ is the first standard basis vector with $1$ in its first component and of appropriate length. Solving this lower triangular system via forward substituition provides a way to recursively evaluate each component of $\mathbf{P}^{(\alpha,\beta)}(x)$ and thus also $\mathbf{P}^{(\alpha,\beta)}(x)^{\mathsf{T}}\mathbf{f}$ if the coefficients of $f(x)$ in this basis are known. Clenshaw’s algorithm is conceptually similar but uses backward substition on the system

[TABLE]

where a is the column vector collecting $a_{0}$ to $a_{N}$. The case for the Jacobi polynomials on the triangle was recently discussed in [36] and on the basis of the recurrence in (6) involves a block triangular system for evaluation at $\mathbf{x}_{*}=(x_{*},y_{*})$ instead:

[TABLE]

where $\mathbb{1}_{k}$ denotes the $k\times k$ identity matrix. As this is not a triangular but a block triangular matrix one cannot use forward substitution without first applying a preconditioner:

[TABLE]

$\tilde{\mathcal{L}}_{N}(\mathbf{x}_{*})$ is then a proper lower triangular matrix and can be used in an analogous system to the ones above to evaluate the polynomials, and thus a function expanded into that polynomial basis, recursively via forward substitution. A preconditioner which satisfies these requirements is the block diagonal matrix whose elements are comprised of a left inverse of the blocks

[TABLE]

such that $B_{n}^{+}B_{n}=\mathbb{1}_{n}$. Clenshaw’s algorithm for the triangle Jacobi polynomials is thus

[TABLE]

This system can be solved via backward substitution in optimal $O(N^{2})$ complexity if one chooses $B_{n}^{+}$ carefully, see [36].

3. A numerical method for Volterra integral equations

3.1. Volterra integrals on the triangle

In this section we describe how to represent Volterra integrals using bivariate orthogonal polynomials on a triangle domain by moving to a view of operators acting on coefficient vectors. The following section extends this method to Volterra integral equations of the first and second kind.

We first describe the idea behind the relevant operators and their use before determining their entries in matrix representation. The first operator we need is the integration operator for a function given as the coefficients of orthogonal polynomials on a triangle. We label this operator $\mathrm{Q}_{y}$ and it acts as

[TABLE]

where $\mathrm{W}_{\mathrm{Q}}$ is a to-be-determined weight function which depends on the used basis. The reason for the limits of integration to be defined in this way for $\mathrm{Q}_{y}$ will become clear once we discuss the explicit form of these operators and how one can make optimal use of the triangle domain’s symmetries. Second, we need an operator $\mathrm{E}_{y}$ which extends a one-dimensional function on $[0,1]$ to one on $T^{2}$, that is:

[TABLE]

Together these two operators can be used to compute integrals of the form

[TABLE]

with function $f$ depending on a single variable. To instead integrate from [math] to $x$ we use a reflection operator. Due to symmetries of the polynomials, particular basis changes in a Jacobi basis obey the simple rule [32, 19]:

[TABLE]

We use $\mathrm{R}$ to refer to the operator that uses the above property to reflect the function on the $[0,1]$ interval via a basis change, i.e.

[TABLE]

$\mathrm{J}_{x}$ and $\mathrm{J}_{y}$ have important commutation relations with the introduced $\mathrm{Q}_{y}$ and $\mathrm{E}_{y}$ operators. As the $\mathrm{Q}_{y}$ operator integrates with respect to $y$ and collapses a bivariate coefficient vector back to a univariate one the multiplication-with-$x$ operator changes from being multiplication-with-$x$ on the triangle ($=\mathrm{J}_{x}$) to being multiplication-with-$x$ on the real interval ($=\mathrm{J}$) when pulled through the $\mathrm{Q}_{y}$ operator. A similar relation holds for similar reasons for $\mathrm{J}_{y}$ and $\mathrm{E}_{y}$:

[TABLE]

We now give the explicit matrix representations for the operators $\mathrm{Q}_{y}$ and $\mathrm{E}_{y}$ and discuss a sensible polynomial basis choice. The explicit form of the Jacobi operators on the real line is known in the literature (e.g. [19, 36]) and thus receives no further discussion here. To determine the explicit form of $\mathrm{Q}_{y}$ we begin by plugging in the polynomial expansion of $f(x,y)$ into the intended integral operation and using the Jacobi polynomials on the triangle domain as seen in (5) for our basis $p_{n,k}$ with $\alpha=\beta=\gamma=0$:

[TABLE]

where a substitution of $\frac{y}{1-x}\rightarrow s$ was made in the last step. As $\tilde{P}_{k}^{(0,0)}$ are just the Legendre polynomials on $[0,1]$ we see that $\int_{0}^{1}\tilde{P}_{k}^{(0,0)}\left(s\right)\mathrm{d}s=0,\forall k>0$ and $\int_{0}^{1}\tilde{P}_{0}^{(0,0)}\left(s\right)\mathrm{d}s=1$, resulting in

[TABLE]

for integration from [math] to $1-x$. Via (9) we further obtain

[TABLE]

for integration from [math] to $x$. This derivation shows that starting in the Jacobi polynomial basis on the triangle $T^{2}$ with $\alpha=\beta=\gamma=0$ for the approximation of $f(x,y)$ results in the following block diagonal structure for the integration from [math] to $1-x$ operator with weight $\mathrm{W}_{\mathrm{Q}}=(1-x)$:

[TABLE]

where the $n$-th block is an $n$-dimensional row vector with $1$ in the first element and [math] in all remaining elements. An additional $(-1)^{n}$ term and change of basis changes this integration to be from [math] to $x$ instead. The expansion operator $\mathrm{E}_{y}$ from the $\mathbf{P}^{(1,0)}(x)$ basis to the canonical triangle Jacobi polynomials where $\alpha=\beta=\gamma=0$ has the block diagonal structure

[TABLE]

where the $n$-th block is an $n$-dimensional column vector whose $j$-th entry is given by

[TABLE]

Importantly, multiplication of $\mathrm{Q}_{y}$ and $\mathrm{E}_{y}$ yields a diagonal matrix whose $n$-th entry can be directly generated without any matrix multiplication being required (compare [32]):

[TABLE]

These observations justify the basis choices as well as the choice of the limits of integration for $\mathrm{Q}_{y}$ from the standpoint of computational efficiency. Defining $\mathrm{Q}_{y}$ as the integration operator from [math] to $x$ does not avoid the reflection step and only results in a less efficient or equivalent placement for it.

3.2. Kernel computations using Clenshaw’s algorithm

Putting all the above observations together means one can save a significant amount of computation time by the use of a recurrence when simultaneously using an operator valued polynomial approximation for the kernel $K(\mathrm{J}_{x},\mathrm{J}_{y})$ and then using the known commutation relations in (10–11). To illustrate the idea behind this approach we first discuss how to do this for a monomial kernel (or equivalently a kernel approximated in a monomial basis) and then show how these ideas can be expanded to arbitrary polynomial bases for the kernel using a variant of Clenshaw’s algorithm.

Assuming a monomial expansion for the kernel, i.e. $K(x,y)=\sum_{n=0}^{\infty}\sum_{j=0}^{n}k_{nj}x^{n-j}y^{j}$, the primary part of the Volterra integration operator has the form

[TABLE]

where we have used the commutation relations in (10–11) to rewrite the summation using the Jacobi operator for the interval Jacobi polynomials. Recalling that $\mathrm{Q}_{y}\mathrm{E}_{y}$ is a diagonal matrix which can be generated without any need to separately compute and multiply $\mathrm{Q}_{y}$ and $\mathrm{E}_{y}$, all that is left to compute are the required combinations of $\mathrm{Q}_{y}\mathrm{E}_{y}$ with the Jacobi operators, which can be built up recursively. This kind of recursive computation of all the required elements for the kernel can save significant computation cost if executed correctly. Since only the coefficients of $K(x,y)$ for this basis actually change across different problems one can in principle also store the basis elements $\mathrm{J}^{n-j}\mathrm{Q}_{y}\mathrm{E}_{y}\mathrm{J}^{j}$ and re-use them making this numerical evaluation of Volterra integrals even faster upon repeated use. This approach differs slightly depending on whether one intends to compute integrals from [math] to $1-x$ or to compute integrals from [math] to $x$. In the case of integrals from [math] to $x$, one is either required to supply $K(1-x,y)$ to the algorithm or alternatively the Jacobi operators on the left can be replaced by $(\mathbb{1}-\mathrm{J})$ to account for the reflection, meaning that the basis elements become $(\mathbb{1}-\mathrm{J})^{n-j}\mathrm{Q}_{y}\mathrm{E}_{y}\mathrm{J}^{j}$. Taking the weight $\mathrm{W}_{\mathrm{Q}}$ into consideration the full Volterra integral operator is then

[TABLE]

This straightforward approach evidently only works if the kernel is of a form that may sensibly be approximated using monomials but it inspires an analogous approach based on expanding the kernel in its own orthogonal polynomial basis which need not be the same as those used to expand the function $f$. We use a variant of the Clenshaw algorithm introduced in section 2.3 to build the kernel in terms of the Jacobi operators. In principle one could compute $K(\mathrm{J}_{x},\mathrm{J}_{y})$ as a full multiplication operator acting on a triangle Jacobi coefficient vector using an operator-valued version of Clenshaw’s algoritm as discussed in [36]. This is not the most efficient way to approach this problem, however, as it would mean losing the diagonal $\mathrm{Q}_{y}\mathrm{E}_{y}$ since for such an operator the multiplication with $K(\mathrm{J}_{x},\mathrm{J}_{y})$ would need to happen between $\mathrm{Q}_{y}$ and $\mathrm{E}_{y}$. Nevertheless, we will briefly discuss how to generate this multiplication by $K(\mathrm{J}_{x},\mathrm{J}_{y})$ operator in order to see which modifications one can make to this approach in order to respect the symmetries of the triangle and end up with recursive basis generation similar to the monomial kernel expansion case.

The multiplication by $K(x,y)$ operator, which we label $\mathrm{M}_{K}$, can be written in an operator Clenshaw approach as (see [36, 33, 46]):

[TABLE]

where $\otimes$ denotes the Kronecker product and $\mathcal{L}$ is defined as

[TABLE]

As discussed for the Clenshaw evaluation method in section 2.3 this system requires preconditioning to become solvable via backward substitution. For this case the preconditioner is

[TABLE]

with the $B_{n}^{+}$ defined as in section 2.3. Using such an operator valued Clenshaw algorithm one can compute $\mathrm{M}_{K}$ and thus obtain $\mathrm{Q}_{y}K(\mathrm{J}_{x},\mathrm{J}_{y})\mathrm{E}_{y}$ via $\mathrm{Q}_{y}\mathrm{M}_{K}\mathrm{E}_{y}$. However, as discussed above, for our purposes of Volterra integral operators this is computationally wasteful and misses the chance to take advantage of the triangle symmetries which allow for $\mathrm{Q}_{y}\mathrm{E}_{y}$ to be directly computable and diagonal. So instead we replace the $\mathbf{K}_{\vartriangle}$ in (12) by $(\mathbf{K}_{\vartriangle}\otimes\mathrm{Q}_{y}\mathrm{E}_{y})$. The relations (10–11) then imply that all $\mathrm{J}_{x}$ operators may be replaced by a left multiplication with $\mathrm{J}$ and all $\mathrm{J}_{y}$ operators may be replaced by a right multiplication with $\mathrm{J}$ (respectively denoted by a $\diamond$ on the appropriate side). The system to solve thus becomes

[TABLE]

with

[TABLE]

After preconditioning as above, this allows the recursive and efficient computation of $\mathrm{Q}_{y}K(\mathrm{J}_{x},\mathrm{J}_{y})\mathrm{E}_{y}$ via an operator valued Clenshaw type algorithm while at the same time taking advantage of the diagonal nature of $\mathrm{Q}_{y}\mathrm{E}_{y}$. As in the monomial case, this approach has to be modified when integrating from [math] to $x$ instead of from [math] to $1-x$. In the [math] to $x$ case one needs to take the reflection into account, which ends up either replacing all the left multiplications with $\mathrm{J}$ by left multiplications with $(\mathbb{1}-\mathrm{J})$ for the same reasons as above, while the right multiplications corresponding to $y$ multiplication remain the same, or requiring that $K(1-x,y)$ be supplied to the algorithm. Finally, this operator still requires left multiplication with the basis dependent weight $\mathrm{W}_{\mathrm{Q}}$ to represent the full Volterra integral operator for this approach.

3.3. Numerical solutions to linear Volterra integral equations

The above described computational method for Volterra integrals has a natural extension to solving Volterra integral equations, which we describe in this section. Most generally a Volterra integral equation is any equation in which the unknown appears at least once as the integrand of a Volterra integral as defined in (1) above. One usually distinguishes between at least two types of Volterra integral equations which are labelled Volterra integral equations of the first and second kind respectively. The Volterra integral equation of the first kind we will be interested in takes the following form:

[TABLE]

where $u(x)$ is the unknown function to be solved for, $K(x,y)$ is a given kernel and $g(x)$ is a given function. Volterra integral equations of the second kind we will be interested in take the following form:

[TABLE]

where once again $u(x)$ is the unknown function and $K(x,y)$ and $g(x)$ are given. While this is not further explored in this paper, there are natural extensions of these methods for other linear Volterra-type integral equations such as the third kind equations discussed in [2, 3, 42].

Whenever we write $\mathrm{Q}_{y}K(\mathbb{1}-\mathrm{J}_{x},\mathrm{J}_{y})\mathrm{E}_{y}$ in the coming sections, we mean to imply that this operator is computed using the Clenshaw approach detailed in section 3.2.

3.3.1. Equations of the first kind

Extending the above methods for Volterra integrals to Volterra integral equations is straightforward, though one needs to be mindful of the appropriate reflections. Using the above notation conventions, one way to write the Volterra integral equation of the first kind is

[TABLE]

The notation $\bar{\mathbf{g}}$ is used to indicate that we are directly supplying the coefficients of the reflected $g(1-x)$ to save an unnecessary additional reflection step, as formally we are solving the equivalent

[TABLE]

All function coefficient vectors in this section are initially expanded in the $\tilde{\mathbf{P}}^{(1,0)}(x)$ basis. This method works in numerical experiments but deriving convergence properties for it proves to be difficult (as is usual for Volterra equations of the first kind). However, under the condition that we can expand the function $q(x)=\frac{g(1-x)}{1-x}$ instead of $g(1-x)$ in $\tilde{\mathbf{P}}^{(1,0)}(x)$, one can find convergence conditions (see section 5 for details). Note that solvability of the Volterra integral equation of the first kind implies that both $g$ and $q$ must vanish when the upper limit of integration vanishes. When using $\mathbf{q}$ to denote the coefficient vector of $q(x)=\frac{g(1-x)}{1-x}$ the method then becomes

[TABLE]

meaning that solving this type of equation for $u(x)$ is as simple as computing the coefficient vectors and operators (see the respective sections above for efficient ways to do so) and then solving a banded system of linear equations.

3.3.2. Equations of the second kind

Using the above-introduced weighted lowering operator $\mathrm{L}_{(1,0)}^{(0,0)}$ which shifts to the $\tilde{\mathbf{P}}^{(0,0)}(x)$ basis while multiplying with $(1-x)$, reflecting the result and then using a raising operator $\mathrm{S}_{(0,0)}^{(1,0)}$ to return to the $\tilde{\mathbf{P}}^{(1,0)}(x)$ basis we can write Volterra integral equations of the second kind as

[TABLE]

which can once again be solved for $u(x)$ using any linear system of equations solver. Reflecting without the lowering and raising operator is not possible (although there are alternative ways to use such operators to accomplish the same goal) as this would result in an inconsistency between the bases used for the two appearances of $\mathbf{u}$.

3.3.3. Different limits of integration

As mentioned above, a similar derivation leads to an analogous method for Volterra integral equations of the first and second kind with different limits of integration:

[TABLE]

This results in an identity operator replacing the reflection and conversion operators in the above solution methods and in fact makes these types of equations even more efficient to solve but limits of integration of this sort are seen less often in applications. In particular, the operator version of Volterra integral equations of the first kind with limits of integration [math] to $1-x$ is:

[TABLE]

where now $\mathbf{q}$ is the coefficient vector of $q(x)=\frac{g(x)}{1-x}$. Equations of the second kind with these limits of integration can be written as:

[TABLE]

We present an implementation of both options for the limits of integration in the next section.

4. Examples and applications

We present three sets of numerical examples to validate our implementation. The first set concerns itself with Volterra integral equations of the first kind, the second with Volterra integral equations of the second kind with kernels of varying oscillatory intensity and the third set discusses a singular Volterra integral equation stemming from a heat conduction problem with mixed boundary conditions. As oscillatory kernels require high orders of polynomials to approximate accurately and the method was not designed for singular kernels, the second and third set are designed to test the method’s stability.

The computations presented in this section have been performed with an implementation of the scheme in the Julia programming language [8] in the framework of ApproxFun.jl and MultivariateOrthogonalPolynomials.jl [34, 33, 43]. The coefficients of the solution have relative accuracy with standard floating point arithmetic, even as they decay below machine precision. Values for absolute errors presented in this section converge beyond the precision of 64-bit floating point numbers because of the rapid convergence of the method and the way ApproxFun.jl implements function approximation (cf. [34, 33, 43])—the only time beyond 64-bit floating point precision numbers (via "BigFloat") were used is in the analytic solutions used as comparisons, as otherwise the convergence of the error would be capped by the precision at which the analytic solution is evaluated.

4.1. Set 1: Volterra integral equations of the first kind

We investigate the numerical solution of the following two example Volterra integral equations of the first kind:

[TABLE]

[TABLE]

The analytic solution to the first equation can be found to be:

[TABLE]

We present the absolute error between the analytic and numerical solution for $u_{1}(x)$ using the orthogonal polynomial method introduced in this paper in Figure 1A for different matrix dimensions $n\times n$ and the absolute error between the numerical solution for $u_{2}(x)$ and a high degree solution computed with $n=5050$ in Figure 1B.

4.2. Set 2: Volterra integral equations of the second kind with oscillatory kernels

We seek numerical solutions $u_{1}$, $u_{2}$ and $u_{3}$ to the following three Volterra integral equations of the second kind with kernels of varying oscillatory intensity:

[TABLE]

Accurate approximation of these kernels on the canonical triangle domain requires coefficient vectors of length exceeding $10^{3}$. We include contour plots of the specified kernels on said domain in Figure 2. One can find an analytic solution to the first equation:

[TABLE]

For the other two equations, we instead compare to a numerical solution of high degree ($n=5050$). We plot the absolute error convergence of the numerical solutions in Figure 3. Due to the oscillatory character of these kernels and the number of coefficients involved, this can be considered a moderate stress test of the Clenshaw approach to the computations of the Volterra integral operator.

4.3. Set 3: Singular Volterra integral equation of the second kind in heat conduction with mixed boundary conditions

Finally we discuss a more application-oriented example discussed in a handful of different variations in [18, 17, 16, 48, 7]:

[TABLE]

To see how equations of this type can result from heat conduction problems of the form $\frac{\partial^{2}u}{\partial x}-\frac{1}{\alpha^{2}}\frac{\partial u}{\partial y}=0$ with mixed boundary conditions, see for example [17]. This equation varies both in its singularity properties as well as its number of solutions depending on the parameter $\mu$. This example equation stemming from an application of Volterra integrals demonstrates that the method developed in this paper has a broader range of applicability and can in some cases extend to certain classes of singular problems as well, despite this not being part of the considerations during the development of the method. For testing purposes we choose the following for $g(x)$:

[TABLE]

The following analytic solutions to these equations can be found for general $\mu$ for $g_{1}$ (e.g. in [48]) and for $\mu=3$ for $g_{2}$:

[TABLE]

As the kernel is separable, the problem can instead be treated as

[TABLE]

which can be solved by appropriately adding multiplications with Jacobi operators or altering the supplied $g(x)$ in the method to solve Volterra integral equations of the second kind. We plot numerical solutions obtained for $g_{1}(x)$ with $\mu=7$ and $g_{2}(x)$ with $\mu=3$ in Figure 4. The naturally more error prone neighborhood of the singularity can be well approximated arbitrarily close to the singularity (though not at the exact point of the singularity itself) using higher values of $n$ if needed. For $g_{2}(x)$ the method shows no instability at the weak singularity of the kernel.

5. Stability and convergence of the method

In this section we make use of the fact that the coefficient space of orthogonal polynomials is equivalent to an infinite-dimensional Banach space (in particular a sequence space). The strategy for the analysis of the method is to show that the operators to be inverted for Volterra integral equations of the second kind can be written as compact perturbations of the identity (compare [33, 41, 28]), i.e. can be written as

[TABLE]

where $\mathcal{K}$ is compact. Operators of this form are either invertible or neither injective nor surjective by the Fredholm alternative, cf. [6, 27]. The assumption of well-posedness for the equation thus guarantees that an operator of this form is invertible and standard convergence results for finite section methods [10] then guarantee convergence. We begin by discussing the solver for Volterra integral equations of the second kind, as the analysis for first kind problems is more involved.

5.1. Equations of the second kind

Definition 5.1**.**

We define the projection operators $\mathcal{P}_{n}:\ell^{2}\rightarrow\ell^{2}$ which map a given coefficient vector to a truncated version of itself with non-zero entries for the first $n$ coefficients only.

Definition 5.2**.**

The analysis operator $\mathcal{E}:L^{2}(0,1)\rightarrow\ell^{2}$ is the inclusion of a square integrable function into the $\ell^{2}$ coefficient space of the complete basis of orthogonal Jacobi polynomials, which is guaranteed to exist by the Stone–Weierstrass theorem and is a bounded operator. The synthesis operator is its inverse $\mathcal{E}^{-1}:\ell^{2}\rightarrow L^{2}(0,1)$, which is also bounded. Note the terms analysis and synthesis are terminology in frame theory [13, 14].

Lemma 5.1**.**

The coefficient space Volterra integral operator $\mathrm{V}_{K}$ is compact, where $\mathrm{V}_{K}:\ell^{2}\rightarrow\ell^{2}$ for a given kernel $K(x,y)\in L^{2}[T^{2}]$ with limits of integration [math] to $x$ acting on the coefficient vector Banach space $\ell^{2}$ of the Jacobi polynomials $\tilde{\mathbf{P}}^{(1,0)}(x)$ is of the form

[TABLE]

with the respective operators defined as in section 3.

Proof.

$\mathrm{V}_{K}=\mathrm{L}_{(1,0)}^{(0,0)}\mathrm{Q}_{y}K(\mathbb{1}-\mathrm{J}_{x},\mathrm{J}_{y})\mathrm{E}_{y}$ follows from the definition of the involved operators, see section 3. To see compactness of $\mathrm{V}_{K}$ we consider the following diagram of functions between Banach spaces which represents the formalized version of the method:

$\mathcal{V}_{K}$ for a kernel $K(x,y)\in L^{2}[T^{2}]$ is the Volterra integral operator for said kernel acting on $L^{2}(0,1)$. It is a classical result of functional analysis that such Volterra integral operators $\mathcal{V}_{K}$ are Hilbert–Schmidt operators and thus compact [31]. It follows that $\mathrm{V}_{K}=\mathcal{E}\circ\mathcal{V}_{K}\circ\mathcal{E}^{-1}$ is a finite composition of bounded and compact operators between Banach spaces and hence itself compact. ∎

Lemma 5.2**.**

For $\mathrm{V}_{K}$ and $\mathcal{P}_{n}$ defined as above, we have

[TABLE]

Proof.

This follows directly from the compactness of $\mathrm{V}_{K}$ and the fact that $\ell^{2}$ is a Hilbert space and thus has the approximation property [27]. ∎

The above lemma justifies referring to the finite-dimensional projections $\mathcal{P}_{n}\mathrm{V}_{K}\mathcal{P}_{n}^{\mathsf{T}}$ of the Volterra operator as approximations.

Lemma 5.3**.**

$\mathrm{S}_{(0,0)}^{(1,0)}\mathrm{R}\mathrm{L}_{(1,0)}^{(0,0)}\mathrm{Q}_{y}K(\mathbb{1}-\mathrm{J}_{x},\mathrm{J}_{y})\mathrm{E}_{y}$* is compact on $\ell^{2}$ and thus Volterra integral equations of the second kind can be written in the form $(\mathbb{1}+\mathcal{K})\mathbf{u}=\mathbf{g}$ with $\mathcal{K}$ compact.*

Proof.

The operators $\mathrm{S}_{(0,0)}^{(1,0)}$ and $\mathrm{R}$ acting on the Banach space $\ell^{2}$ can both readily be seen to be bounded operators from their definitions from the Jacobi polynomial’s recurrence relationships [32, 18.9.5]. The result then follows from the observation that the Volterra integral operator $\mathrm{L}_{(1,0)}^{(0,0)}\mathrm{Q}_{y}K(\mathbb{1}-\mathrm{J}_{x},\mathrm{J}_{y})\mathrm{E}_{y}$ was shown to be compact and composition of bounded operators with a compact operator yields a compact operator. ∎

An analogous chain of arguments immediately establishes:

Lemma 5.4**.**

The Volterra integral operator for the limits [math] to $1-x$ is compact and can be written as

[TABLE]

The method is thus also of the form in (24).

Corollary 5.5**.**

The method described in section 3.3 converges like $\|\mathbf{u}-\mathcal{P}_{n}\mathbf{u}\|\rightarrow 0$ as $n\rightarrow\infty$ for well-posed Volterra integral equations of the second kind.

Proof.

As the method is of the form in (24), i.e. $(\mathbb{1}+\mathcal{K})\mathbf{u}=\mathbf{g}$ with $\mathcal{K}$ compact, the result is a corollary of the above results combined with the known invertibility and convergence properties for problems of this form in finite section methods, see e.g. [10]. ∎

5.2. Equations of the first kind

The Fredholm alternative and Neumann series arguments underlying the proofs above break down for first kind problems as the Volterra operator $\mathrm{V}_{K}:\ell^{2}\rightarrow\ell^{2}$ is compact on the infinite dimensional Banach space $\ell^{2}$ and therefore is strictly singular, cf. [6]. Thus, while the finite dimensional approximations $\mathrm{V}_{n}$ of the Volterra operator may have an inverse $\mathrm{V}_{n}^{-1}$, it is not obvious that $\mathbf{u}_{n}=\mathrm{V}_{n}^{-1}\mathbf{q}$ converges to $\mathbf{u}$ in the limit. The problem can be made well-posed, however, if one considers the Volterra operator as a map between two different appropriately chosen Banach spaces. Under sufficient continuity assumptions as well as the assumption that a given Volterra integral equation of the first kind has a solution, this problem may then be salvaged by finding a preconditioner which allows us to rewrite it as a problem involving operators which are compact perturbations of Toeplitz operators. We begin by assuming a polynomial kernel from where an extension argument directly yields that it also applies for the non-polynomial case. Note that in this section we will prove convergence of the method only for the case of limits of integration [math] to $1-x$. This is not a limitation for the case of integral equations of the first kind, since solving

[TABLE]

and

[TABLE]

are formally equivalent, as solving one automatically solves the other with $t=1-x$. The reason for the particular choice for our proofs is that some arguments are more clear in this variant. Furthermore, as the monomial expansion and Clenshaw algorithm based Volterra operators are exactly the same for polynomial kernels the analysis will make use of the simpler structure of the former.

To discuss invertibility for equations of the first kind we need to reframe the Volterra operator as a a map between two different Banach spaces, which are similar in spirit to Sobolev spaces.

Definition 5.3**.**

Let $\ell^{2}_{\lambda}$ with $\lambda\geq 0$ denote the Banach space with norm

[TABLE]

Any $\mathbf{u}\in\ell^{2}_{\lambda}$ corresponds uniquely to a $\mathbf{u}\in\ell^{2}$ so we have $\ell^{2}_{\lambda}\subset\ell^{2}$ whereas the converse is clearly not the case.

Lemma 5.6**.**

Let $\mathrm{V}_{K}:\ell^{2}\rightarrow\ell^{2}_{1}$ denote the Volterra operator in coefficient space of $\tilde{\mathbf{P}}^{(1,0)}(x)$ with limits of integration [math] to $1-x$ for a given polynomial kernel

[TABLE]

Then

[TABLE]

with $\mathrm{D}=\mathrm{Q}_{y}\mathrm{E}_{y}$, $\mathrm{D}:\ell^{2}\rightarrow\ell^{2}_{1}$ and $\mathrm{D}^{-1}:\ell^{2}_{1}\rightarrow\ell^{2}$.

Proof.

That $\mathrm{D}=\mathrm{Q}_{y}\mathrm{E}_{y}$ is diagonal with entries $\frac{(-1)^{n+1}}{n}$ is due to properties of the Jacobi polynomials, see section 3 as well as [32, 18.6.1 and 18.17.1]. The important observation to make is that $\mathrm{D}$ can be thought of as $\mathrm{D}:\ell^{2}\rightarrow\ell^{2}_{1}$, which makes $\mathrm{D}$ a bounded and invertible operator with $\mathrm{D}^{-1}:\ell^{2}_{1}\rightarrow\ell^{2}$. With $\mathrm{V}_{K}$ and $K(x,y)$ as above, we thus have

[TABLE]

via Section 3.2. ∎

Definition 5.4**.**

When solving Volterra integral equations of the first kind with the method described in Section 3.3, it is useful to distinguish the operator without the weight $(1-x)$ which is to be inverted from the full Volterra operator. We will denote this operator $\tilde{\mathrm{V}}_{K}:\ell^{2}\rightarrow\ell^{2}_{1}$, where

[TABLE]

We furthermore see that

[TABLE]

as an immediate corollary of Lemma 5.6.

Lemma 5.7**.**

$\tilde{\mathrm{V}}_{K}$* may be written as*

[TABLE]

where $\mathrm{T}[f]$ is a Toeplitz operator with symbol $f$ and $\mathcal{K}$ is compact. Furthermore, the symbol is uniquely determined by the coefficients of the polynomial kernel $K(x,y)=\sum_{n=0}^{M}\sum_{j=0}^{n}k_{nj}x^{n-j}y^{j}$ to be

[TABLE]

Proof.

From the Lemma 5.6 we see that the first statement is equivalent to the claim that

[TABLE]

is of the form $\mathrm{T}+\mathcal{K}$ and thus asymptotically Toeplitz. To show this we need two observations: First, under sufficient continuity assumptions for the kernel, which are satisfied due to the kernel being polynomial, we have that

[TABLE]

and in particular

[TABLE]

where $H[a]$, $H[\bar{a}]$ and $H[\bar{b}]$ are compact Hankel operators [9]. Thus any asymptotically Toeplitz operator (of sufficiently continuous symbol) raised to a finite power is again an asymptotically Toeplitz operator, as $(T+\mathcal{K})^{2}=T^{2}+T\mathcal{K}+\mathcal{K}T+\mathcal{K}^{2}$ and $T^{2}$ is again Toeplitz plus something compact via the above relation. The composition of bounded operators with compact operators is compact making $T\mathcal{K}+\mathcal{K}T+\mathcal{K}^{2}$ compact. An induction argument demonstrates that this is true for any power $n\in\mathbb{N}$. In particular, since it is known that $\mathrm{J}$ is a compact perturbation of a Toeplitz operator [32] we know that $\mathrm{J}^{j}$ is a compact perturbation of a Toeplitz operator as well. The second observation is that for the banded operator $\mathrm{J}^{n-j}$, the operator $\mathrm{D}^{-1}\mathrm{J}^{n-j}\mathrm{D}$ is also a compact perturbation of a Toeplitz operator and in fact we have that $\mathrm{J}^{n-j}$ and $\mathrm{D}^{-1}\mathrm{J}^{n-j}\mathrm{D}$ differ only in their compact part, i.e. have the same Toeplitz component. Via (25) we thus have that $\sum_{n=0}^{M}\sum_{j=0}^{n}k_{nj}\mathrm{D}^{-1}\mathrm{J}^{n-j}\mathrm{D}\mathrm{J}^{j}$ is of the form $(\mathrm{T}+\mathcal{K})$ and thus asymptotically Toeplitz.

Along with the above observations, Equation (25) tells us that we can compute the symbol of the Toeplitz part of a product of operators which are compact perturbations of Toeplitz operators if we know the symbols of the individual Toeplitz components. Due to bandedness it is straightforward to confirm that the symbol of the Toeplitz part of the multiplication operator $\mathrm{J}$ is $(\frac{1}{2}+\frac{z}{4}+\frac{\bar{z}}{4})=\cos^{2}\left(\frac{\theta}{2}\right)$ for the Jacobi polynomials $\tilde{\mathbf{P}}^{(1,0)}(x)$, which is thus also the symbol of the Toeplitz part of $\mathrm{D}^{-1}\mathrm{J}\mathrm{D}$. Note at this point that

[TABLE]

due to the outer operators cancelling. Given these tools as well as the linearity of the Fourier series it follows that the symbol of the Toeplitz part of the Volterra operator $\tilde{\mathrm{V}}_{K}$ is the linear combination

[TABLE]

∎

Theorem 5.8**.**

The method described in Section 3.3 converges for well-posed Volterra integral equations of the first kind with limits of integration [math] to $1-x$

[TABLE]

rewritten as

[TABLE]

with $q(x)=\frac{g(x)}{1-x}$ for a polynomial kernel $K(x,y)\in L^{2}[T^{2}]$ and with $\mathbf{q}\in\ell^{2}_{1}$, subject to the symbol of the Toeplitz part of $\tilde{\mathrm{V}}_{K}$ not vanishing on the complex unit circle. This condition is fulfilled if and only if $\forall x\in[0,1]:K(x,x)\neq 0$.

Proof.

The requirement $\mathbf{q}\in\ell^{2}_{1}$ arises formally due to the need to first invert $\mathrm{D}$ and can be understood as stemming from the inverse integration being a differentiation. The invertibility conditions of asymptotically Toeplitz operators of the form $(\mathrm{T}+\mathcal{K})$ are known in the literature (see e.g. [22, 10] and the references therein): A compactly perturbed Toeplitz operator on $\ell^{2}$ is invertible if it is a Fredholm operator, its index is [math] and it has a trivial kernel [21, 10, 22]. Furthermore, a compactly perturbed Toeplitz operator is Fredholm if its symbol (which is just the symbol of the Toeplitz part) does not vanish anywhere on the complex unit circle.

In general, it holds that the index of a Toeplitz operator which is Fredholm is the sign-flipped winding number of its symbol on the complex unit disk [10]. Since the symbol of the Toeplitz part of the unweighted Volterra operator is real-valued and continuous its index is thus 0 if and only if it does not vanish anywhere on the complex unit circle, which is a necessary condition for it to be Fredholm in the first place. Since $\cos^{2}\left(\frac{\theta}{2}\right)\in[0,1]$, the symbol vanishes at some point $\theta\in[0,2\pi]$, i.e.

[TABLE]

if and only if for some $x\in[0,1]$ we have

[TABLE]

This in turn is precisely the condition that $K(x,x)=0$, since

[TABLE]

Conversely, if $\forall x\in[0,1]:K(x,x)\neq 0$ then the Volterra operator is Fredholm because the symbol of its Toeplitz part has no roots on the unit circle and as this symbol is real valued its winding number and thus index is 0. This necessary condition for invertibility of the operator becomes a sufficient condition if in addition to this we have $\text{ker}(\mathrm{T}+\mathcal{K})=\{0\}$, as this yields injectivity and via the index formula [10]:

[TABLE]

with $\text{ind}(\mathrm{T}+\mathcal{K})=0$ also implies surjectivity. $\text{ker}(\mathrm{T}+\mathcal{K})=\{0\}$ is a consequence of the classical result that the Volterra integral operator has no non-zero eigenvalues. The convergence of the method is then a consequence of known results in the theory of finite section methods, see e.g. [22]. ∎

Remark: The motivation for solving $\tilde{\mathrm{V}}_{K}\mathbf{u}=\mathbf{q}$ with $q(x)=\frac{g(x)}{1-x}$ instead of $\mathrm{V}_{K}\mathbf{u}=\mathbf{g}$ directly can be understood at this point, since for $\mathrm{V}_{K}$ the symbol of the Toeplitz part is instead found to be

[TABLE]

which always has a root on the complex unit circle at $\theta=0$ and thus its induced Toeplitz operator is not Fredholm and not invertible. Therefore the presented proof strategy only succeeds if $q(x)=\frac{g(x)}{1-x}$ may be used instead to get rid of the additional sine terms. The symbol of the Toeplitz part of $\tilde{\mathrm{V}}_{K}$ is comparably very well-behaved for a variety of kernels.

So far we have only been working with polynomial kernels of order $M$, henceforth denoted $K_{M}$, when it comes to Volterra equations of the first kind. We will need the following theorem (see [4, 44]) which we restate without proof for the extension of the above arguments to a non-polynomial kernel:

Theorem 5.9**.**

Let $X$ and $Y$ be normed linear spaces with one or both being Banach spaces and let $\mathcal{T}:X\rightarrow Y$ be a bounded and invertible operator with $\mathcal{T}^{-1}:Y\rightarrow X$. Then if the bounded operator $\mathcal{M}:X\rightarrow Y$ satisfies

[TABLE]

it follows that $\mathcal{M}$ is also invertible with bounded inverse operator $\mathcal{M}^{-1}:Y\rightarrow X$ and

[TABLE]

[TABLE]

Lemma 5.10**.**

Given that

[TABLE]

for a sequence of Volterra operators induced by polynomial kernels $K_{M}(x,y)$ and a not necessarily polynomial kernel $K(x,y)$, we have

[TABLE]

where $\mathbf{u}_{M}$ is the solution to the approximated problem

[TABLE]

Proof.

The method can be extended to more general $K=K(x,y)$ if $K_{M}$ is interpreted as the polynomial approximation of order $M$ of the full kernel $K$. To show that the method can be extended sensibly to non-polynomial kernels what remains to be shown is that $\|\mathbf{u}_{M}-\mathbf{u}\|\xrightarrow[M\rightarrow\infty]{}0.$ This can be achieved by use of Theorem 5.9: The assumptions of the theorem are satisfied when setting $\mathcal{T}=\tilde{\mathrm{V}}_{K}$ and $\mathcal{M}=\tilde{\mathrm{V}}_{K_{M}}$ since if $\|\tilde{\mathrm{V}}_{K_{M}}-\tilde{\mathrm{V}}_{K}\|\xrightarrow[M\rightarrow\infty]{}0$ then for some $M$ all subsequent $\tilde{\mathrm{V}}_{K_{M}}$ satisfy

[TABLE]

This immediately yields invertibility of $\tilde{\mathrm{V}}_{K_{M}}$ and more importantly the desired result that

[TABLE]

which justifies calling the solution $\mathbf{u}_{M}=\tilde{\mathrm{V}}_{K_{M}}^{-1}\mathbf{q}$ an approximation to $\mathbf{u}=\tilde{\mathrm{V}}_{K}^{-1}\mathbf{q}$. ∎

6. Discussion

The method proposed in this paper can efficiently compute Volterra integrals as well as solve Volterra integral equations of the first and second kind with high accuracy using bivariate orthogonal polynomials to resolve the kernel along with an operator valued Clenshaw algorithm and is not restricted to convolution kernels. Numerical experiments suggest it can even be applicable to certain singular equations. Our approach takes advantage of the sparsity of the required integration and extension operators which are due to the symmetries of the Jacobi polynomial basis on the triangle domain. The method was shown to converge for well-posed Volterra integral equations of the first and second kind, using a link to compact perturbations of Toeplitz operators.

Extensions of this approach to various so-called integro-differential equations of Volterra type, where both differentiation and Volterra operators act on the unknown function, as well as extensions to non-linear Volterra equations, where the unknown function can appear in non-linear fashion in the Volterra integral, while non-trivial are conceivable and will be addressed in future works.

Acknowledgments

We thank Nick Hale and Kuan Xu for crucial help on Volterra integral equations at the initial stages of this project. We thank Mikael Slevinsky for thoroughly reading a draft and providing detailed comments.