Efficient rigorous numerics for higher-dimensional PDEs via one ...

Efficient rigorous numerics for higher-dimensional PDEs via one-dimensional estimates Marcio Gameiro

∗

Jean-Philippe Lessard

†

Abstract We present an efficient rigorous computational method which is an extension of the work Analytic estimates and rigorous continuation for equilibria of higher-dimensional PDEs (M. Gameiro and J.-P. Lessard, J. Differential Equations, 249(9):2237–2268, 2010). The idea is to generate sharp one-dimensional estimates using interval arithmetic which are then used to produce high-dimensional estimates. These estimates are used to construct the radii polynomials which provide an efficient way of determining a domain on which the contraction mapping theorem is applicable. Computing the equilibria using a finite dimensional projection, the method verifies that the numerically produced equilibrium for the projection can be used to explicitly define a set which contains a unique equilibrium for the PDE. A new construction of the polynomials is presented where the nonlinearities are bounded by products of one-dimensional estimates as opposed to using FFT with large inputs. It is demonstrated that with this approach it is much cheaper to prove that the numerical output is correct than to recompute at a finer resolution. We apply this method to PDEs defined on 3D and 4D spatial domains.

1

Introduction

Computing numerical approximations of solutions of partial differential equations (PDEs) defined on spatial domains of dimension greater than two is inevitably affected by the socalled curse of dimensionality. In other words, given a fixed grid size, the number of spatial discretization points required to realize that size increases exponentially as the dimension of the domain grows. As a consequence, verifying the correctness of the numerical outputs can be a real challenge. For instance, in the context of a three-dimensional PDE, the standard approach of assessing the correctness of the numerical result simply based upon its reproducibility at different levels of refinement may be impractical. More explicitly, suppose that a numerical method to solve a three-dimensional model involves an algorithm of computational complexity n3 (e.g., computing the LU decomposition of an n × n matrix). Since the input of the algorithm is a discretization of a function defined on a three-dimensional domain, its size is n = m3 , where m represents the number of points in each dimension. A naive attempt to reproduce the result by doubling the number of mesh points in each dimension would increase the computational cost by a factor of 512, since then n3 = 512m9 . ∗ Instituto

de Ciˆ encias Matem´ aticas e de Computa¸c˜ ao, Universidade de S˜ ao Paulo, Caixa Postal 668, 13560-970, S˜ ao Carlos, SP, Brazil ([email protected]). † Universit´ e Laval, D´ epartement de Math´ ematiques et de Statistique, 1045 avenue de la M´ edecine, Qu´ ebec, QC, G1V0A6, Canada and BCAM - Basque Center for Applied Mathematics, Bizkaia Technology Park, 48160, Derio, Bizkaia, Spain ([email protected]).

1

Using this approach a one week numerical simulation could take about ten years to be verified. In an attempt to address this fundamental issue in the context of PDEs, we propose here a rigorous computational method which provides approximate solutions together with precise bounds within which exact solutions are guaranteed to exist in the mathematically rigorous sense. While the complexity of our algorithm grows exponentially as the dimension of the domain grows, the important feature of our method is that most of the extra steps of the verification does not require computing on a finer resolution than the resolution used to compute the numerical approximate solution. It is worth mentioning that several methods have been developed to address the problem of the curse of dimensionality arising from computing numerical approximations of highdimensional PDEs. For example, sparse grid approximations [1, 2, 3], sparse tensor product approximations [4, 5], tensor-product wavelet bases [6] and the Monte Carlo method [7, 8] are standard tools used in computing on high-dimensional problems. We refer to [9] and the references therein for a more extensive review. Our goal here is not to compare our approach to these more standard methods, but instead to minimize the extra computational cost needed to obtain computer-assisted proofs of existence of solutions. One of the key ingredients of the proposed approach is to use interval arithmetic [10] to generate sharp one-dimensional analytic estimates (see Appendix A) which are then used to produce high-dimensional estimates. The purpose of the high-dimensional estimates is twofold. First, they are used to construct rigorous upper bounds for the convolution terms involved in the nonlinearities of the PDE. Second, they are used to estimate the truncation error terms introduced by computing on a finite dimensional approximation of the original problem. Furthermore, since these high-dimensional estimates are solely based on products of one-dimensional estimates, their rigorous computer-assisted construction does not increase in complexity as the dimension of the domain increases. It is well known that evaluating nonlinearities of high-dimensional PDEs can be very expensive. Efficient algorithms like the fast Fourier transform (FFT) (see e.g. [11]) have been successfully applied to speed up such computations, but as mentioned earlier, when done at the frontier of computability, such algorithms can be hard to apply (see Example 1.2). In order to address this difficulty, the present method proposes to use sharp analytic estimates to bound the nonlinearities as opposed to compute them using the FFT algorithm. As it is discussed in details in Section 3.3, with this new approach, it is dramatically cheaper to prove that the numerical output is correct than to recompute at a finer resolution. Note that this result is similar in spirit to the result of [12], where it is demonstrated that the computational cost of the so-called validated continuation method is less than twice the cost of the standard continuation method alone. However only one-dimensional examples were presented in [12] and the method there did not include control of the floating point errors involved in the computations by means of interval arithmetic, meaning that the cost comparison there is a result regarding non-rigorous computations. Before proceeding further, let us discuss some of the limitations and the potential of the proposed approach. First of all, it is important to note that for the moment, our method can be applied to the class of autonomous PDEs with polynomial nonlinearities posed on rectangular domains. The reason for this restriction is because we need to know explicitly the Laplacian eigenvalues and eigenfunctions. Another restriction of our approach is that the boundary conditions of the PDEs have to be either Neumann or periodic since we are using spectral Fourier methods. However, since some 3D fluid models fall in this class of problems, we believe that the present work can be used to answer important questions in fluid dynamics. Another potential of the present method is to develop rigorous computational tools for time-periodic solutions of PDEs, where time is considered as an extra dimension.

2

For instance, the computational complexity involved in a proof of existence of a time-periodic solution of a 3D PDE would be the same than the computational complexity involved in a proof of existence of a steady state of a 4D PDE. Since in the present work, we rigorously compute steady states of a 4D PDE, one could potentially apply the proposed approach to prove existence of a time-periodic solution of a 3D fluid model. To the best of our knowledge, this is the first time that a rigorous computation of a steady state of a PDE defined on a spatial domain of dimension larger than three is obtained. The method introduced in the present work is based on the rigorous continuation method of [13], developed to prove existence of equilibria of parameter dependent PDEs of the form ut = L(u, λ) +

p X

qn un ,

(1)

n=2

in a rectangular domain Ω ⊂ Rd , where λ ∈ R is a parameter, L(·, λ) is a linear operator, and qn = qn (λ) ∈ R are the coefficients of the polynomial nonlinearities. It is demonstrated in [13] that under certain regularity conditions on the solutions, finding equilibria of (1) is equivalent to finding solutions of f (x, λ) = 0, (2) for x = {xk }k∈Zd in a Banach space X s of coefficients decaying algebraically at least as fast as { ω1s }k∈Zd , where {ωks }k∈Zd are weights with decay rate s. The map f := {fk }k∈Zd is k given component-wise by fk (x, λ) := µk (λ)xk +

p X

n=2

qn

X

xk1 · · · xkn , k ∈ Zd ,

(3)

k1 +···+kn =k kj ∈Zd

