    # Computing with functions in the ball

A collection of algorithms in object-oriented MATLAB is described for numerically computing with smooth functions defined on the unit ball in the Chebfun software. Functions are numerically and adaptively resolved to essentially machine precision by using a three-dimensional analogue of the double Fourier sphere method to form "ballfun" objects. Operations such as function evaluation, differentiation, integration, fast rotation by an Euler angle, and a Helmholtz solver are designed. Our algorithms are particularly efficient for vector calculus operations, and we describe how to compute the poloidal-toroidal and Helmholtz–Hodge decomposition of a vector field defined on the ball.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Three-dimensional spherical geometries are common in computational science and engineering, arising in weather forecasting [coiffier2011fundamentals], geophysics [flyer2009radial, orszag1974fourier, spotz1998fast, wright2010hybrid], hydrodynamics [hollerbach1994magnetohydrodynamic, kuang1999numerical, zhang2000magnetohydrodynamics, zhang1989convection], and computational fluid dynamics [kerswell2005recent, serre2001three]. In each of these applications, it is routine to derive models that are continuous, even though one immediately discretizes them to compute a solution. Ballfun is a software system written in MATLAB that exploits object-oriented programming to allow users to compute with scalar- and vector-valued functions defined on the three-dimensional unit ball while being oblivious to our underlying discretizations. Ballfun is the first extension of Chebfun to three-dimensional spherical geometries [Driscoll2014] and follows the development of Spherefun [townsend2016computing] and Diskfun [wilber2017computing] for computing with functions in the sphere and the unit disk. Software systems in Dedalus [burns2019dedalus] (written in Python) and Approxfun [olver2019approxfun] (written in Julia) for computations on the ball may follow soon. Dedalus and Approxfun already have excellent functionality for computing on the 2-sphere and disks [olver2019approxfun, vasil2019tensor].

For computations with functions defined on the unit ball, a standard approach is to employ spherical coordinates , where , , and

denote the radial, azimuthal, and polar variables, respectively. Thus, computations on the unit ball can be conveniently related to analogous tasks involving functions defined on a cuboid, which allows for efficient algorithms based on tensor-product structure. Unfortunately, this simple coordinate transform comes with several significant disadvantages due to the artificial pole singularities introduced by the transform.

In this paper, we employ a technique known as the double Fourier sphere (DFS) method [fornberg1997comparison, merilees1973pseudospectral, orszag1974fourier] in conjunction with tensor-product expansions of functions. More precisely, we use a three-dimensional analogue of the DFS method that extends ideas from the disk and sphere [townsend2016computing, wilber2017computing] (implemented in Spherefun and Diskfun, which are part of Chebfun), while preserving the additional structure that is present in the 3D ball (see Definition 1). The DFS method alleviates some of the computational difficulties with spherical coordinates while having approximants that have an underlying tensor-product structure for efficient algorithms and FFT-based fast transforms. We use DFS approximants to develop a collection of algorithms for performing everyday computational tasks on scalar- and vector-valued functions defined on the unit ball and thus provide a convenient computational environment for scientific explorations.

Our algorithms are designed, whenever possible, to compute each operation on a function to essentially machine precision by using data-driven techniques from approximation theory. Our codes are also designed to have no required user-defined algorithmic parameters and to be as intuitive as possible for MATLAB users. For example, sum(v) returns the sum of the entries of a vector v in MATLAB, while sum3(f) returns the integral of a function f over the unit ball. Moreover, v.*w performs entry-by-entry vector multiplication, while f.*g returns a function representing the multiplication of f and g in Ballfun. During the operation f.*g, our algorithm automatically selects the discretization of the output so that the result is approximated to essentially machine precision. We repeat this idea in the one hundred or so Ballfun commands by constantly expanding and pruning underlying discretizations to represent functions as efficiently as possible.

There are several existing approaches for computing with functions on the unit ball and we seriously considered two other approaches:

[wide=]

Spherical harmonic expansions:

Spherical harmonic expansions of a function are given by , where is a spherical harmonic. These expansions can be thought of as the ball analogue of trigonometric expansions for periodic functions. When truncated, they provide essentially uniform resolution of a function over the ball. They have major applications in geophysics [lowrie2007fundamentals] and the numerical solution of separable elliptic equations.

Orthogonal polynomials on the ball:

Given an appropriate weight function on the ball, one can derive various families of orthogonal polynomials that are built from ultraspherical polynomials [dunkl2014orthogonal, Sec. 5.1]. Expanding functions in any one of these bases provides excellent resolution properties, along with fast evaluation, differentiation, and integration of the expansions. Unlike spherical harmonics expansions, they are rarely employed in practice.

We require a representation for functions on the ball that can be adaptively computed, as we would like to achieve an accuracy close to machine precision. While there are optimal-complexity spherical harmonic transforms [slevinsky2017fast], it is highly desirable to have the most computationally efficient fast transform associated with an expansion. The DFS method offers a simple and computationally efficient fast transform based on the FFT (see Section 2.2). Unlike spherical harmonic and orthogonal polynomial expansions, the DFS method does not guarantee that an expansion is infinitely differentiable on the ball, even when the original function is. For this reason, our algorithms must strictly preserve a structure in the DFS expansions to ensure that they represent a smooth function on the ball (see Definition 1).

Using DFS approximants, we develop a variety of algorithmic tools to provide a convenient computational environment for integrating, differentiating, and solving partial differential equations (see

Section 4), as well as representing vector-valued functions. This allows us to develop a set of algorithms for performing vector calculus (see Section 5), including computing the Helmholtz–Hodge and poloidal-toroidal decomposition.

The paper is organized as follows. We briefly introduce the software that accompanies this paper in Section 1.1. Then, in Section 2, we explain the methods used to discretize smooth functions on the ball. Next, in Section 3, we discuss some of the operations implemented in the software such as integration, differentiation, and a fast rotation algorithm. Following this, in Section 4, we describe a fast and spectrally accurate Helmholtz solver for solving equations with Dirichlet or Neumann boundary conditions. Finally, Section 5 consists of a description of the vector calculus algorithms, including the poloidal-toroidal and Helmholtz–Hodge decompositions.

### 1.1 Software

Ballfun is part of Chebfun [Driscoll2014], which is a software system for computing with functions and solving differential equations on an interval [battles2004extension], rectangle [townsend2013extension], cuboid [hashemi2017chebfun], disk [wilber2017computing], and the surface of a sphere [townsend2016computing]. Accompanying this paper is the publicly available MATLAB code in Chebfun [Driscoll2014] with two new classes called ballfun and ballfunv. We encourage the reader to explore this paper with the latest version of Chebfun downloaded and ready for interactive exploration. On the Chebfun website, we provide documentation in the form of a chapter of the Chebfun Guide [Driscoll2014] as well as several examples.111The Ballfun examples are available at http://www.chebfun.org/examples/sphere/. Functions on the ball can be easily constructed in the software by calling the appropriate command. For instance, f = ballfun(@(x,y,z) sin(cos(y))) defines the function . Underneath, Ballfun adaptively resolves the function to machine precision and represents it using the DFS method. For example,

  f = ballfun(@(x,y,z) sin(cos(y))) % ballfun representing sin(cos(y))
ballfun object:
Ψ    domain           r    lambda    theta
unit ball        21      45       41


where , , and are discretization parameters that Ballfun automatically determined necessary to resolve to machine precision. The Ballfun software is highly adaptive and automatically truncates the expansion to resolve functions on the ball to machine precision after each operation. After its construction (see Section 2), a function can be manipulated and analyzed through the nearly one hundred operations implemented in the package (see Table 1 and Table 2).

## 2 The Ballfun constructor

In this section, we explain how smooth functions, expressed in Cartesian or spherical coordinates, are discretized and constructed in our software. Smooth functions on the ball expressed in the spherical coordinate system can potentially introduce artificial boundaries at the origin or poles, as well as the loss of periodicity in the polar variable

. To overcome this issue, we first sample functions on a tensor-product grid in spherical coordinates. Then, we compute a Chebyshev–Fourier–Fourier expansion that interpolates the samples, using the ball analogue for the double Fourier sphere (DFS) method.

### 2.1 The double Fourier sphere method in the ball

The DFS method for the sphere was originally proposed by Merilees [merilees1973pseudospectral] and is used to construct spherefun objects in Chebfun. It naturally extends to the 3D settings and maps a function defined on a ball onto a 3D cuboid so that the origin and poles of the ball are not treated as artificial boundaries and the polar variable can be represented in a Fourier series [boyd1978choice, fornberg1995pseudospectral, orszag1974fourier, yee1980studies]. The method can also be applied to disks, cylinders, and ellipsoids [townsend2016computing, wilber2017computing].

