Computer-assisted proofs for radially symmetric solutions of PDEs

We obtain radially symmetric solutions of some nonlinear (geometric) partial differential equations via a rigorous computer-assisted method. We introduce all main ideas through examples, accessible to non-experts. The proofs are obtained by solving for the coefficients of the Taylor series of the solutions in a Banach space of geometrically decaying sequences. The tool that allows us to advance from numerical simulations to mathematical proofs is the Banach contraction theorem.


1.
Introduction. In the last decade powerful techniques for turning numerical simulations of ordinary and partial differential equations into mathematically rigorous statements have developed rapidly (see e.g. [10,16,28,31,34,36]). These methods bridge the divide between scientific computing and abstract mathematical analysis. The foundations for the field were laid in [12,13,17]. These methods have received renewed attention due to contemporary computational power as well as recent algorithmic advances. With a clever hybrid approach one can off-load the verification of intricate computationally intensive estimates to the computer to prove existence to infinite dimensional continuum problems near numerical approximations. This permits using classical ideas (such as fixed point methods) to generate results not accessible by any other means. Thus, going beyond simulations, we can get our hands on solutions of nonlinear (geometric) partial differential equations (PDEs) using the computer and still argue about them with all the rigour of a classical pencil and paper proof.
In this paper we present examples of this technique which we believe are valuable because we prove existence of particular solutions of nontrivial nonlinear PDEs, while the technicalities are minimal (and the associated coding is easily manageable). In particular, we consider two instances of the nonlinear elliptic problem ∆u + f (u) = 0. (1) In the first type we consider u = u(x) to be scalar with the independent variable x lying in the 2-sphere. In the second type, the unknown u = u(x) is a vector in R 2 and x lies in the unit disc in R 2 (with Dirichlet boundary conditions). The specific choices for the nonlinearities f will be made explicit in Sections 3 and 4. When restricting attention to solutions which are invariant under rotation (radially symmetric solutions), these PDEs reduce to ordinary differential equations (ODEs). The boundary value problems for these nonlinear ODEs are amenable to an approach based on power series. Numerically one can obtain approximations for finitely many of the Taylor coefficients. To prove that such a truncated series corresponds to a solution of the differential equation, we will apply a Newton-Kantorovich type argument, see Section 2. This not only proves the existence of a solution, but also yields a rigorous bound on the truncation and approximation errors.
Finally, one ingredient that is needed in the analysis in order to control rounding errors, is computation using interval arithmetic. We refer to [21,24,29] for an introduction to interval arithmetic. The results in this paper make use of the package Intlab [23] for MATLAB. All MATLAB scripts are available at [4].
This paper, and in particular the examples discussed, resulted from a three week summer school for undergraduate students at Simon Fraser University in the summer of 2015. During that school, there were exploration sessions for students to investigate the ideas of rigorous computing through exercises and research problems. A small fraction of the students managed to obtain novel results about open problems, and these results constitute the content of the present paper.