where µk (λ) are the eigenvalues of the linear operator L(·, λ). At λ = λ0 , the rigorous continuation is based on a contraction mapping argument applied to a Newton-like operator T , which depends on an approximate inverse of the Frechet derivative Df (¯ x, λ0 ) of f with respect to x, where x ¯ is a numerical approximation of a finite dimensional Galerkin approximation of (2). The method focuses on efficiently determining balls B(¯ x, r) := x ¯ +B(r) centered at x ¯ and of positive radius r in X s on which T is a contraction mapping. In the end, the proof of existence of solutions of (2) is obtained by verifying a finite number of polynomial inequalities, the so-called radii polynomials {pk (r)}k∈Zd , which provide sufficient conditions to have that T : B(¯ x, r) → B(¯ x, r) is a contraction. Since the extra cost of the proof of the rigorous continuation method of [13] is the construction of the coefficients of these polynomials, the strategy here is to reduce the cost involved in their computation. In order to present this strategy, we review the idea and the ingredients involved in the construction of the radii polynomials. For more details, we refer to [13]. The radii polynomials {pk (r)}k∈Zd are upper bounds satisfying     r x) − x ¯ k + sup DT (¯ x + b)c k − s ≤ pk (r), for every k ∈ Zd . (4) T (¯ ω b,c∈B(r) k

Assuming that one can find a radius r > 0 such that pk (r) < 0 for all k ∈ Zd , then by the contraction mapping theorem, the operator T has a unique fixed point within the ball B(¯ x, r) ⊂ X s . In order to avoid having to verify an infinite number of polynomial inequalities, one can construct a polynomial p˜M (r), independent of k, such that pk (r) = p˜M (r)

r , for k 6∈ FM , ωks 3

(5)

where M = (M1 , . . . , Md ) ∈ Nd is a computational parameter refered to as the verification dimension and FM := {k ∈ Zd | |k| < M }, where k < M and |k| denote component-wise inequalities and absolute values, respectively. The verification dimension provides the size of the finite dimensional system on which the hypotheses of the contraction mapping theorem will be verified. More explicitly, if N is the cardinality of FM , one computes N finite radii polynomials {pk (r)}k∈FM satisfying (4) and a single tail radii polynomial satisfying (5). To prove existence of steady states of (1), the following result from [13] is helpful. Lemma 1.1. Fix λ = λ0 . Consider the finite radii polynomials {pk }k∈FM satisfying (4) and the tail radii polynomial p˜M satisfying (5). If there exists r > 0 such that pk (r) < 0 for all k ∈ FM and p˜M (r) < 0, then there exists a unique x ˜ ∈ B(¯ x, r) such that f (˜ x, λ0 ) = 0. Since the finite radii polynomials {pk }k∈FM and the tail radii polynomial p˜M encode the upper bounds (4), one has a certain freedom in how to construct them. There are however two fundamental constraints. First, the verification dimension M has to be chosen large enough so that the tail radii polynomial satisfies p˜M (r) < 0 for some r > 0. Second, the finite radii polynomials have to be constructed so that their evaluation is not too expensive. In [13], the radii polynomials are constructed in a way that these two constraints compete against each other. Indeed, the nonlinearities of fk given by (3) are split as X X X xk 1 · · · xk n = xk 1 · · · xk n + xk 1 · · · xk n , (6) k1 +···+kn =k kj ∈Zd

k1 +···+kn =k kj ∈FM

k1 +···+kn =k

{k1 ,...,kn }6⊂FM

and the first finite sum is evaluated using the FFT algorithm while the second infinite sum is bounded using analytic estimates. In order to eliminate the aliasing error involved in the FFT computation of the nonlinearity, the size of the FFT inputs is larger than the verification dimension M . Hence, this approach is very expensive as the dimension of the domain grows. In order to make that important point clear, we give an explicit example. Example 1.2. Suppose that the dimension of the domain is d = 3, let M = (M, M, M ) 3 and a vector x = {xk }k∈FM ∈ R(2M −1) . To compute a rigorous bound for the cubic sum X xk 1 xk 2 xk 3 (7) k1 +k2 +k3 =k kj ∈FM

using FFT, one needs to enlarge the vector x to control the aliasing error (see e.g. [14]). Let M ∗ the smallest power of 2 such that M ∗ ≥ 4M − 1. Define M ∗ = (M ∗ , M ∗ , M ∗ ) and ∗ 3 define x ˜ = {˜ xk }k∈FM ∗ ∈ R(M ) component-wise by ( xk , if k ∈ FM x ˜k := 0, if k ∈ FM ∗ \ FM . An rigorous enclosure of (7) can be obtained by computing the Fourier transform of x ˜∈ ∗ 3 R(M ) . This computational task can be done using the FFT algorithm, which is a powerful way to compute high-dimensional convolutions (see [11]). However, in case of a 3D model, this algorithm may be hard to apply, since the value of (M ∗ )3 may be large. For example, in [13], in order to prove existence of non trivial equilibria for the 3D Cahn-Hilliard PDE, cubic convolutions of the form (7) are computed with M = 218, where M is the smallest integer so that the tail radii polynomial (5) could be successfully solved (see Figure 4 in [13]). In this case M ∗ = 210 = 1024 ≥ 4M − 1 = 871. Hence, a rigorous enclosure of (7) involves computing a FFT with inputs of size 109 . This is a serious computational task. 4

Based on the previous example, we now underline three important differences between [13] and the present work. The first and most significant improvement is a strategy to essentially avoid computing using a finer resolution during the verification. In this regard, a new construction of the polynomials is presented where the nonlinearities (7) are bounded by products of one-dimensional estimates as opposed to performing higher resolution FFT computations. Second, an explicit construction of the radii polynomials for general polynomial nonlinearities is introduced. In [13], the presentation of the polynomials is introduced only for cubic nonlinearities. Finally, since this new approach depends heavily on sharp estimates, we improve the one-dimensional estimates of [13]. The estimates are presented in Appendix A. In order to demonstrate the sharpness of the new estimates, we present in Section 3.4 a result about the existence of several equilibria of the one-dimensional Allen-Cahn equation which improves dramatically a result of [14]. The paper is organized as follows. In Section 2, the rigorous numerical method for PDEs via one-dimensional estimates is introduced. In Section 2.1, the general formulae of the radii polynomials are introduced. In Section 2.2, an analysis of the computations required to evaluate the radii polynomials is presented. In Section 3, we present some applications and the computational aspects of the method. The method is applied in Section 3.1 to prove existence (and local uniqueness) of several equilibrium solutions of a pattern formation model defined on a 4D domain. In Section 3.2, we present the cost involved in computing the numerical approximation (Section 3.2.1) and the cost involved in the rigorous verification (Section 3.2.2). In a context of a 3D PDE, it is demonstrated in Section 3.3 that it is significantly cheaper to prove that the numerical output is correct than to recompute at a finer resolution. It is shown in Section 3.4 that the new one-dimensional estimates provide a significant improvement over the method introduced in [14]. The mathematical justification of the radii polynomials is presented in Section 4 and we conclude the paper in Section 5. The Appendix contains the derivation of the one-dimensional estimates.

2

Rigorous numerics via one-dimensional estimates

In recent years, efficient algorithms from numerical analysis and scientific computing were incorporated within rigorous numerics. For example, path-following algorithms were used to handle parameter dependent equations [12, 15, 16, 17, 18] and the FFT algorithm was combined with interval arithmetic to bound nonlinearities of PDEs [14, 19]. However, like most numerical methods, the cost of the rigorous computational methods increases dramatically as the dimension of the domain grows. The goal here is to reduce the size of the inputs required for the rigorous verification method. More precisely, we focus our efforts on minimizing the computational cost required to build the radii polynomials. In this regard, we present an efficient construction of the radii polynomials whose computation essentially does not require computing on a finer resolution. The general formulae of the radii polynomials are introduced in Section 2.1 and an analysis of their computation is presented in Section 2.2. Their explicit derivation is postponed to Section 4.

2.1

General formulae for the radii polynomials

In order to perform the necessary computations one needs a finite dimensional approximation of (2). For this purpose, fix a computational parameter m = (m1 , . . . , md ) ∈ Nd , which is referred to as the numerics dimension and which represents the Galerkin projection dimension. Let Fm = {k ∈ Zd | |k| < m}. The m-dimensional Galerkin projection

5

(m)

}k∈Fm , with fk

(m)

(xFm , λ) := fk ((xFm , 0Im ), λ) = µk (λ)ck +

f (m) := {fk

fk

(m)

: Rm1 ···md × R → Rm1 ···md is given by p X

n=2

qn

X

ck 1 · · · ck n ,

(8)

k1 +···+kn =k kj ∈Fm

where xFm := {xk }k∈Fm and xIm := {xk }k6∈Fm are the finite part of size m and the corresponding infinite part of x = {xk }k∈Zd , respectively. Assume that a numerical zero x ¯Fm of f (m) at a given parameter value λ0 is computed, (m) −1 that is, f (¯ xFm , λ0 ) ≈ 0. Let Jm be a numerical approximation for the inverse of the (m) Jacobian matrix Df (¯ xFm , λ0 ), which is assumed to be invertible. Fix a computational parameter M = (M1 , . . . , Md ), which is referred to as the verification dimension. More explicitly, the parameter M determines the size of the finite radii polynomials, that is there are N = M1 M2 · · · Md finite radii polynomials. Recalling the value µk (λ) in (3), let  −1 (m)  (¯ xFm , λ0 )| k , if k ∈ Fm |Jm f    Yk := |µk (λ)−1 fk (¯ x, λ0 )|, if k ∈ FM \ Fm    0, if k ∈ 6 FM .

The one-dimensional weights

ωks

:=

(

1, if k = 0 s |k| , if k 6= 0,

(9)

(10)

are used to define the d−dimensional weights ωks :=

d Y

s

ωkjj > 0,

(11)

j=1

where s = (s1 , . . . , sd ) is the decay rate. These weights are used to define the norm kxks = sup ωks |xk |.

(12)

k∈Zd

:= {1/ωks }k∈Fm and Using (11), define ω −s := {1/ωks }k∈Zd , ωF−s m h i −s e(1) := I − J −1 Df (m) (¯ Z x , λ ) ω , Fm 0 m Fm k k

(13)

where | · | denotes component-wise absolute value and I denotes the identity matrix. Using (n) the definition of the one-dimensional estimates αk given by (34) in Section A.1 and (35) in Section A.2 we define the high-dimensional estimates by (n)

(n)

αk = αk (s, M ) :=

d Y

(n)

αkj (sj , Mj ) > 0.

(14)

j=1

The estimates (14) are used to bound convolution sums of the form (6) by terms of the (n) form αk /ωks , provided x = {xk }k∈Zd decays to zero as 1/ωks . We refer to Figure 1 for a geometrical interpretation of the decay of the high-dimensional estimates (14).

6

5000

(3)

(3)

α k 1 ,k 2

4000

α k 1 ,k 2

4000

3000

3000 2000 2000 1000

1000 60

0 0

k1

20

40 60

0

60

0 0

40 20

40 20

k2

k1

20

40 0

60

k2

4

15000

x 10

(3)

α k 1 ,k 2

9

(3)

α k 1 ,k 2

10000

6 5000

3 0 0 20

k1

60 60

20 0

40 20

40

40

60

0 0

k1

k2

20

40 60

0

k2

Figure 1: The case d = 2: two dimensional estimates αk(3) = αk(3) as product of one-dimensional 1 ,k2

estimates for M = (50, 50) and s = (s, s). The cases s = 2 (top–left), s = 3 (top-right), s = 4 (bottom-left) and s = 5 (bottom-right) are depicted. (n)

Using the definition of the one-dimensional estimates εk given by (36) in Section A.3 one defines another set of high-dimensional estimates by   sj (n)   ω α kj (n) (n) (n) εk = εk (s, M ) := ks max ε (s , M ) . (15) j j  ωk j=1,...,d  α(n) (sj , Mj ) kj kj Consider now a computational parameter ¯ 1, . . . , M ¯ d ), ¯ := (M M

such that

¯ ≤ M, m≤M

(16)

which is referred to as the FFT dimension. Using (14) and (15), let

(1)

Zk

 p p   X X  (n) −s −s n−1  n|q | |¯ x | ∗ (ω − ω ) + n|qn |k¯ xkn−1 εk , for k ∈ Fm  n s Im IM  ¯  k  n=2 n=2     p p   X X (n) n|qn | |¯ x|n−1 ∗ ωF−s + n|qn |k¯ xkn−1 εk , for k ∈ FM := ¯ \ Fm s ¯ M k   n=2 n=2     p (n)  X  αk   n|qn |k¯ xkn−1 , for k ∈ 6 FM ¯  s ωks n=2 (17)

7

and for ` ∈ {1, . . . , p − 1}, let  p    X  n − 1 h n−1−` `+1   ) n|qn | |¯ x| ∗ (ωF−s  ¯ M  ` k  n=`+1  i (n) (`+1) n−1−` + k¯ xks (` + 1)εk , for k ∈ FM Zk := ¯     p (n)  X  n − 1 αk   , for k 6∈ FM n|qn |k¯ xksn−1−` ¯   ωs `

(18)

k

n=`+1

where for sake of simplicity of the presentation, we identify xFM¯ with (xFM¯ , 0IM¯ ), xIM¯ with (0FM¯ , xIM¯ ) and xFm with (xFm , 0Im ), and we use the discrete convolution notation n times

z }| { (x )k = (x ∗ · · · ∗ x)k = n

X

xk 1 · · · xk n .

k1 +···+kn =k kj ∈Zd

Combining (13), (17) and (18), let p−1  X   (1)  −1 (`+1) Jm  Z r`+1 + Zek r,  F  m  k `=0 Zk (r) := p−1  X  1  (`+1) `+1   Z r ,  |µk (λ)| k `=0

for k ∈ Fm (19) for k ∈ FM \ Fm .

Assume that we can find µ ˜M (λ), independent of k, such that µ ˜M (λ) ≤ |µk (λ)|,

for all

k 6∈ FM .

(20)

As in [13], for M ∈ N, with M ≥ 6 and s ≥ 2 we define the one-dimensional asymptotic estimate n o (n) (n) (n) α ˜M = α ˜ M (s, M ) := max αk (s, M ) | k = 0, . . . , M , which is used to define the high-dimensional asymptotic estimate     d   Y (n) (n) (n) (n) αMj (sj0 , Mj0 ) α ˜ Mj (sj , Mj ) . α ˜M = α ˜ M (s, M ) := max 0  j0 =1,...,d    j=1

(21)

j6=j0

Define

p−1

Z˜M (r) :=

X 1 µ ˜M (λ)

p X

n|qn |k¯ xkn−1−` s `=0 n=max{`+1,2}



 n − 1 (n) ` α ˜M r . `