The ball analogue of the DFS method is obtained by constructing a Chebyshev–Fourier–Fourier expansion of a function defined on instead of . A continuous function on the ball is first written in spherical coordinates as

 f(r,λ,θ)=f(rcosλsinθ,rcosλcosθ,rcosθ),(r,λ,θ)∈[0,1]×[−π,π]×[0,π].

The function is not periodic in . Under the DFS mapping, it is recovered by “doubling up” the polar variable to in the sense that is sampled twice. The radial variable is also doubled to remove the artificial boundary at . Using this ideas, we extend the function to a new function , defined on . The function can be expressed as

 ~f(r,λ,θ)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩g(r,λ+π,θ),(r,λ,θ)∈[0,1]×[−π,0]×[0,π],h(r,λ,θ),(r,λ,θ)∈[0,1]×[0,π]×[0,π],g(−r,λ,π−θ),(r,λ,θ)∈[−1,0]×[0,π]×[0,π],h(−r,λ+π,π−θ),(r,λ,θ)∈[−1,0]×[−π,0]×[0,π],h(r,λ+π,−θ),(r,λ,θ)∈[0,1]×[−π,0]×[−π,0],g(r,λ,−θ),(r,λ,θ)∈[0,1]×[0,π]×[−π,0],h(−r,λ,π+θ),(r,λ,θ)∈[−1,0]×[0,π]×[−π,0],g(−r,λ+π,π+θ),(r,λ,θ)∈[−1,0]×[−π,0]×[−π,0], (1)

where

 g(r,λ,θ)=f(r,λ−π,θ),h(r,λ,θ)=f(r,λ,θ),(r,λ,θ)∈[0,1]×[0,π]×[0,π].

Functions that satisfy Eq. 1 are said to be block-mirror-centrosymmetric (BMC) [wilber2017computing]. A more intuitive description is given by the visualization

 (2)

where flip1 (resp. flip3) refers to the MATLAB command that reverses the order of the first (resp. third) component of a tensor.

In addition to satisfying the BMC structure, must be constant at as well as and , corresponding to the origin and the poles. We call these functions BMC-III functions. (BMC-I and BMC-II functions are defined in [townsend2016computing] and [wilber2017computing], respectively.)

###### Definition 1 (BMC-III function)

A function is a BMC-III (Type-III BMC) function if it is a BMC function, , and, for any , and , where and only depend on such that , for some constant .

Figure 1 shows the DFS method applied to the earth and the type-III BMC structure.

There are two salient features of the DFS method that make it attractive for developing a package for computing with functions on the ball. First, tensor product expansions of Fourier and Chebyshev bases can be used to represent . If is a function in Cartesian coordinates on the ball, then after applying the DFS method, we have a function defined on the cuboid that can be approximated as

 ~f(r,λ,θ)≈n−1∑i=0n/2−1∑j=−n/2n/2−1∑k=−n/2αijkTi(r)eijλeikθ, (3)

where are spherical coordinates, denotes the Chebyshev polynomial of the first kind of degree , and is an even integer that is determined by the adaptive procedure described in Section 2.3.222For simplicity, we have chosen the same number of terms in each sum. Ballfun uses a generalization of Eq. 3 with a different number of terms in each sum. This representation allows us to use fast transforms as well as 1D and 2D algorithms for Chebyshev and trigonometric expansions. The second feature is that the DFS mapping of a function leads to a BMC structure (see Fig. 1) that, if preserved, ensures smoothness of the solution throughout the entire domain [townsend2016computing, wilber2017computing]. The BMC symmetry is imposed exactly by evaluating the function on and extending it to using Eq. 2. This ensures that the resulting function is smooth on the ball, i.e., at least continuous and differentiable. There are representations of functions on the ball that preserve full regularity [vasil2019tensor], however there are less appropriate in our settings since they do not allow efficient FFT-based transforms.

### 2.2 Computing the Chebyshev–Fourier–Fourier coefficients

Once the BMC-III function is found, it is approximated by a truncated Chebyshev–Fourier–Fourier (CFF) series [mason2002chebyshev, trefethen2000spectral, trefethen2013approximation]. For some even integer , the CFF coefficients are stored as an tensor and the entries are computed in operations, as follows:

1. [wide=]

2. The function is evaluated over at the tensor-product grid:

 ⎛⎜ ⎜⎝cos⎛⎜ ⎜⎝(n2−1−i)πn−1⎞⎟ ⎟⎠,2jπn,2kπn⎞⎟ ⎟⎠,0≤i≤n2−1,−n2≤j≤n2−1,0≤k≤n2. (4)