2.
A Newton-Kantorovich type theorem. When trying to solve numerically a finite dimensional nonlinear zero-finding for F : R n → R n , the classical approach is to apply Newton's algorithm, where one iterates the map until the residual is sufficiently small. Often this numerical outcome x num is either assumed "good enough" (roughly, in applied mathematics) or "not rigorous" (roughly, in pure mathematics). Our goal is to merge these perspectives by providing a rigorous quantitative statement on how good this numerical outcome is, i.e., how close it is to a solution of F (x) = 0. The crucial observation is that Newton's method works so well because the map T is a contraction with a very small contraction constant on a neighborhood of a zero x sol of F , provided that the Jacobian DF (x sol ) is invertible, i.e., x sol is a non-degenerate zero of F (the generic situation). This allows us to prove that the numerical approximation x num is close to the true solution x sol .
Rather than going into detail about how to obtain such a computer-assisted proof, we first move to the infinite dimensional setting. We first construct the framework in general Banach spaces in Theorem 2.1 then, in Remark 1 and Example 2.2 present the (relatively straightforward) application to finite dimensional problems.
Consider now F : X → X a smooth map between two Banach spaces. A standard approach of mathematical analysis is to turn the zero finding problem F (x) = 0 into a fixed point problem. The Newton operator itself is usually impractical because the inverse of the derivative of an infinite dimensional map is hard to work with. Instead, one may choose an injective linear map A ∈ L(X , X) and study the fixed point problem The main difficulty is choosing A such that T is a contraction on some neighborhood of the (unknown) fixed point that we are looking for. In the context of a computerassisted proof, the approach is the following. First, using your favorite method from scientific computing, choose a finite dimensional 'projection' F num of F , solve it numerically to find x num , and reinterpret x num as an element of the infinite dimensional space, which we denote byx (hence F num (x num ) ≈ 0 and F (x) ≈ 0). We expect the solution to be close tox, hence we would like to choose A an approximate inverse of DF (x).
How can we determine that A is an appropriately accurate approximate inverse? And, on which neighborhood ofx do we have that T is a contraction? The following theorem is one way to make this precise. It uses, as an intermediate tool, and approximation A † ∈ L(X, X ) of DF (x).
Theorem 2.1. [Radii polynomial approach in Banach spaces] Let X and X be Banach spaces. Denote the norm on X by · X . Consider bounded linear operators A † ∈ L(X, X ) and A ∈ L(X , X). Assume F : X → X is C 1 , that A is injective and that AF : X → X. (2) Consider an approximate solutionx of F (x) = 0 (usually obtained using Newton's method on a finite dimensional projection). Let Y 0 , Z 0 , Z 1 , and Z 2 be positive constants satisfying where · B(X) is the induced operator norm for bounded linear operators from X to itself. Define the radii polynomial If there exists r 0 > 0 such that p(r 0 ) < 0, Proof. Using the mean value theorem, it is not hard to show that T maps B r0 (x) into itself, and that T is a contraction on that ball with contraction constant κ ≤ Z 0 + Z 1 + Z 2 r 0 < 1. The result then follows directly from the Banach contraction theorem. We refer to [15] for additional details.
This theorem fits in a long tradition of quantitative Newton-Kantorovich type theorems. The bounds are parametrized in terms of r, so that an appropriate r need not be assumed in advance but rather can be determined as the final step of the proof process. Moreover, we obtain balls for an interval of values for the radius. Small values of r 0 give us the tightest control on the distance between the solution and the numerical approximation, while large values of r 0 provide us with the best information about the isolation of the solution.
Before we show, in Sections 3 and 4 how to apply Theorem 2.1 in practice in infinite dimensional settings, we first consider to finite dimensional problems. Remark 1. In finite dimensions Theorem 2.1 can be implemented very easily and generally, if one has an interval arithmetic package, and if one has explicit formulas for D j F i and D 2 jk F i for 1 ≤ i, j, k ≤ n. We use the ∞-norm |x| = max 1≤i≤n |x i |. Let x = x num be the numerical approximation of a solution (e.g. found using Newton iterations). Let A be a numerical computed (hence approximate) inverse of the numerically computed DF (x num ). Now compute with interval arithmetic A † = DF (x num ), which implies we may set Z 1 = 0. Then evaluate, again with interval arithmetic, the following three computable expressions: where absolute values are taken component-wise, and sup denotes the supremum of the interval obtained. 2. the matrix norm By checking that Z 0 < 1, one verifies the hypothesis in Theorem 2.1 requiring A be injective.
3. the second derivative estimate (which provides a bound (6) via the mean value theorem) Here we choose a loose a priori upper bound r * on the value of r, and we bound the second derivative uniformly on this ball of radius r * .
Since Y 0 and Z 0 are usually near machine precision, the quadratic formula then gives a very small r 0 for which p(r 0 ) < 0. After checking that r 0 ≤ r * we can then invoke Theorem 2.1 to prove that a (unique) zero x sol of F lies within distance r 0 of x num . Example 2.2. We consider the circular restricted four body problem (CR4BP), where three bodies (with masses (m 1 , m 2 , m 3 ), normalized so that m 1 +m 2 +m 3 = 1 and m 1 ≥ m 2 ≥ m 3 ) move in circular periodic orbits around their center of mass in a triangular configuration that is fixed in the co-rotating frame. A fourth massless satellite moves in the effective potential (in this co-rotating frame) Here (x, y) is the position of the satellite in the plane of the triangle, and the fixed positions (x i , y i ) of the three bodies can be expressed in terms of their masses: [5,9].
The equilibria of the system are given by the critical points of the effective potential Ω: (8) It is known that the number of equilibrium points varies from 8 to 10 when the masses are varied (e.g. see [18,6,7] and [26]).
Using the general bounds introduced in Remark 1 for the finite dimensional case, we applied the radii polynomial approach to prove the existence of several solutions of (8), hence yielding rigorous bounds for relative equilibria of (CR4BP) Let us present a sample result. In case of equal masses m 1 = m 2 = m 3 = 1/3, the routine script equilibria.m, available at [4], computes (using Newton's method) x num = (−0.467592983336122 , 0.809894804400869), and yields the bounds Y 0 = 1.775 × 10 −15 , Z 0 = 1.23 × 10 −14 and Z 2 = 12.6987 with the choice of r * = 0.02. In this case, we obtain p(r 0 ) < 0 for any r 0 ∈ [1.78 × 10 −15 , 0.02], with p the radii polynomial defined in (7). In Figure 1 one can find several sample results, where each point has been rigorously validated using the computer program.
We now turn our attention to the elliptic PDE problems, where we apply Theorem 2.1 in an infinite dimensional setting. 3. Radially symmetric solutions of a nonlinear Laplace-Beltrami operator on the sphere. We consider the partial differential equation posed on the sphere S 2 ⊂ R 3 , where ∆ is the Laplace-Beltrami operator on the manifold (the natural geometric generalization of the Laplace operator). Here λ ≥ 0 is a parameter. The PDE (9) describes a classical nonlinear elliptic problem [19], often studied on the unit ball in arbitrary dimension and with a variety of nonlinearities.
Here we restrict attention to a quadratic nonlinearity, and we pose the problem on a sphere, cf. [11,33]. Letting x = r cos φ sin θ, y = r sin φ sin θ and z = r cos θ, letting u(x, y, z) = u(r, φ, θ), We look for solutions of (9) that are radially symmetric (with respect to rotations around the z-axis) and symmetric in about equator. This reduces the PDE to an ODE and leads to the following boundary value problem (BVP) for u = u(θ): The goal is to prove existence of solutions of (10) via the radii polynomial approach (see Theorem 2.1). The first step in doing so is to introduce a zero finding problem of the form F (x) = 0 on a Banach space.
3.1. The zero finding problem for the Laplace-Beltrami problem. Our approach is based on Taylor series, and is best posed (see Remark 2 below) on the domain [0, 1]. Hence we rescale the independent variable θ = π 2 ϑ. The algebra simplifies if we also scale the dependent variable, as well as the parameter λ: λ.

COMPUTER-ASSISTED PROOFS FOR SOLUTIONS OF PDES 67
The BVP (10) in the new variables becomes For the sake of presentation, here we have anticipated that it will be convenient to split off a factor ϑ in the second term of the differential equation, as it will allow us to deal with smooth functions that are even in ϑ only. We search for v as a power series of ϑ around zero: v(ϑ) = ∞ n=0 a n ϑ n . Let us assume for the moment that the radius of convergence of our power series is larger than 1, then the coefficients are in a Banach space of geometrically decaying coefficients. More precisely, for ν > 1, denote Given a, c ∈ 1 ν , denote by a * c the Cauchy product given component-wise by An important property of the space 1 ν is that it is a Banach algebra under the Cauchy product: The expansion for the cotangent is is the Riemann zeta function. Hence we get where the b j are the Taylor coefficients defined by The decay rate of these coefficients shows that we should restrict attention to 1 < ν < 2. After expanding all terms in (11) as Taylor series, using the Cauchy product, and equating powers, we arrive at the operator F (a) = (F n (a)) n≥0 defined as ja j for n = 1, n(n − 1)a n + (Ja * b) n +λa n−2 + (a * a) n−2 for n ≥ 2.
Here the multiplication operator J on sequence spaces is defined by Remark 2. Without the rescaling of the independent variable, the formula for F 1 corresponding to the boundary condition at θ = π 2 would have become j≥1 ja j ( π 2 ) j . Evaluating this sum (say up to some finite N ) is numerically unstable because of the high powers of π 2 . Hence the rescaling. The following result demonstrates that solving F (a) = 0 provides a solution of (11).
Proof. Let ν ∈ (1, 2) and a ∈ 1 ν . Then the power series v(ϑ) = n≥0 a n ϑ n converges uniformly for ϑ ∈ [0, 1], and similarly for the derivatives of v. The first two equations F 0 (a) = 0 and F 1 (a) = 0 imply that v satisfies the boundary conditions in (11), whereas the remaining equations F n (a) = 0, n ≥ 2 imply that v satisfies the differential equation in (11). Now that we have identified the zero finding problem F (a) = 0 to be solved in the Banach space X = 1 ν with ν ∈ (1, 2) to be chosen later, we are ready to apply the radii polynomial approach as introduced in Theorem 2.1. The first ingredient is a numerical approximationā of F (a) = 0. Given N ∈ N and a = (a n ) n≥0 ∈ X = 1 ν , denote by a (N ) = (a n ) N n=0 ∈ R N +1 the finite dimensional projection of a, and by Assume that a numerical approximationā (N ) ∈ R N +1 has been computed. We abuse slightly the notation by identifyingā (N ) ∈ R N +1 with a = (ā 0 ,ā 1 , · · · ,ā N , 0, 0, 0, · · · ) ∈ X = 1 ν . Denote by DF (N ) (ā) the Jacobian of F (N ) atā. The radii polynomial approach as introduced in Theorem 2.1 requires defining the operators A † and A. Let where the diagonal tail is chosen in view of the dominant term in (13) for large n.
Consider an (N + 1) It follows that A is injective whenever the matrix A (N ) is injective, which we verify (see Section 3.3) in order to check the injectivity assumption in Theorem 2.1. One way to visualize the operator A is as The following lemma states that condition (2) of Theorem 2.1 holds.
We leave the proof to the reader. We now introduce the formulas for the bounds Y 0 , Z 0 , Z 1 and Z 2 satisfying respectively (3), (4), (5) and (6). These bounds are derived in Sections 3.2-3.5. We first introduce an auxiliary result, used for the Z 0 and Z 2 bounds. We again leave the proof to the reader (see e.g. [15]).

3.2.
The Y 0 bound. We look for a bound Y 0 satisfying AF (ā) 1,ν ≤ Y 0 . We note that sinceā n = 0 for n > N , we have (ā * ā) n = 0 for all n > 2N . Hence, recalling the definition of A in (14), When calculating the finite sums in this expression, computing any b l involves evaluating the zeta function, which is itself an infinite series. We approximate this series by summing finitely many terms and control the remainder via a straightforward integral estimate. Concerning the final term in the expression above, since

ISTVÁN BALÁZS ET AL.
We thus set We remark that the tails of A and A † are exact inverses, hence B m,n = 0 when m > N or n > N . Letting we get from Lemma 3.3 that I − AA † B( 1 ν ) ≤ Z 0 . We note that, as in the finite dimensional case in Remark 1, it suffices to check that Z 0 < 1 to infer that the matrix A (N ) is injective, which in turn implies implies that the linear operator A is injective. In particular, if one finds an r 0 > 0 for which the radii polynomial p(r) is negative, then we see from (7) that Z 0 < 1, hence the injectivity assumption in Theorem 2.1 is then "automatically" fulfilled.

The Z 1 bound. We look for a bound A DF
For the finite part (0 ≤ n ≤ N ) we see that z 0 = 0, z 1 = j≥N +1 jh j , and z n = 0 for n = 2, . . . N.
It is not entirely straightforward to find initial guesses for solutions, i.e., starting points for applying Newton's method to the finite truncation F (N ) . To find such approximate solutions, note that the trivial solution v = 0 undergoes transcritical bifurcations atλ = π 2 n(n + 1 2 ) for n ∈ N. The solution branches that bifurcate at these parameter values can then be followed numerically (using branch following techniques) to other parameter values.
The main restriction on the Taylor series based approach presented here, is that not all solutions have an analytic extension to the complex ball of radius 1. These solutions, although real analytic on [0, 1], cannot be described by a single Taylor series around the origin. One way to overcome this is via domain decomposition (i.e. matching together several power series), but we will not pursue that here, and it is also by no means the only option.  4. Radially symmetric equilibria of the Swift-Hohenberg equation on the 3D unit ball. We consider the Swift-Hohenberg equation [27] with Dirichlet boundary conditions: Here D 1 ⊂ R 3 is the unit ball, and λ ∈ R is a parameter. The parabolic PDE (20) is a popular deterministic model for pattern formation, see e.g. [20]. It has been well studied, analytically in one spatial dimension and predominantly numerically in two spatial dimensions. Here we consider time-independent solutions in three spatial dimensions. Indeed, we will focus on radially symmetric equilibrium solutions of (20). Letting v = (∆ − 1)u, these solutions also correspond to radially symmetric equilibria of the reaction-diffusion system on the unit ball with Dirichlet boundary conditions. The method works very generally for radially symmetric equilibrium solutions in reaction-diffusion systems (cf. [25]), which are ubiquitous in models in the life sciences. This motivates us to work with the system (21) rather than with the equivalent (at the level of equilibria) scalar equation (20). As an additional benefit, the analysis below illustrates how the method of radii polynomials extends naturally to systems of equations. Looking for radially symmetric equilibria of (21), i.e., time independent solutions of the form u(x, y, z) = u(s) = u( x 2 + y 2 + z 2 ), leads to a coupled systems of ODEs: We expand the functions u and v as power series u(s) = ∞ n=0 a n s n and v(s) = ∞ n=0 b n s n . Define the coefficient sequences as a = (a n ) n≥0 and b = (b n ) n≥0 . Consider the Banach space endowed with the norm x X = max { a 1,ν , b 1,ν }. The equations for the Taylor coefficients are F 1 (x) = ((F 1 (x)) n ) n≥0 = 0 and F 2 (x) = ((F 2 (x)) n ) n≥0 = 0, given component-wise by k≥0 a k for n = 1, n(n + 1)a n − a n−2 − b n−2 for n ≥ 2, where (a 3 ) n = (a * a * a) n = n 1 +n 2 +n 3 =n n 1 ,n 2 ,n 3 ≥0 a n1 a n2 a n3 = n n1=0 a n1 n−n1 n2=0 a n2 a n−n1−n2 .
It is not difficult to derive that all odd coefficients will vanish, but we do not exploit that here. Denoting F = (F 1 , F 2 ), the problem is to find x = (a, b) ∈ X = 1 ν × 1 ν for some ν > 1 such that F (x) = 0. To achieve this, we use the radii polynomial approach as introduced in Theorem 2.1. Given N ∈ N, denote x (N ) = ((a n ) 0≤n≤N , (b n ) 0≤n≤N ) ∈ R 2N +2 and consider the finite dimensional projection Givenx ∈ R 2N +2 a numerical approximation of F (N ) (x) = 0, denote by DF (N ) (x) the Jacobian of F (N ) atx, and let us write it as The radii polynomial approach requires defining the operators A † and A. Let whose action on an element for 0 ≤ n ≤ N, δ i,1 n(n + 1)(h 1 ) n for n > N, where δ i,j is the Kronecker δ. Consider now a matrix A (N ) ∈ R (2N +2)×(2N +2) computed so that A (N ) ≈ DF (N ) (x) −1 . We decompose it into four (N +1)×(N +1) blocks: .

COMPUTER-ASSISTED PROOFS FOR SOLUTIONS OF PDES 75
Thus we define A as whose action on an element h = (h 1 , h 2 ) ∈ X is defined by (Ah) i = A i,1 h 1 + A i,2 h 2 , for i = 1, 2. The action of A i,j is defined as As in Section 3.1, to conclude that A is injective it suffices to check that A (N ) is injective. The latter is implied by Z 0 < 1 (see Section 4.2), which is automatically fulfilled when the radii polynomial is negative for some r 0 > 0. Finally, we set T (x) = x − AF (x), which indeed maps X into itself. As in Section 3, the next step is to derive explicit, computable expressions for the bounds (3), (4), (5) and (6). 4.1. The Y 0 bound. Observe first that the nonlinear term of F 2 (x) is the Cauchy product (ā * ā * ā) n−2 , which vanishes for n ≥ 3N + 3. This implies that (F 1 (x)) n = (F 2 (x)) n = 0 for all n ≥ 3N + 3. For i = 1, 2, we set which is a collection of finite sums that can be evaluated with interval arithmetic. We get and we set 4.2. The Z 0 bound. We look for a bound of the form I − AA † B(X) ≤ Z 0 . Recalling the definitions of A and A † given in (24) and (23), let B def = I − AA † the bounded linear operator represented as We remark that (B i,j ) n1,n2 = 0 for any i, j = 1, 2 whenever n 1 > N or n 2 > N . Hence we can compute the norms B i,j B( 1 ν ) using Lemma 3.3.
Then, assuming a loose a priori bound r ≤ 1 on the radius, we set where A B(X) is computed using Lemma 3.3, see also (26).

4.5.
Computer-assisted proofs. With a numerical continuation algorithm, we continued a branch of solutions of (22) that bifurcates from the zero solution at λ = (π 2 + 1) 2 . Combining the bounds Y 0 , Z 0 , Z 1 and Z 2 given in (25), (26), (27) and (28), respectively, we define the radii polynomial p(r) as in (7). We prove the existence of six solutions of (22)  We observe from these data that as we increase the parameter λ, the dimension of the projection N needs to increase while the decay parameter ν needs to decrease. This is due to the fact that the Taylor coefficients of the solutions decay slower as λ increases. Moreover, notice that the values of r min and r max are approaching each other, meaning that the proofs are getting harder and harder to obtain. This suggests that for larger parameter values a single Taylor expansion is not enough. The corresponding solutions can be seen in Figure 5. We plot in Figure 6 the solution (u, v) of (20) at λ = 500.

Conclusion.
We have seen that some nontrivial boundary value problems originating from nonlinear (geometric) PDEs can be solved in a Taylor series setting, one that requires relatively little technical machinery. Cutting off the Taylor series at some finite order and solving the associated finite dimensional algebraic system numerically, leads to an approximate solution, and we have proven that the true solution lies nearby. Indeed, based on Theorem 2.1 and computable bounds (using interval arithmetic), we can estimate the distance between approximate and true solution rigorously and explicitly. This turns the numerical computation into a mathematical statement about the PDE. There are limitations to the particular choice of Taylor series as our means of describing a solution, i.e., using monomials as our basis functions. In particular, we have seen that this puts limits on the parameter range where we can apply this approach. More generally, there is a large variety of functional analytic forumulations and numerical algorithms, adapted to the particular problem under study, which fit into the general framework of Theorem 2.1. Successful examples include domain decomposition, Fourier series, Chebyshev series, splines, finite elements, as well as combinations of these. With these tools one is able, using the paradigm illustrated by the two examples in Sections 3 and 4, to solve eigenvalue problems, find periodic and connecting orbits, continue solutions in parameters, and analyze bifurcations. Due to the nonlinear nature of these problems, you usually simply cannot get your hands on such solutions without the help of a computer. For further reading we refer to [1,2,8,14,22,30,32,35] and the references therein. Of particular note is [3], where the authors apply a related method to find asymmetric solutions of a variant of (1) with u = u(x) scalar and x in the unit disc in R 2 .