(22)

Definition 2.1. Recall the definition of Yk and Zk (r) satisfying (9) and (19), respectively. The finite radii polynomials {pk }k∈FM are defined by r pk (r) = Yk + Zk (r) − s , for k ∈ FM . (23) ωk Recall Z˜M (r) given by (22). The tail radii polynomial is defined by p˜M (r) = Z˜M (r) − 1.

(24)

Since the computer-assisted proofs of existence of solutions follow by verifying the hypotheses of Lemma 1.1, the extra cost of the proof is the construction of the coefficients of the radii polynomials. Let us now analyze the cost involved in their computation. 8

2.2

Analysis of the computation of the radii polynomials

¯ satisfying m ≤ M ¯ ≤ M is the parameter that governs the cost to The FFT dimension M compute the finite radii polynomials defined in (23). The convolution sums in (17) and (18) ¯ . Hence, the closer M ¯ is to are evaluated using the FFT algorithm with inputs of size M the numerics dimension m, the faster the computation of the polynomials is. The smallest ¯ and m (so that the hypotheses of Lemma 1.1 are satisfied) possible distance between M (n) is determined by the sharpness of the estimates εk given by (15). The sharper the one(n) (n) dimensional estimates αk and εk of Appendix A are, the sharper the high-dimensional (n) (n) estimates εk are. The sharper the high-dimensional estimates εk are, the smaller the ¯ is and the more efficient the computation of the finite radii polynomials FFT dimension M are. This fundamental importance of the one-dimensional estimates is the reason why we present a new improved version of the estimates from [13] in Appendix A. ¯ = m, then the cost of the FFT computations involved in the construction of the If M radii polynomials is the same as the cost of the FFT computations required to compute the numerical approximation itself. In this case, most of the extra steps of the verification method does not require computing on a finer resolution. The evaluation of Y in (9) is the ¯ 6= m, then the only computation done with a finer resolution (see Section 3.2.2). If M ¯ measures the extra level of refinement required for the proof. distance between m and M (1) (d) Finally, the computation of Zk , . . . , Zk in (17) and (18) for the cases k 6∈ FM ¯ is not affected by the curse of dimensionality. Indeed, their computation only involves evaluating the (n) (n) estimates αk , which are products of the one-dimensional estimates αk from Section A.2. ¯ = m, we have an efficient construction. As a matter of fact, Therefore, in case we can set M ¯ = m, it is demonstrated in the context of a three-dimensional in Section 3.3, setting M PDE that it is about significantly cheaper to verify rigorously the numerical outputs with interval arithmetic than to recompute at a finer resolution.

3

Applications and computational aspects

In this section, we present some applications and analyze the computational cost involved in the method. In Section 3.1, we prove existence of several equilibria of the Swift-Hohenberg PDE defined on a 4D domain. In Section 3.2, we study the computational cost of the two main parts of our method, namely the numerics and the verification. The cost involved in computing the numerical approximation is studied in Section 3.2.1 while the cost involved in the verification is studied in Section 3.2.2. In the context of a 3D domain, it is demonstrated in Section 3.3 that it is much cheaper to obtain a computer-assisted proof of existence of solutions than to recompute at a finer resolution. In Section 3.4, we test the new onedimensional estimates by studying the 1D Allen-Cahn equation.

3.1

Rigorous computation of equilibria of 4D Swift-Hohenberg

We consider a pattern formation model defined on a four-dimensional spatial domain, namely the Swift-Hohenberg equation ut = λu − (1 + ∆)2 u − u3

(25)

with even periodic boundary conditions on a rectangular bounded domain Ω = [0, 2π/L1 ] × [0, 2π/L2 ] × [0, 2π/L3 ] × [0, 2π/L4 ] ⊂ R4 . More precisely, we are interested in periodic solutions that satisfy the symmetry conditions u(y, t) = u(|y|, t), where |y| := (|y1 |, |y2 |, |y3 |, |y4 |). 9

Equation (25) describes the onset of Rayleigh-B´enard convection, and is widely used as a model for pattern formation. The parameter λ > 0, the reduced Rayleigh number, is the continuation parameter. In the context of the general PDE (1), the linear part here is given by L(u, λ) = λu − (1 + ∆)2 u and q3 = −1 is the only non zero coefficient of the polynomial nonlinearity. The equilibria of (25) can be expanded using a cosine basis {ψk }k∈N4 given Q4 P by ψk (y) := j=1 cos(kj Lj yj ). Using the expansion u = k∈Z4 xk ψk with the assumption that x|k| = xk for k ∈ Z4 , then steady state solutions of (25) correspond to solutions of fk (x, λ) := µk (λ)xk −