3. The samples of are doubled-up (see Eq. 2). This extends them to be samples of on at the tensor-product grid:

 (cos(iπn−1),2jπn,2kπn),0≤i≤n−1,−n2≤j≤n2−1,−n2≤k≤n2−1

without any additional evaluations of .

4. The CFF coefficients are computed using the discrete Chebyshev tranform (DCT) [gentleman1972implementing, gentleman1972implementing2, mason2002chebyshev] and the FFT [cooley1965algorithm].

There is also the inverse procedure, which evaluates at the grid in (4) in operations. This operation is particularly important in our plotting commands and is achieved by reversing steps 1-3, using the inverse DCT and FFT.

### 2.3 Determination of the discretization size

To construct a ballfun object to represent a given function , we first sample at a CFF grid and compute the corresponding CFF coefficients (see Section 2.2). We then successively increase the grid size independently in each variable from to to , and so on, until we deem the function to be resolved in each variable. We use these samples to compute the CFF coefficients corresponding to an , where and then gauge the resolution in each variable by creating vectors of the absolute maximum of the coefficients along each variable, i.e.,

 Colsi=maxj,k∣∣αijk∣∣,Rowsj=maxi,k∣∣αijk∣∣,Tubesk=maxi,j∣∣αijk∣∣.

One can now inspect these vectors to identify whether or not the function is resolved to machine precision in each variable, relative to the magnitude of on  [aurentz2017chopping, hashemi2017chebfun, wright2015extension]. One can identify a near-optimal discretization size in that variable by recording the last entry in each vector above machine precision [aurentz2017chopping].

The constructor typically terminates when it encounters vectors Cols, Rows, and Tubes as shown in Fig. 2 for . In particular, for , the Ballfun constructor selected a CFF series of size to represent to essentially machine precision over the ball. Once the function is represented in a CFF expansion, the approximant is stored as a ballfun object, ready for further computation. For the rest of this paper, we will assume that the functions are represented by an tensor (though our software can deal with rectangular discretizations).

## 3 Algorithms for numerical computations with functions on the ball

Once a ballfun object is computed, there are many operations that can be performed on it. In fact, many of the operations can be decomposed into a sequence of 1D operations, which are particularly efficient for approximants of the form Eq. 3. This includes pointwise evaluation (see Section 3.1), integration (see Section 3.2), differentiation (see Section 3.3), and fast rotation (see Section 3.4).

### 3.1 Pointwise evaluation

The evaluation of a function at a point in the ball can be computed in operations. It follows from the CFF approximation of in Eq. 3 that one has

 ~f(r∗,λ∗,θ∗) =n−1∑i=0n/2−1∑j=−n/2n/2−1∑k=−n/2αijkTi(r∗)eijλ∗eikθ∗ =n/2−1∑k=−n/2⎛⎝n/2−1∑j=−n/2(n−1∑i=0αijkTi(r∗))eijλ∗⎞⎠eikθ∗.

Therefore, can be computed by first evaluating using Clenshaw’s algorithm [trefethen2013approximation, Chap. 19], which returns an matrix of values. Then, one can compute the summand over the index using Horner’s scheme [wright2015extension], which returns an vector, before finally computing the summand over the index using Horner’s scheme. This algorithm is implemented in the feval command and returns a scalar for . It is also possible to evaluate functions in Cartesian coordinates and Ballfun does a change of variables to spherical coordinates in this case.

### 3.2 Integration