X

xk1 xk2 xk3 = 0, k ∈ Z4 ,

(26)

k1 +k2 +k3 =k kj ∈Z4

where µk (λ) = µ(k1 ,k2 ,k3 ,k4 ) (λ) = λ − [1 − (k12 L21 + k22 L22 + k32 L23 + k42 L24 )]2 . We now compute rigorously steady states of (25) with L1 = 1, L2 = 1.01, L3 = 1.02 and L4 = 1.03. We chose this perturbation of the 4D hypercube in order not to have degenerate bifurcations. To find the first non trivial solution, we consider a bifurcation from the trivial solution corresponding to the mode k = (1, 1, 1, 1), that is, to the eigenfunction ψk (y) = cos(L1 y1 ) cos(L2 y2 ) cos(L3 y3 ) cos(L4 y4 ). This bifurcation occurs when µk (λ) = 0, that is for λ ≈ 9.743138. We computed a branch of equilibria using a continuation method and for each point on the branch, we constructed the radii polynomials from Definition 2.1 and verified rigorously the hypotheses of Lemma 1.1. We refer to Figure 2 for the parameters used in the proofs. 0.6 0.5

kuk

0.4 0.3 0.2 0.1 0

9.76

9.78

9.8

λ

9.82

9.84

9.86

Figure 2: Branch of equilibria for Swift-Hohenberg (25) in the 4D domain Ω = [0, 2π]×[0, 2π/1.01]× [0, 2π/1.02] × [0, 2π/1.03]. For all 42 points on the plot the verification method was successful in proving the existence of a unique equilibria near the numerically computed solution. We set the ¯,M ¯,M ¯,M ¯ ), the verification ¯ = (M numerics dimension m = (m, m, m, m), the FFT dimension M dimension M = (M, M, M, M ) and the decay rate s = (2, 2, 2, 2). For the first points on the ¯ = 5 and M = 15 and these parameters increased as needed along the branch to branch, m = 5, M ¯ reach m = 5, M = 11 and M = 19 for the last point. The total running time for the whole branch is 121.17 minutes meaning that the average time per proof is about 2.88 minutes.

3.2

Computational aspects of the algorithms

We now analyze the computational cost arising from the two distinct parts of the algorithm. The first part, denoted as the numerics, computes the approximate solution of the Galerkin projection (8), while the second part, denoted as the verification, computes the coefficients of the radii polynomials from Definition 2.1 and verifies the hypotheses of Lemma 1.1.

10

3.2.1

Algorithm and computational cost involved in the numerics

To compute a numerical approximate equilibrium of the PDE model (25), we consider the Galerkin projection (8) of dimension m on which a continuation algorithm is performed. This method involves a predictor and a corrector step: given, within a prescribed tolerance, a solution x0 at parameter value λ0 , the predictor step produces an approximate equilibrium x ˜1 at a nearby parameter value λ1 = λ0 + ∆λ, and the corrector step takes x ˜1 as its input and produces, once again within the prescribed tolerance, an equilibrium x1 at λ1 . More explicitly, the predictor step requires computing a tangent vector x˙ 0 satisfying Dx f (m) (x0 , λ0 )x˙ 0 = −Dλ f (m) (x0 , λ0 ) and the predictor is given by x ˜1 = x0 + (∆λ)x˙ 0 . Since for (25), Dλ f (m) (x0 , λ0 ) = x0 this computation is trivial. Hence, the cost of the predictor step is given by computing the LU decomposition of Dx f (m) (x0 , λ0 ) in order to get the action of its inverse. The corrector step is based on the Newton-like scheme (n+1)

(n)

(0)

−1

(n)

(0)

x1 = x1 − Dx f (m) (x1 , λ1 ) f (m) (x1 , λ1 ), with x1 = x ˜1 . The computational cost of the corrector is given by the evaluations of the function f (m) at every iteration (this is (0) done with the FFT algorithm) and finding the LU decomposition of Dx f (m) (x1 , λ1 ) in order to get the action of its inverse. Since the algorithm to compute the numerical approximation requires two LU decompositions and several FFT evaluations, it is clear that the cost involved in this part of our method grows exponentially as the size of the domain grows. Let us now discuss the cost of the verification. 3.2.2

Algorithm and computational cost involved in the verification

The main computational cost involved in the verification comes from the computation of the coefficients of the radii polynomials from Definition 2.1. First note that the one-dimensional estimates are defined in terms of simple sums of size M (hence are very fast to compute) and the computation of the high-dimensional estimates are just component-wise product of one-dimensional estimates. Therefore, the computation of the estimates is not influenced −1 by the curse of dimensionality. Note that a numerical approximation Jm for the inverse of (1) (m) the Jacobian matrix Df (¯ xFm , λ0 ) is needed in order to compute Yk in (9), Zek in (13) and Zk (r) in (19). Computing this inverse is a computationally expensive step, however the continuation algorithm presented in Section 3.2.1 already computes the LU decomposition −1 of Df (m) (¯ xFm , λ0 ) when computing the tangent vector. Hence, Jm does not need to be computed during the verification part of the algorithm. The computation of the Zk (r) in ¯ . The computation of the Yk in (9) requires (19) requires computing several FFT of size M computing fk (¯ x, λ0 ) for all k ∈ FM . This is done with one FFT computation of size M . ¯ = m, the computation of the Y bound is the only extra step involved in the In case M verification that requires computing on a finer resolution. This justifies the claim in Section 1 that the important feature of our method is that most of the extra steps of the verification does not require computing on a finer resolution than the resolution used to compute the approximation.

3.3

Cost comparison for Swift-Hohenberg in 3D

We now consider the Swift-Hohenberg equation (25) defined on a three-dimensional spatial domain and present a cost comparison between the numerical computation of the equilibria and the verification method to prove that the numerical solution represents a true solution of the PDE. This equation was considered in [13], where computations in 3D could be performed only for relatively small ranges of parameters due to high computational costs of the proofs for large parameter values. Figure 3 shows the diagram of the computed solutions 11

of (25), where the range of parameters is much larger than the one considered in [13]. A plot of the level surfaces of the last point on the branch in Figure 3 is presented in Figure 3.

10

!u!

8 6 4 2

10

20

30

ν

40

50

Figure 3: Left: Branch of equilibria for the Swift-Hohenberg PDE in the 3D domain Ω = [0, 2π] × [0, 2π/1.1] × [0, 2π/1.2]. For all the points on the plot the verification method was successful in proving the existence of a unique equilibria of (25) near the numerically computed solution. We used m = (m, m, m) and s = (2, 2, 2). The continuation started by considering a nontrivial equilibrium coming out of a bifurcation. For the first point on the branch, we choose m = 8 arbitrarily. Then, we keep the same m as long as the proof is successful. If the proof fails, we increase m by one and attempt again. Using this strategy, one can prove all points on the branch and m increases as needed along the branch to reach m = 12 for the last point. Right: Level surfaces for the solution corresponding to the last point on the branch on the left at λ = 53.824306721432478. The proof ¯ = (12, 12, 12) and M = (58, 58, 58). was performed using m = M

In Figure 4, we show a cost comparison between the numerical computation (Numerics) of the equilibria and the verification method (Verification). The verification time measures the cost involved in computing the radii polynomials and verifying the hypotheses of Lemma 1.1, that is, it represents the extra computational cost involved in performing a computer-assisted proof of existence and local uniqueness of equilibria of (25). The running times correspond to a single solution on the branch of Figure 3 at λ = 7.4093499072734623, which corresponds to the fourth non trivial point rigorously computed on the branch. We recomputed the same solution for several values of m and ran the verification algorithm on that solution for each value of m. As we can see from Figure 4 the cost of proving the existence of a solution near the numerical result is comparable to the cost of computing the numerical solution itself. Hence it is much cheaper to prove the correctness of the output than to recompute at a finer resolution, since the computational cost increases dramatically with m. Let us give a more quantitative description of that claim. Consider the case m = 9 in Figure 4. It took about 5.84624 seconds to compute a numerical approximation of an equilibrium solution of (25). As mentioned in Section 1, a standard approach of assessing the correctness of a numerical result is based upon its reproducibility at a finer level of refinement. In that case, such a standard approach would require taking m = 18 to assess the correctness of the output, meaning that in the context of computing equilibria of (25), the extra computational cost to verify the numerical output would be 3084.07 seconds (see Figure 4). On the other hand, the extra computational cost of our new proposed rigorous verification method is 15.73299 seconds. This implies that the ratio between the cost of the standard validating method versus our rigorous method is 3084.07/15.73299 = 196.025676. Therefore, in this context, it is about 200 times cheaper to prove that the numerical output is correct than to recompute at a twice as fine resolution. While the above ratio (≈ 196) was computed with a Galerkin projection dimension m = (m, m, m), where m = 9, we believe that this ratio should increase as one increases m. 12

While the above ratio (⇡ 196) was computed with a Galerkin projection dime m = (m, m, m), where m = 9, we believe that this ratio should increase as one increas

3500 3000

Numerics Verification

m 8 9 11 14 18

2500

t

2000 1500 1000 500 0 8

10

12

m

14

16

Numerics 1.78967 5.84624 37.8836 323.518 3084.07

Verification 11.94395 15.73299 35.01097 317.60105 2539.81061

18

1: Running in seconds. (b)atimes Figure 4: Left: Running times for the numerics andTable the verification for solution in Figure 3. The horizontal axis represents the projection dimension m and the vertical axis refers to the time in seconds needed for the computations. The times shown refer to the same solution that was recomputed for several values of m. This shows that it is much cheaper to verify the correctness of A itdramatic improvement for Allen-Cahn the solution than to3.3 recompute with a larger value of m. Right: Running times in seconds. 1D

As mentioned in Section 1, in order to demonstrate the efficiency of the new proposed me and the sharpness of for the new estimates of 3.4 A dramatic improvement Allen-Cahn 1DAppendix A, we now present a comparison the rigorous computational method introduced in [14]. There, the construction of the As mentioned in Section 1, in order to demonstrate the efficiency of the new proposed method polynomials is done with di↵erent estimates. In particular, with the method is applied t and the sharpness of the new estimates of Appendix A, we now present a comparison one-dimensional Allen-Cahn equation the rigorous computational method introduced in [14]. There, the construction of the radii polynomials is done with different estimates.(In particular, the method is applied to the ut = ⇡ 2 / uxx + u u3 , in ⌦ = [0, 1] one-dimensional Allen-Cahn equation (
ux 3= 0, on @⌦, 2 ut = π /λ uxx + u − u , in Ω = [0, 1] ux = 0, on ∂Ω,

(27)

where > 0 is used as the continuation parameter. For (27) the Fourier basis {cos(k k= 1, continuation 2, . . .} is used and thenFor the(27) solutions of (27) represented using the expans where λ > 0 is used as 0, the parameter. the Fourier basisare {cos(kπx) | k = 0, 1, 2, . . .} is used and then the solutions of (27) are represented using the expansion 1

X

∞ u(x, t) = u0 (t) + 2 uk (t) cos(k⇡x). X uk (t) cos(kπx). u(x, t) = u0 (t) + 2 k=1 k=1

So (27) takes the form

So (27) takes the form

duk = µk (λ)uk − dt

du Xk = uµkk1(uk)u k 2 uk3 , dt

k1 +k2 +k3 =k

X

uk1 uk2(28) u k3 ,

k1 +k2 +k3 =k

where µk (λ) = 1 − π 2 k 2 /λ thethe linear operator inof(27). The radii where µk (are) the = 1eigenvalues ⇡ 2 k 2 / ofare eigenvalues the linear operator in (27). The polynomials introduced in Definition 2.1 were constructed in the context of (28) and several polynomials introduced in Definition 2.1 were constructed in the context of (28) and se equilibria of (27) were proved to exist by verifying the hypotheses of Lemma 1.1. In Figure 5, equilibria of (27) were proved to exist by verifying the hypotheses of Lemma 1.1. In Figu we present a comparison of the branches computed with the construction of [14] and the weusing present a comparison of the branches computed withusing the the construction of [14] and branches computed the new method introduced in this paper. Note that branches using the for newthe method introduced this paper. Note that usin present method, we did not computed encounter any failure method. This is a in significant present method, we[14], did where not encounter any failure for the method. This is a signifi improvement compare to the result of all computations eventually failed at improvement compare to the result of [14], where all computations eventually faile fairly small parameter values. Another improvement of theparameter method can be seen by comparing the minimal numerics fairly small values. dimension (Galerkin projection dimension) m is required verify of Another improvement that of the methodtocan be the seenhypotheses by comparing the minimal Gal Lemma 1.1 using the construction of the radii polynomials from [14] and the construction projection dimension m that is required to verify the hypotheses of Lemma 1.1 usin of the radii polynomials done in this paper. We refer to Figure 6 for a comparison. Notice

construction of the radii polynomials from [14] and the construction of the radii polyno done in this paper. Notice that using the present method, the finite dimensional Gal projection can be taken13105 times smaller than in [14] while still producing a proof. T a dramatic improvement. 13

1.4

1.5

1.2 1

kuk

kuk

1

0.8 0.6

0.5

0.4 0.2 0 0

200

400

λ

600

800

0 0

1000

400

800

1200

λ

1600

Figure 5: Branches from [14] (left) and some computed using the proposed method (right). that using the present method, the finite dimensional Galerkin projection can be taken 105 times smaller than in [14] while still producing a proof. This is a dramatic improvement. 5

6

x 10

60 50

4

40

m

m

5

3

30

2

20

1

10

0 0

200

400

λ

600

800

0 0

1000

400

800

λ

1200

1600

Figure 6: The minimal required dimension of the Galerkin projection m as a function of the continuation parameter λ along: (left) the blue dotted branch of equilibria of Allen-Cahn from Figure 5 (left) which was proved using the method from [14]; (right) the branch (3) of equilibria of Allen-Cahn from Figure 5 (right) which was proved using the present method.

4

Justification of the formulae of the radii polynomials

In this section, we justify and describe the construction of the radii polynomials that are defined in Section 2.1. The first step in the construction of the radii polynomials is to find a numerical zero for f (m) at a given parameter value λ0 , that is x ¯Fm such that f (m) (¯ xFm , λ0 ) ≈ 0. Defining x ¯ := (¯ xFm , 0Im ) one expects that f (¯ x, λ0 ) ≈ 0 provided x ¯Fm is computed using a sufficiently large numerics dimension (Galerkin projection dimension) m. This assumption is mathematically justified by the fact that since equilibria of parabolic PDEs are solutions which exist globally in time, the linear term L in (1) regularizes the solutions and hence equilibria are very smooth (see e.g. [20]). Hence, the coefficients of the Fourier expansion should decay very fast (exponentially) and so by increasing the dimension of the Galerkin projection, one expects to approximate a true solution of the infinite dimensional PDE. Using the numerical solution x ¯, the equation f (x, λ) = 0, where f is given componentwise by (3), is now transformed into an equivalent fixed point problem for a Newton-like −1 operator about x ¯. For this purpose, let Jm be a numerical approximation for the inverse

14

of the Jacobian matrix Df (m) (¯ xFm , λ0 ), which is assumed to be invertible. Recalling the definition of the high-dimensional weights in (11), we define the Banach space X s = {x | kxks < ∞} ,

(29)

consisting of sequences with algebraically decaying tails according to the rate s. The linear operator J −1 on X s is defined by i h −1   Jm (xFm ) , if k ∈ Fm h i k J −1 (x) := 1  k  xk , if k 6∈ Fm . µk (λ)

Notice that J −1 is an approximation for the inverse of Df (¯ x, λ0 ). Defining T (x) := x − J −1 f (x, λ0 ),

one can readily see that finding zeros of f is equivalent to finding fixed points of T (see [13]). The idea is to uniquely enclose fixed points of T into closed balls B(¯ x, r) in X s centered s at x ¯. One can easily check that the closed ball of radius r in X , centered at the origin, is given by Y  r r  B(r) := B(0, r) = − s, s . ωk ωk d k∈Z

The closed ball of radius r centered at x ¯ is given by B(¯ x, r) = x ¯ + B(r). As proved in Lemma 3.3 in [13], to show that T : B(¯ x, r) → B(¯ x, r) is a contraction mapping in X s , one needs bounds Yk and Zk satisfying   x) − x ¯ k ≤ Yk , for every k ∈ Zd , (30) T (¯   x + b)c k ≤ Zk (r), for every k ∈ Zd . sup DT (¯

(31)

b,c∈B(r)

Recalling that the radii polynomials must satisfy (4), it is sufficient to compute Yk and Zk satisfying (30) and (31) to construct them.

4.1

Computation of Yk

In order to compute the upper bounds Yk , it is important to note that, since x ¯ is such that x ¯k = 0 for k 6∈ Fm , we have (¯ xn )k = 0 for every n ≥ 1 and k ∈ Zd such that kj ≥ n(mj − 1) + 1 for some 1 ≤ j ≤ d. Hence, choosing M such that Mj ≥ p(mj − 1) + 1 for all 1 ≤ j ≤ d, we have that fk (¯ x, λ0 ) = 0 for every k ≥ M . Therefore, since T (¯ x) − x ¯= −J −1 f (¯ x, λ0 ), we can define Y = {Yk }k∈Zd as  −1 (m)  |Jm f (¯ xFm , λ0 )| k , if k ∈ Fm    Yk := |µk (λ)−1 fk (¯ x, λ0 )|, if k ∈ FM \ Fm    0, if k ∈ 6 FM .

15

4.2

Computation of Zk

In order to compute Zk we denote J˜m := Df (m) (¯ xFm , λ0 ) and introduce the linear operator on X s h i  J˜m (xF ) , if k ∈ Fm h i m ˜ k J(x) := µ (λ)x , k if k 6∈ Fm , k k

:= {1/ωks }k∈Fm . which is an approximate inverse to J −1 . Define ω −s := {1/ωks }k∈Zd and ωF−s m We write     x + b, λ0 ) − J˜ c, (32) DT (¯ x + b)c = I − J −1 J˜ c − J −1 Df (¯

and notice that the first term on the right hand side of (32) is zero for k 6∈ Fm , and is very small for k ∈ Fm admitting the following upper bound h i h
i (1) −1 =: rZek , Df (m) (¯ xFm , λ0 ) ωF−s I − J −1 J˜ c ≤ r I − Jm m k

k

where | · | denotes component-wise absolute value. Notice that

p n−1 X X n − 1  
Df (¯ x + b, λ0 )c k = µk (λ)ck + nqn x ¯n−1−` ∗ b` ∗ c k ` n=2 `=0

and   ˜ Jc = k

 p X
  µk (λ)ck + nqn x ¯n−1 ∗ cFm k , for k ∈ Fm    µk (λ)ck ,

n=2

for k 6∈ Fm .

Set b = ru and c = rv, with u, v ∈ B(1). Now, for the sake of simplicity of presentation, we identify uFM¯ with (uFM¯ , 0IM¯ ) and uIM¯ with (0FM¯ , uIM¯ ), for u ∈ X s . Hence, we have p h  i X (j) Df (¯ x + b, λ0 ) − J˜ c = Ck rj , k

where (1)

Ck :=

 p 
X  nqn x ¯n−1 ∗ vIm k , for k ∈ Fm    n=2

p  X
   nqn x ¯n−1 ∗ v k ,   n=2

and for ` ∈ {1, . . . , p − 1}, (`+1)

Ck

:=

p X

n=`+1 (j)

j=1

nqn

for k 6∈ Fm ,

 
n−1 x ¯n−1−` ∗ u` ∗ v k , ` (j)

(j)

for ` = 1, . . . , p − 1.

Next upper bounds Zk so that |Ck | ≤ Zk , for j = 1, . . . , p, are constructed. For the (1) (p) cases k ∈ FM ¯kj = 0 for ¯ , we split the sums involved in the Ck , . . . , Ck . Recalling that x

16

kj 6∈ Fm we can write  p p X X 

 n−1  nqn x ¯ nqn x ¯n−1 ∗ vIM¯ k , for k ∈ Fm ∗ (vIm − vIM¯ ) k +    n=2 n=2 (1) Ck = p p  X X

   nqn x nqn x for k ∈ FM ¯n−1 ∗ vFM¯ k + ¯n−1 ∗ vIM¯ k , ¯ \ Fm .   n=2

n=2

Now, for a given ` ∈ {1, . . . , p − 1} and k ∈ FM ¯ , consider the splitting   p  X n − 1  n−1−` (`+1) Ck = nqn x ¯ ∗ u`FM¯ ∗ vFM¯ ` k n=`+1   p i X n − 1  n−1−` h ` + . nqn x ¯ ∗ u ∗ v − u`FM¯ ∗ vFM¯ ` k n=`+1

Hence, for general k ∈ Zd , one has the following component-wise upper bounds  p   X  −s   n|qn | |¯ x|n−1 ∗ (ωI−s − ω )  I ¯ m M  k n=2   p     X  −s n−1 −s n−1  , for k ∈ Fm + n|q |k¯ x k (ω ) ∗ ω  n s IM  ¯  k  n=2     p X   (1) n|qn | |¯ x|n−1 ∗ ωF−s Ck ≤ ¯ M  k   n=2 p    X   −s n−1 −s n−1 + n|q |k¯ x k (ω ) ∗ ω , for k ∈ FM ¯ \ Fm  n s IM ¯  k   n=2    X p 
   n|qn |k¯ xkn−1 (ω −s )n k , for k 6∈ FM ¯  s n=2

and for ` ∈ {1, . . . , p − 1},

(`+1)

Ck

 p    X n − 1  n−1−`  −s `+1   n|q | |¯ x | ∗ (ω ) n  FM ¯  ` k   n=`+1     p  X  n−1   + n|qn |k¯ xksn−1−` × ` ≤ n=`+1 i  h   −s `+1 −s n−1−` −s `+1  (ω ) ∗ (ω ) − (ω ) ,  FM  ¯  k    p  
n−1  X  n|qn |k¯ xkn−1−` (ω −s )n k ,  s  ` n=`+1

for k ∈ FM ¯ for k 6∈ FM ¯.

The fundamental step required to finalize the construction of the radii polynomials is to find efficient ways to bound the nonlinear convolutions. This is precisely where the sharp one-dimensional estimates presented in Appendix A become fundamental. Fix ` ∈ {0, . . . , p}, n ∈ {max{` + 1, 2}, . . . , p} and k ∈ FM ¯ . Using the definition of the one(n) dimensional estimates αk from Appendix A we define the high-dimensional estimates (n)

(n)

αk = αk (s, M ) :=

d Y

j=1

17

(n)

αkj (sj , Mj ).

(n)

Using the definition of the one-dimensional estimates εk given by (36) in Appendix A we define another set of high-dimensional estimates   s (n)   ωkjj αk (n) (n) (n) εk = εk (s, M ) := s max ε (s , M ) . j j k  ωk j=1,...,d  α(n) (sj , Mj ) j kj

From Lemma 2.2 in [13], we have that   `+1  

−s −s n−1−` −s `+1 ω ∗ ω − ωFM¯ = k

X

(n)

−s ωk−s 1 · · · ωkn ≤ (` + 1)εk .

k1 +···+kn =k

{k1 ,...,k`+1 }6⊂FM ¯

Similarly, given n ∈ {2, . . . , p} and k ∈ / FM ¯ , we can use the Lemma 2.1 from [13] to get that (ω −s )n As defined in (17) and (18), letting

(1)

Zk




(n)

k

≤

αk . ωks

 p p   X X  (n) −s −s n−1  − ω ) + n|qn |k¯ xkn−1 εk , for k ∈ Fm n|q | |¯ x | ∗ (ω  n s IM Im  ¯  k n=2 n=2    X p p    X (n) n|qn | |¯ x|n−1 ∗ ωF−s + n|qn |k¯ xkn−1 εk , for k ∈ FM := ¯ \ Fm s ¯ M k   n=2 n=2     p (n)  X  αk   n|qn |k¯ xkn−1 , for k ∈ 6 FM ¯  s ωks n=2

and, for ` ∈ {1, . . . , p − 1},  p    X  n − 1 h n−1−` −s `+1   n|q | |¯ x | ∗ (ω ) n  FM ¯  ` k   n=`+1 i (n) (`+1) n−1−` + k¯ xk s (` + 1)εk , for k ∈ FM Zk := ¯     (n) p  X  n − 1 αk   n|qn |k¯ xkn−1−` , for k 6∈ FM ¯,  s  ` ωs k

n=`+1

we recall the definition of Zk (r) in (19) to obtain  p−1   X −1 (`+1)  `+1   Jm Z  Z (r) = r  Fm  k k   `=0    (1) +Z˜k r, sup DT (¯ x + b)c k ≤  b,c∈B(r)   p−1  X  (`+1) `+1   Z (r) = |µ−1 r ,  k |Zk  k `=0

Now assume that we can find µ ˜M (λ), independent of k, such that µ ˜M (λ) ≤ |µk (λ)|, 18

for all

k 6∈ FM .

for k ∈ Fm for k ∈ FM \ Fm .

Then, for k 6∈ FM , we have   x + b)c k sup DT (¯

b,c∈B(r)

≤ =

p−1 X `=0

p−1 X `=0

≤ =

(`+1) `+1

|µk (λ)−1 |Zk |µk (λ)−1 |

p X

r

n|qn |k¯ xkn−1−` s

n=max{`+1,2}

r ˜ ZM (r) ωks  p−1 X r  1 ωks µ ˜M (λ)

p X



 (n) n − 1 αk `+1 r ` ωks

n|qn |k¯ xksn−1−`

`=0 n=max{`+1,2}



  n − 1 (n) `  α ˜M r . `

We have all the bounds required to define the radii polynomials given by (23) and (24).

5

Conclusion

In this paper, a new construction of the radii polynomials is proposed where the nonlinearities are bounded by products of one-dimensional estimates as opposed to using FFT with large inputs. An explicit construction of the radii polynomials for general polynomials nonlinearities is presented. The method is applied to prove existence of equilibria of PDEs defined on 3D and 4D spatial domains. It is demonstrated that with this new approach it is much cheaper to prove that the numerical output is correct than to recompute at a finer resolution. We believe that the present method can be used to develop rigorous computational tools for time-periodic solutions of PDEs. Since in the present work, we rigorously compute steady states of a 4D PDE, one could potentially apply the proposed approach to prove existence of time-periodic solutions of 3D fluid models. Hence, the present work has the potential to be used to answer important questions in fluid dynamics.

6

Acknowledgements

We would like to thank the anonymous referees for helpful suggestions that improved the presentation of the paper.

A

Appendix: sharper one-dimensional estimates