The triple definite integral of a function on the unit ball can be written as follows:

 ∫B(0,1)~f(r,λ,θ)dV =∫10∫π0∫π−π~f(r,λ,θ)r2sinθdλdθdr =n−1∑i=0n/2−1∑j=−n/2n/2−1∑k=−n/2αijk∫10r2Ti(r)dr∫π−πeijλdλ∫π0sinθeikθdθ, =2πn−1∑i=0n/2−1∑k=−n/2αi0k(∫10r2Ti(r)dr)(∫π0sinθeikθdθ) =2πn−1∑i=0n/2−1∑k=−n/2αi0kνi1+eiπk1−k2,νi={2i2−6i4−10i2+9,i even,0,i.

Here, the last equality follows by calculating the integrals in explicitly (see, for example, [townsend2016computing, (4.3)]). Therefore, the integral of reduces to a basic task in linear algebra and can be computed in operations. This is implemented in Ballfun in the sum3 command. For example, the function has an integral of over the ball and can be computed in Ballfun by

  f = ballfun(@(x,y,z) x.^2); % ballfun representing x^2
sum3(f)                     % Integrate over the ball
ans =
0.837758040957278


The error is computed as abs(sum3(f) - 4*pi/15) and is given by .

### 3.3 Differentiation

Differentiation of a function on the ball with respect to spherical coordinates in , , and may introduce singularities at the poles and origin. For instance, consider the smooth function . The derivative of with respect to is

, which is not smooth at the poles. However, we are interested in computing derivatives that arise in vector calculus, such as the gradient, the divergence, the curl, or the Laplacian. All these operations can be expressed as partial derivatives in the Cartesian coordinates system. Therefore, our default is to allow for ballfun objects to be differentiated in the Cartesian coordinate system.

We follow the same approach as Spherefun [townsend2016computing] and express the partial derivatives in , , and in terms of the spherical coordinates , , and as follows:

 ∂∂x =cosλsinθ∂∂r−sinλrsinθ∂∂λ+cosλcosθr∂∂θ, (5) ∂∂y =sinλsinθ∂∂r+cosλrsinθ∂∂λ+sinλcosθr∂∂θ, (6) ∂∂z =cosθ∂∂r+sinθr∂∂θ. (7)

Then, Eq. 5, Eq. 6, and Eq. 7 involve operations on the tensor of CFF coefficients representing . For example, the derivative of with respect to can be expressed as

 ∂~f∂λ=n−1∑i=0n/2−1∑j=−n/2n/2−1∑k=−n/2αijkijTi(r)eijλeikθ.

Multiplications and divisions by , , , and in Eqs. 7, 6 and 5 are computed by multiplying the tensor of CFF coefficients by the corresponding matrices of linear operators, expressed in the Fourier basis. For example, we write , where satisfies

 B(:,j,:)=A(:,j,:)M−⊤sin,−n2≤j≤n2−1,Msin=i2⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣01−1⋱⋱⋱⋱1−10⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

Here, is the matrix of multiplication by in the Fourier basis. It is nonsingular if we choose

to be even (in this case the eigenvalues are

, ). Moreover, we write , where satisfies

 B(:,j,:)=M−1rA(:,j,:),−n2≤j≤n2−1,Mr= ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣012101212⋱⋱⋱⋱1212012120\par⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

Here, stands for the matrix of multiplication by in the Chebyshev basis. This matrix is invertible for even since its determinant is equal to .

Working directly on coefficients allows us to circumvent potential singularity issues at and the poles, while the standard technique uses a “shifted grid” procedure in the physical space [cheong2000application, fornberg1995pseudospectral, heinrichs2004spectral]. This procedure shifts the grid of sampled points in the latitude and radial directions, which avoids evaluation at the poles and the origin but can be numerically inaccurate near these points.

These operations are implemented in Ballfun in the diff command. For example, the derivative of with respect to is also represented as a ballfun object and can be computed as

  f = ballfun(@(x,y,z) cos(x.*y));% ballfun representing cos(xy)
diff(f, 1)                      % Compute ballfun representing df/dx
ans =
ballfun object
domain          r    lambda    theta
unit ball         24        43       40


Ballfun calls the constructor after each operation to readjust the grid sizes (see Section 2.3). Here, a discretization size of was determined necessary to resolve while is represented by a CFF series.

### 3.4 Fast rotation algorithm using a nonuniform Fourier transform

Rotating functions defined on the ball has applications in many fields, including quantum mechanics, inverse scattering, and geophysics. Ballfun has a rotate command to efficiently perform rigid-body rotations of functions. Every rigid-body rotation can be specified by an Euler angle in the Z-X-Z convention [arfken1985mathematical], which corresponds to rotating first by around the -axis then rotating by around the (original) -axis and then, finally, rotating by around the new -axis. All the angles are given in radians. The algorithm to achieve this rotation requires a nontrivial computation because the rotated function must be represented by an approximant in the original coordinate system.

The classical algorithm for computing the rotation of a function on the ball is to first express in terms of a spherical harmonic expansion and then to use the fact that the spherical harmonics form a basis of  [gelfand2018representations]. Since Ballfun does not represent functions using spherical harmonic expansions, we use an algorithm based on the DFS method and the 2D nonuniform FFT [ruiz2018nonuniform]. We do this by taking the CFF grid in (4) and rotating it by Euler angle . Then, we evaluate the function at this rotated grid and call the Ballfun constructor. Since the rotated grid is almost always non-uniform in the and variables of the doubled-up spherical coordinates and a Chebyshev grid in , the evaluation is done in operations with a 2D nonuniform FFT in and and a DCT in .

For example, the following code snippet calculates the rotation of by Euler angle (see Fig. 3).

f = ballfun(@(x,y,z) sin(50*z) - x.^2) % ballfun for sin(50z)-x^2
f =
ballfun object
domain          r    lambda    theta
unit ball       90         5        179
g = rotate(f, -pi/4, pi/2, pi/8)       % Rotate by (-pi/4,pi/2,pi/8)
g =
ballfun object
domain          r    lambda    theta
unit ball       91      180       182


As one can see, the rotate command is also adaptive and selects the appropriate discretization to resolve the rotated function.

## 4 Fast spectral method for solving Helmholtz equation

In this section, we describe a fast algorithm for solving the Helmholtz equation on the ball with Neumann boundary conditions. An optimal-complexity algorithm for solving Helmholtz equation with Dirichlet conditions on the ball is described in [fortunato2017fast], though it cannot immediately be generalized to the situation with Neumann conditions. Helmholtz solvers are useful in computational fluid dynamics as well as the computation of vector decompositions such as the poloidal-toroidal and Helmholtz–Hodge decompositions [backus1986poloidal, bhatia2013helmholtz].

### 4.1 Discretization of Helmholtz equation

Consider the Helmholtz equation on the ball, i.e., with Neumann boundary conditions on the sphere and a real wave number . The change of variables given by , where , , and , transforms the equation into

 1r2∂∂r(r2∂u∂r)+1r2sinθ∂∂θ(sinθ∂u∂θ)+1r2sin2θ∂2u∂λ2+K2u=f. (8)

One can multiply Eq. 8 by to remove the singularities in the variable coefficients at the origin and at the poles of the ball. We then use the DFS method (see Section 2.1) to represent over the domain . This allows us to solve Eq. 8 by seeking a tensor of Chebyshev–Fourier–Fourier coefficients for (see Eq. 3).

Let and be tensors of CFF coefficients of and , respectively. Since Eq. 8 decouples in the azimuthal variable , the following equation holds for :

 sin2θ∂∂r(r2∂uj∂r)+sinθ∂∂θ(sinθ∂uj∂θ)+(K2r2sin2θ−j2)uj=fj, (9)

where the functions and are defined by

We discretize Eq. 9 in the radial variable using the ultraspherical spectral method [olver2013fast], and in the polar variable using the Fourier spectral method. Partial derivatives in the polar variable and multiplication by are represented by sparse and banded matrices in the Fourier basis. The ultraspherical spectral method results in sparse and banded matrices of operators (such as differentiation or multiplication by ) between Chebyshev and ultraspherical polynomials. This allows us to write Eq. 9 in the form of a generalized Sylvester equation [gardiner1992solution] in the unknown matrix :

 LrU(:,j,:)M⊤sin2+S02U(:,j,:)L⊤θ=S02F(:,j,:), (10)

where is the matrix representing the operator from the Chebyshev basis to the ultraspherical basis and is the conversion matrix between these bases [olver2013fast]. The matrices and represent the multiplication by and the operator in the Fourier basis, respectively.

### 4.2 Imposing Neumann boundary conditions when K≠0

It is essential to slightly modify Eq. 10 to impose Neumann boundary conditions on , i.e., . The first step is to double-up the smooth function in the variable using the DFS method [townsend2016computing] and define its Fourier–Fourier matrix of coefficients . Since the radial variable of Eq. 8 has been doubled-up, we need to impose a Neumann condition at and . The matrix of coefficients of the boundary condition at , , can be deduced from (see Section 2.1) and takes the form

 G−(j,:)=(−1)jG+(j,:),−n2≤j≤n2−1.

The Neumann operators and are represented in the Chebyshev basis by the matrices

 B+(k)=k2,B−(k)=(−1)k+1k2,0≤k≤m−1,

respectively. The Neumann conditions also decouple in the variable and can be written as

 (B+B−)U(:,j,:)=(G+(j,:)G−(j,:)),−n2≤j≤n2. (11)

Therefore, a Helmholtz’s solver Eq. 8 with Neumann boundary conditions is realized by solving the following generalized Sylvester equation with linear constraints:

 LrU(:,j,:)M⊤sin2+S02U(:,j,:)L⊤θ=S02F(:,j,:), (12) (B+B−)U(:,j,:)=(G+(j,:)G−(j,:)), (13)

where . The constraints Eq. 13

can be used to remove degrees of freedom in

and transform Eq. 12 into a generalized Sylvester equation with a unique solution without constraints [townsend2015automatic], i.e.,

 ~LrXjM⊤sin2+~S02XjL⊤θ=~S02~F(:,j,:),−n2≤j≤n2−1. (14)

Fig. 4 shows the sparsity structure of the matrices and in Eq. 14. We obtain decoupled Sylvester matrix equations, where each one can be solved in operations using the Bartels–Stewart algorithm [bartels1972solution, gardiner1992solution]. Once has been computed, the matrix of coefficients can be recovered using the linear constraints. Thus, the total complexity is operations.

### 4.3 Imposing Neumann boundary conditions when K=0

We consider the zeroth Fourier mode of Eq. 10 with (Poisson equation). The solution to this equation with Neumann boundary conditions is unique only up to a constant. However, this additional constraint cannot be imposed on a Sylvester matrix equation. Therefore, we transform Eq. 10 into the Chebyshev–Legendre basis to decouple this Sylvester equation in the polar variable . The function , defined in Section 4.1, satisfies the following equation:

 ∂∂r(r2∂u0∂r)+1sinθ∂∂θ(sinθ∂u0∂θ)=r2f0. (15)

The functions and can be expressed in the Chebyshev–Legendre (CL) basis as

 u0(r,θ)=m−1∑i=0p−1∑k=0~ui0kTi(r)Pk(cosθ),r2f0(r,θ)=m−1∑i=0p−1∑k=0~fi0kTi(r)Pk(cosθ).

The zeroth Fourier modes of the Neumann boundary conditions at , , and at , , can also be written as Legendre series

 g+0(θ)=p−1∑k=0~g+0kPk(cosθ),g−0(θ)=p−1∑k=0~g−0kPk(cosθ).

The orthogonality of the Legendre basis allows us to decouple Eq. 15 in the polar variable as ordinary differential equations with Neumann boundary conditions:

 ∂∂r(r2∂~u0k∂r)−k(k+1)~u0k =~f0k, (16) ∂~u0k∂r∣∣∣r=1=~g+0k,∂~u0k∂r∣∣∣r=−1 =~g−0k, (17)

for . The functions and are defined by

 ~u0k(r)=m−1∑i=0~ui0kTi(r),~f0k(r)=m−1∑i=0~fi0kTi(r),r∈[−1,1].

For each Eqs. 17 and 16 can be solved in operations using the ultraspherical spectral method [olver2013fast]. The case is solved by the same technique with the additional linear constraint to impose uniqueness of the global solution .

Once the matrix of Chebyshev–Legendre coefficients has been computed, it can be converted to the Chebyshev–Fourier basis in operations using the Legendre–Chebyshev tranform [potts1998fast].

### 4.4 Numerical examples

In Fig. 5 (a) we plot a solution to the Helmholtz equation

 ∇2u+20u=−80sin(10x), (18)

with Neumann boundary conditions . The error between the computed and the exact solution to Eq. 18 is shown in Fig. 5 (c) and confirms the spectral convergence of our method. The computed solution is resolved to machine precision for . Our Helmholtz solver on the ball can be invoked in Ballfun via the helmholtz command.

  rhs = ballfun(@(x,y,z)-80*sin(10*x));    % Right-hand side
bc = @(x,y,z)10*x.*cos(10*x);            % Boundary conditions
K = sqrt(20);                            % Wave number
n = 50;                                  % Spectral discretization
u = helmholtz(rhs, K, bc, n, ’neumann’); % Helmholtz solver


The execution times333Timings were performed on a 3.3 GHz Intel Core i5 CPU using MATLAB 2018a without explicit parallelization. to solve Eq. 18 for different discretization sizes (see Fig. 5 (d)). We can then solve Helmholtz’s equation on a ball with one million degrees of freedom in a few seconds on a standard CPU.

As a second example we consider the advection-diffusion equation on the unit ball

 ∂c∂t=D∇2c−v⋅∇c, (19)

where is the diffusion coefficient and is a divergence-free vector field. We choose and to satisfy the no-slip condition. The no-flux condition for reduces to at the boundary. We impose the initial condition and solve Eq. 19 by using the implicit-explicit backward differentiation of order one (IMEX-BDF1) scheme. This yields the following Helmholtz equation:

 ∇2cn+1+K2cn+1=K2cn+1Dv⋅∇cn,∂c∂→n∣∣∣∂B(0,1)=0,

where denotes the solution at time , is the time step, and . The solution to Eq. 19 at different times is illustrated in Fig. 6.

## 5 Vector-valued functions on the ball

Ballfun is also designed to work with vector-valued functions defined on the unit ball as well as scalar-valued ones. Expressing vector-valued functions in spherical coordinates is inconvenient since the unit vectors in this coordinate system are singular at the poles of the ball [swarztrauber1981approximation]. Therefore, we express vector-valued functions in Cartesian coordinates as the components of the vector field are then themselves smooth functions. After using this convention, vector-valued functions introduce few complications from the point-of-view of approximation as each component is represented as an independent scalar function. A vector-valued function can be constructed in Ballfun as follows:

V = ballfunv(@(x,y,z) sin(x), @(x,y,z) x.*y, @(x,y,z) cos(z))
ballfunv object containing
ballfun object:
domain          r    lambda    theta
unit ball      14      27       27
ballfun object:
domain          r    lambda    theta
unit ball       3       5        5
ballfun object:
domain          r    lambda    theta
unit ball      15       1       29


It can be seen that a vector field is stored as a ballfunv object, which consists of three ballfun objects corresponding to the three components in Cartesian coordinates. Each ballfun object has its own discretization in , and .

### 5.1 Vector calculus on the ball

The more interesting side of vector-valued functions in Ballfun is the set of operations that can be implemented, which are potentially useful for applications. The standard operations for vector calculus such as the curl, the gradient or the divergence are implemented in Ballfun in the curl, gradient, and divergence commands, respectively. Due to the way we represent vector-valued functions, we compute these operations in Cartesian coordinate system. For example, the curl of a vector-valued function can be written as

 ∇×V=[∂Vz∂y−∂Vy∂z,∂Vx∂z−∂Vz∂x,∂Vy∂x−∂Vx∂y]T.

The curl of can be computed and displayed using the quiver command (see Fig. 7):

  W = curl( V ); % Compute curl of V
quiver( W )    % Plot vector field W using quiver


One can also verify basic vector calculus identities. For example, the divergence theorem asserts that the vector-valued function satisfies

 ∭B(0,1)(∇⋅V)dV=\oiint∂BV⋅→ndS,

where denotes the unit normal vector to . This theorem can be verified by executing the following code:

  lhs = sum3(divergence(V));  % Compute volume-integral of div( V )
nhat = spherefunv.unormal;  % Unit normal vector to surface of ball
Vbc = V(1,:,:, ’spherical’);% Restrict V to the bdy of ball
rhs = sum2(dot(Vbc, nhat)); % Compute dot-product & surface-integral


The error is .

### 5.2 Poloidal-toroidal decomposition

The poloidal-toroidal (PT) decomposition of a smooth divergence-free vector field expresses the field as the sum of two orthogonal fields. The PT decomposition is a well-known tool in fluid dynamics [horn2015toroidal] and magnetohydrodynamics [benton1983rapid, boronski2007poloidal, boronski2007poloidal2] simulations to analytically impose an incompressibility condition on flows in cylindrical and spherical geometries. In this section, we describe an algorithm for computing the PT decomposition of a smooth vector field in the ball.

Given a smooth divergence-free vector field, , defined on the ball, the PT decomposition [backus1986poloidal] writes as an orthogonal sum of a poloidal and toroidal field, i.e.,

 V=P+T,∫B(0,1)P⋅TdV=0.

Here, there exist two poloidal and toroidal scalar-valued functions and such that and , where is the unit radial vector. It can be shown that and are unique up to the addition of an arbitrary function of  [backus1986poloidal].

A vector field whose components are expressed in the Cartesian coordinate system can be converted in spherical coordinates using the following identities:

 Vr =sinθ(cosλVx+sinλVy)+cosθVz, Vλ =−sinλVx+cosλVy, Vθ =cosθ(cosλVx+sinλVy)−sinθVz.

Then, any poloidal and toroidal scalars for satisfy the following relations [backus1986poloidal]:

 ∇21Φ =−rVr, (20) ∇21Ψ <