In this section some improvements of the one-dimensional estimates of [13] are introduced. The reason for presenting these estimates is because the one-dimensional estimates need to be as sharp as possible since they play a fundamental role in the construction of the radii polynomials. We present in Figure 7 a comparison between the estimates (for a cubic nonlinearity) presented here and the ones presented in [13]. Consider a decay rate s ≥ 2, a computational parameter M ≥ 6 and define, for k ≥ 3, 

k γk = γk (s) := 2 k−1

s



4 ln(k − 2) π 2 − 6 + + k 3

19



2 1 + k 2

s−2

.

(33)

Lemma A.1. For s ≥ 2 and k ≥ M ≥ 6 we have k−1 X

k1 =1

A.1

k−1 X ks ωks = ≤ γk ≤ γM . s ωks1 ωk−k k1s (k − k1 )s 1 k1 =1

Sharper quadratic estimates (2)

(2)

For k ∈ Z, we define αk = αk (s, M ) by

(2)

αk :=

 M X  2 1   1 + 2  2s + M 2s−1 (2s − 1) ,  ω   k1 =1 k1    M  X  2ωks 2ωks   +  s s  ωk1 ωk+k1 (k + M + 1)s M s−1 (s − 1)

for k = 0

k1 =1

k−1 X  ωks    , +2 +  s  ω ωs   k1 =1 k1 k−k1   M  X  1 2   + s−1 + γM ,   2+2 ωks1 M (s − 1)

(34) for 1 ≤ k ≤ M − 1 for k ≥ M,

k1 =1

and for k < 0,

(2)

(2)

αk := α|k| . Lemma A.2 (Quadratic estimates). Given s ≥ 2 and M ≥ 6. Then, for any k ∈ Z, (2)

X

α 1 ≤ ks . s s ωk1 ωk2 ωk =k

k1 +k2 kj ∈Z

Proof. For k = 0, X

1 s ωs ω k1 k2 =0

=

M X 1 1+2 +2 ωk2s1 k1 =1

k1 +k2 kj ∈Z

≤

∞ X

k1 =M +1

Z ∞ M X 1 1 dx ≤1+2 + 2s ωk2s1 ωk2s1 x M k1 =1

M (2) X 1 2 α + 2s−1 = 0s . 1+2 2s ωk1 M (2s − 1) ω0 k1 =1

For 1 ≤ k ≤ M − 1, and recalling the one-dimensional weights (10), " M # ∞ k−1 X X X 1 2ωks ωks 1 X 2ωks 2 = + + s+ s s s ωks1 ωks2 ωks ωks1 ωk+k ωks1 ωk+k ω0 ωks1 ωk−k 1 1 1 k +k =k 1

k1 =1

2

kj ∈Z

1 ≤ s ωk 1 ≤ s ωk

" "

M X

k1 =1 M X

k1 =M +1

2ωks 2ωks + s s ωk1 ωk+k1 (k + M + 1)s

2ωks ωs ωs k1 =1 k1 k+k1

+

(k + M +

Z

∞

M

k1 =1

k−1 X dx ωks + 2 + s s xs ωk1 ωk−k 1

2ωks 1)s M s−1 (s 20

k1 =1

− 1)

+2+

k−1 X

k1 =1

#

ωks s ωks1 ωk−k 1

#

(2)

=

αk . ωks

Finally, for k ≥ M , one gets from Lemma A.1 that " ∞ # k−1 X X X 1 1 2 ωks ωks = 2 + s+ s s ωks1 ωks2 ωks ωks1 ωk+k ω0 ωks1 ωk−k 1 1 k +k =k 1

k1 =1

2

kj ∈Z

" M # ∞ X 1 X 1 2 1 +2 + s + γk 2 ωks ωs ωs ω0 k1 =M +1 k1 k1 =1 k1 " M # Z ∞ X 1 dx 2 1 +2 + s + γM 2 s ωks ωs ω0 M x k1 =1 k1 " M # (2) αk 1 X 2 2 + .  + 2 + γ = M ωks ωks1 M s−1 (s − 1) ωks

≤ ≤ ≤

A.2

k1 =1

k1 =1

Shaper general estimates (n)

(n)

For n ≥ 3, define αk = αk (s, M ) by  M −1 (n−1) (n−1) X  αk1 2αM  (n−1)   α + 2 + , for k = 0 0   ωk2s1 (M − 1)2s−1 (2s − 1)   k =1 1    M −k (n−1) s k−1 (n−1)  X X αk(n−1) ωks αk+k1 ωk αM ωks  1   + +  s ωs s ωs s (M − k)s−1 (s − 1)  ω (M + 1) ω   k1 =1 k1 k+k1 k1 =1 k1 k−k1   (n−1) s M (n−1) s  X  α ω αM ωk  (n−1) (n−1) k k1  + αk + α0 , + +  s ωs s M s−1 (s − 1) (n) ω (M + k + 1) αk := k1 =1 k1 k+k1   for 1 ≤ k ≤ M − 1        M (n−1)  X  2αM 1 (n−1)   α + + Σ∗  M s s−1 (s − 1)  ω M   k1 =1 k1   (n−1)  M  X αk1  (n−1) (n−1)   + + αM + α0 , for k ≥ M   ωks1 k1 =1

(35)

and for k < 0,

(n)

αk

(n)

:= α|k| .

Lemma A.3 (General estimates). Given s ≥ 2 and M ≥ 6. Then, for any k ∈ Z, X

k1 +···+kn

kj ∈Z

(n)

αk 1 s · · · ωs ≤ ωs . ω k1 k kn =k

Proof. In what follows, the estimates are obtained similarly as in the proof of Lemma A.2 (n−1) (n−1) with the difference that we often use the fact αk ≤ αM , for all k ≥ M (see e.g. Remark A.1 in [13]). For k = 0, X

k1 +···+kn kj ∈Z

1 s · · · ωs ω k1 kn =0

≤

(n−1)

α0

+2

M −1 X k1 =1

21

(n−1)

αk1 ωk2s1

(n−1)

+

(n)

2αM α = 0s . (M − 1)2s−1 (2s − 1) ω0

For k > 0, X

k1 +···+kn

kj ∈Z

1 s ωk1 · · · ωksn =k

∞ X

≤

k1 =1

"

(n−1)

1 αk+k1 s ωks1 ωk+k 1 (n−1)

+

1 αk ω0s ωks

k1 =1

(n−1)

≤ αM

"M −k (n−1) X αk+k ωks

(n−1)

αk+k1 1 ≤ s s ωks1 ωk+k ω k 1

1

k1 =1

"

(n−1)

1 αk−k1 s ωks1 ωk−k 1

#

+

∞ X

k1 =1

"

1 s ωk+k 1

(n−1)

αk1 ωks1

#

(n−1)

(n−1)

k1 =1

+

k−1 X

1 α0 ωks ω0s

+

Consider k ∈ {1, . . . , M − 1}. Since αk ∞ X

#

s ωks1 ωk+k 1

. , for all k ≥ M , we have

# (n−1) αM ωks + . (M + 1)s (M − k)s−1 (s − 1)

Similarly, ∞ X

k1 =1

(n−1)

αk1 1 ≤ s s ωks1 ωk+k ω k 1 (n)

From the definition of αk

"

# (n−1) M (n−1) X αk1 ωks αM ωks + . s ωks1 ωk+k (M + k + 1)s M s−1 (s − 1) 1

k1 =1

for k ∈ {1, . . . , M − 1}, X

k1 +···+kn

kj ∈Z

For k ≥ M ,

∞ X

k1 =1

(n)

α 1 ≤ ks . s s ωk1 · · · ωkn ωk =k

" # (n−1) M (n−1) X αk+k1 αM 1 1 (n−1) ≤ s αM + s−1 . s ωks1 ωk+k ωk ωks1 M (s − 1) 1 k1 =1

Using Lemma A.1, k−1 X

k1 =1

(n−1)

αk1 s ωks1 ωk−k 1

=

M −1 X k1 =1

≤ ≤ ≤

(n−1) (n−1) k−1 αk1 1 X ωks αk1 + s s ωks1 ωk−k ωks ωks1 ωk−k 1 1 k1 =M

(n−1) αk1 1
s ωks ω s 1 − kk1 k1 =1 k1 M −1 X

1 ωks 1 ωks

"M −1 X k1 =1

"M −1 X k1 =1

(n−1)

As defined in [13], let α ˜M

(n−1)

+

αM ωks

ωks s ω ωs k1 =M k1 k−k1

(n−1)

αk1

ωks1 1 − (n−1)

αk1


k1 s M

Ms

ωks1 M − k1

+

k−1 X

(n−1) αM

(n−1)


s + αM

M −1

k−1 X

X ωks ωks − s s s ω ω ω ωs k1 =1 k1 k−k1 k1 =1 k1 k−k1 !# M −1 X 1 1 γM − =: s Σ∗a . ωks1 ωk k1 =1

n o (n−1) = max αk (s, M ) | k = 0, . . . , M . Hence,

k−1 X

k1 =1

(n−1)

(n−1) αk1 α ˜M 1 ≤ γM =: s Σ∗b . s s ωks1 ωk−k ω ω k k 1

22

!#

Letting Σ∗ := min {Σ∗a , Σ∗b }, one gets that ∞ X

(n−1)

k1 =1

αk1 1 ≤ s s s ωk1 ωk+k1 ωk

"

k−1 X

k1 =1

(n−1)

αk1 1 ≤ s Σ∗ . Also, s ωks1 ωk−k ω k 1

# (n−1) M (n−1) X αk1 αM + s−1 . ωks1 M (s − 1)

k1 =1

Combining the above inequalities, we get, for the case k ≥ M , " M (n−1) X X 2α 1 1 1 (n−1) ≤ s αM + s−1M + Σ∗ s s s ω · · · ω ω ω M (s − 1) k1 kn k k1 +···+kn =k k1 =1 k1 kj ∈Z # (n−1) M (n) X αk1 αk (n−1) (n−1) + + α + α .  = 0 M ωks1 ωks 100

80

80

60

60

αk

(3)

100

αk

(3)

k1 =1

40

20

40

20

0 0

50

100

150

k

0 0

200

50

100

k

150

200

Figure 7: Comparison of the αk(3) (in red) defined in the present paper and the αk(3) (in blue) presented in [13] with M = 200. The cases s = 2 (left) and s = 3 (right) are depicted.

A.3

(n)

New formulas for εk

¯ ≥ 6 we define, for 0 ≤ k ≤ M ¯ − 1, Given s ≥ 2 and M ≥ M (n)

εk

(n)

¯ , M ) := = εk (s, M

M −k X

¯ k1 =M

+

M +k X

¯ k1 =M

and for k < 0

(n−1)

αk+k1 s ωks1 ωk+k 1

(36)

  (n−1) (n−1) αk1 −k αM 1 1 + + ωks1 ωks1 −k (M + 1)s (s − 1) (M − k)s−1 (M + k)s−1 (n) ¯ , M ) := ε(n) (s, M ¯ , M ). εk (s, M |k|

¯ ≥ 6, for n ≥ 3 and 0 ≤ |k| ≤ M ¯ − 1 we have that Lemma A.4. Given s ≥ 2 and M ≥ M ! n X Y (1) (n) (n) ck1 · · · ckn ≤ ` Ai ε k . k1 +···+kn =k i=1 max{|k ¯ |,...,|k |}≥M 1

`

23

Proof. We have that X

ωks1

k1 +···+kn =k

¯ max{|k1 |,...,|k` |}≥M

1 ≤` · · · ωksn k

X

1 +···+kn

1 s · · · ωs , ω k1 kn =k

¯ |k1 |≥M

and X

k1 +···+kn

¯ − M X 1 1 s · · · ωs = s ω ω k1 kn k1 =k k1 =−∞

¯ |k1 |≥M

X

k2 +···+kn =k−k1 ∞ X

+

¯ k1 =M

≤ ≤

¯ k1 =M

"

M −k X

∞ X αk+k1 (n−1) + α M s ωks1 ωk+k 1

∞ X

¯ k1 =M

+

≤

(n−1)

k1 =M −k+1

(n−1) M +k X αk1 −k ωks1 ωks1 −k ¯ k1 =M

+

+

X

ωs k2 +···+kn =k−k1 k2

#

(n−1)

(n−1) αk+k1 s ωks1 ωk+k 1 ¯ k1 =M M −k X

(n−1)

αk+k1 α + s k1 −k s s ωk1 ωk+k1 ωk1 ωks1 −k

1 ωks1

1 ωks2 · · · ωksn

(n−1)

+ αM

1 · · · ωksn

1 s ωks1 ωk+k 1

∞ X

k1 =M +k+1

1 ωks1 ωks1 −k

(n−1) αk1 −k ωks1 ωks1 −k ¯ k1 =M M +k X

  (n−1) αM 1 1 + . (M + 1)s (s − 1) (M − k)s−1 (M + k)s−1 (n)

The result then follows from the definition of εk .



References [1] Volker Barthelmann, Erich Novak, and Klaus Ritter. High dimensional polynomial interpolation on sparse grids. Adv. Comput. Math., 12(4):273–288, 2000. Multivariate polynomial interpolation. [2] F. Nobile, R. Tempone, and C. G. Webster. An anisotropic sparse grid stochastic collocation method for partial differential equations with random input data. SIAM J. Numer. Anal., 46(5):2411–2442, 2008. [3] F. Nobile, R. Tempone, and C. G. Webster. A sparse grid stochastic collocation method for partial differential equations with random input data. SIAM J. Numer. Anal., 46(5):2309–2345, 2008. [4] Christoph Schwab and Rob Stevenson. Adaptive wavelet algorithms for elliptic PDE’s on product domains. Math. Comp., 77(261):71–92 (electronic), 2008. [5] Tobias von Petersdorff and Christoph Schwab. Numerical solution of parabolic equations in high dimensions. M2AN Math. Model. Numer. Anal., 38(1):93–127, 2004. 24

[6] Michael Griebel and Stephan Knapek. Optimized tensor-product approximation spaces. Constr. Approx., 16(4):525–540, 2000. [7] Paul Glasserman. Monte Carlo methods in financial engineering, volume 53 of Applications of Mathematics (New York). Springer-Verlag, New York, 2004. Stochastic Modelling and Applied Probability. [8] George S. Fishman. Monte Carlo. Springer Series in Operations Research. SpringerVerlag, New York, 1996. Concepts, algorithms, and applications. [9] Hans-Joachim Bungartz and Michael Griebel. Sparse grids. Acta Numer., 13:147–269, 2004. [10] Ramon E. Moore. Interval analysis. Prentice-Hall Inc., Englewood Cliffs, N.J., 1966. [11] Charles Van Loan. Computational frameworks for the fast Fourier transform, volume 10 of Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992. [12] Sarah Day, Jean-Philippe Lessard, and Konstantin Mischaikow. Validated continuation for equilibria of PDEs. SIAM J. Numer. Anal., 45(4):1398–1424 (electronic), 2007. [13] Marcio Gameiro and Jean-Philippe Lessard. Analytic estimates and rigorous continuation for equilibria of higher-dimensional PDEs. J. Differential Equations, 249(9):2237– 2268, 2010. [14] Marcio Gameiro, Jean-Philippe Lessard, and Konstantin Mischaikow. Validated continuation over large parameter ranges for equilibria of PDEs. Math. Comput. Simulation, 79(4):1368–1382, 2008. [15] Jan Bouwe van den Berg, Jean-Philippe Lessard, and Konstantin Mischaikow. Global smooth solution curves using rigorous branch following. Math. Comp., 79(271):1565– 1584, 2010. [16] Sarah Day, Yasuaki Hiraoka, Konstantin Mischaikow, and Toshi Ogawa. Rigorous numerics for global dynamics: a study of the Swift-Hohenberg equation. SIAM J. Appl. Dyn. Syst., 4(1):1–31 (electronic), 2005. [17] Stanislaus Maier-Paape, Ulrich Miller, Konstantin Mischaikow, and Thomas Wanner. Rigorous numerics for the Cahn-Hilliard equation on the unit square. Rev. Mat. Complut., 21(2):351–426, 2008. [18] Marcio Gameiro and Jean-Philippe Lessard. Rigorous computation of smooth branches of equilibria for the three dimensional Cahn-Hilliard equation. Numer. Math., 117(4):753–778, 2011. [19] Yasuaki Hiraoka and Toshiyuki Ogawa. An efficient estimate based on FFT in topological verification method. J. Comput. Appl. Math., 199(2):238–244, 2007. [20] Roger Temam. Infinite-dimensional dynamical systems in mechanics and physics, volume 68 of Applied Mathematical Sciences. Springer-Verlag, New York, second edition, 1997.

25