A Rigorous Integrator and Global Existence for Higher-Dimensional Semilinear Parabolic PDEs via Semigroup Theory

Studying global dynamics and global existence (forward in time) of solutions of higher-dimensional semilinear parabolic partial differential equations (PDEs) is a central problem in mathematics. Perhaps the most striking example is the famous millennium prize problem which raises the question of global existence of solutions of initial value problems (IVPs) in the three-dimensional Navier–Stokes equations [1]. From the perspective of dynamical systems and as made clear from the pioneering work of Poincaré on the three-body problem [2], to understand the global dynamics of a given nonlinear differential equation, it is crucial to grasp the existence of asymptotic objects such as equilibria, periodic orbits and connections between them. This is especially true for parabolic PDEs modeling phenomena in material science, spatial ecology and fluids mechanics, where patterns arise as asymptotic limits, often within global attractors, of solutions of initial value problems. For all its importance and popularity, the problem consisting of understanding the global dynamics of a given PDE is a notoriously difficult task since the model is nonlinear and naturally lead to the notion of an infinite dimensional dynamical system.

In this paper, we introduce a constructive method to solve rigorously initial value problems of scalar semilinear parabolic PDEs of the form

where \(p \in \{0,1\}\) and N is a polynomial satisfying both \(N(0)=0\) and its Fréchet derivative \(DN(0)=0\), and we provide a computational framework to show that solutions exist globally forward in time. While in this work we focus on the scalar case, that is \(u=u(t,x) \in \mathbb {R}\), extending our method to systems (i.e. when \(u=u(t,x) \in \mathbb {R}^n\)) should be rather straightforward, the main limitation essentially being the computational cost. Given a dimension \(d \in \{1,2,3\}\) and domain sizes \(L_1,\dots ,L_d\), we assume that the domain of the Eq. (1.1) is given by the hyper-rectangular domain

which we supplement with periodic boundary conditions. In other words, the geometry of the domain on which we solve the PDE (1.1) is the d-torus \(\mathbb {T}^d\). It is worth mentioning that one case of the Navier–Stokes millennium prize problem considers the model defined on the domain \(\mathbb {T}^3\). The parameters \(\lambda _0,\lambda _1,\lambda _2 \in \mathbb {R}\) in (1.1) are chosen so that the equation is parabolic. The assumption of the PDE being semilinear implies that the degree p of the Laplacian in front of the nonlinear term N is less than the one of the differential linear operator \(\lambda _0 + \lambda _1 \Delta + \lambda _2 \Delta ^2\). While slightly restrictive, the class of equations (1.1) includes many known models including the Swift–Hohenberg, Cahn-Hilliard, Ohta–Kawasaki and phase-field-crystal (PFC) equations. Moreover, to simplify the presentation and the computations, we impose the even symmetry \(u(t, -x) = u(t, x)\) on the solutions. Note that our approach can readily be adapted to the case of odd boundary conditions \(u(t, -x) = -u(t, x)\) or to the general case of periodic boundary conditions without symmetries.

The general method that we introduce in this paper falls in the category of computer-assisted proofs (CAPs) in nonlinear analysis and dynamical systems.

As exemplified by the early pioneering works on the Feigenbaum conjectures [3] and on the existence of chaos and global attractor in the Lorenz equations [4,5,6], and by the more recent works on Jones’ and Wright’s conjectures in delay equations [7, 8], chaos in the 1D Kuramoto-Sivashinsky PDE [9], instability proof in Poiseuille flow [10], bifurcating solutions for 3D Rayleigh–Bénard problems [11], equilibria in the 3D Navier–Stokes (NS) equations [12], solutions of NS on unbounded strips with obstacle [13], 3D gyroids patterns in materials [14], periodic orbits in NS [15], blowup in Euler equations on the cylinder [16] and imploding solutions for 3D compressible fluids [17], it is undeniable that the role of CAPs is starting to play a major role in the analysis of differential equations, dynamical systems and PDEs. We refer the interested reader to the book [18] and the survey papers [19,20,21] for more details.

Due to the interest in understanding global dynamics in semilinear parabolic PDEs and perhaps motivated by the Navier–Stokes millenium problem, it is not surprising that developing CAPs techniques for initial value problems has been a rather active field in the last twenty years. Examples consist of the topological method based on covering relations and self-consistent bounds [9, 22,23,24,25,26], the \(C^1\) rigorous integrator of [27], the evolution operator and semigroup approach of [28,29,30,31], Nakao’s projection method [32], the finite element and self-consistent bounds approach [33], the fully spectral Fourier–Chebyshev approach [34], the Chebyshev interpolation in time method [35] and the finite element discretization based approach of [18, 36].

To the best of our knowledge, all the above methods (expect [28] which considers 2D examples) have focused on PDEs defined on one-dimensional domains. Our contribution here is the presentation of a general constructive method to compute solutions of IVPs for higher-dimensional semilinear parabolic equations via semigroup theory, a method which can then be used to prove rigorously global existence (in time) of solutions. Our approach begins by using the periodic boundary conditions on \(\Omega \) and expand solutions of (1.1) with Fourier expansions in space leading to an equivalent infinite-dimensional system of ordinary differential equations of the form

whose phase space is given by a Banach algebra of Fourier coefficients endowed with analytic norms. For a given time step \(\tau >0\), solving an IVP on a time interval \(J\mathop {=}\limits ^\mathrm{{def}}(0,\tau ]\) with initial condition \(a(0) = \varphi \) reduces to find a zero of the map

Using a finite dimensional projection of the map F, we compute a numerical approximation \(\bar{a}(t)\) of the IVP using Chebyshev series expansions in time, that is \(F(\bar{a}(t)) \approx 0\). We then develop estimates to obtain CAPs for the existence of the solutions of the linearized equations about the approximation \(\bar{a}(t)\) from which a solution map operator is obtained. While semigroup theory (cf. [37], e.g.) guarantees the existence of the solution operator, our computational and rigorous construction provides an explicit control over this operator, a feature which is new in the context of higher-dimensional PDEs. Using the solution operator, the linear operator \(DF(\bar{a}(t))\) is inverted “by hand" and then used to define the Newton-like operator \(T(a(t)) = DF(\bar{a}(t))^{-1}(DF(\bar{a}(t))a(t)-F(a(t)))\). Constructing explicit and computable estimates, and using the time step \(\tau \) to control the contraction rate, we derive a sufficient condition in the form of a polynomial inequality to demonstrate rigorously that T is a self-map and a contraction on a ball centered at \(\bar{a}(t)\) in an appropriate Banach space of time-dependent Fourier coefficients. The unique fixed point together with its rigorous bounds provide the local inclusion of the solution. Applying this technique for multiple time steps leads to constructive proofs of existence of solutions over long time intervals and, when combined with explicit constructions of trapping regions, can be used to prove global existence. While our computer-assisted approach is an extension of the works [28, 30, 31, 38] to the case of higher-dimensional semilinear evolution equations, we believe it provides nontrivial contributions.

First, we believe this is the first time that CAPs techniques are used to prove global existence of solutions of three-dimensional PDEs with initial conditions far from equilibrium. Second, the evolution operator approach presented here is generalized to PDEs with derivatives in the nonlinear terms, which is made possible thanks to the introduction of geometric weights in the Banach algebra. Third, the smoothing property of the evolution operator controls wrapping effects, allowing the method to be applied for multiple time steps to get existence of solutions over long time intervals. Last by not least, the semigroup approach provides a cost-effective way to generalize our approach to multi-steps (or domain decompositions). More precisely, solving the linearized problems at each step can be made a-priori and independently (i.e. the process is naturally parallelizable), which implies that the computational cost is additive in time rather than multiplicative.

The paper is organized as follows. In Sect. 2, we introduce the fixed point operator T whose fixed point provides a solution to a given initial value problem. The operator T requires inverting the linearization about the numerical solution, which is done in Sect. 3, where we construct the evolution operator via a solution map of the linearized problem. In Sect. 4, we generalize our method to a multi-step approach by considering a coupled system of the zero-finding problem over multiple time steps. In Sect. 5, we propose a strategy to demonstrate global existence of solutions to IVPs via the mechanism of convergence towards an asymptotically stable equilibrium solution. This approach is exemplified by the application of the Swift–Hohenberg equation, illustrating the efficacy of computer-assisted proofs of global existence for higher dimensional PDEs. Section 6 is devoted to presenting the result using our rigorous integrator for IVPs of the 2D Ohta–Kawasaki equation, in particular addressing the handling of derivatives in the nonlinear term.

In this section, we consider an initial value problem (IVP) associated to (1.1), that is

where \(u_0(x)\) is a given initial data. Supplementing the PDE with periodic boundary conditions with the even symmetry \(u(t, -x) = u(t, x)\) implies that the unknown solution of the IVP is expanded by the Fourier series in space variables

where \({\varvec{k}}=(k_1,\dots ,k_d)\) and \(x=(x_1,\dots ,x_d)\). The symmetry assumption on u implies that the function u is represented by the cosine series

with \(\alpha _{{\varvec{k}}}=\alpha _{k_1,\dots ,k_d} \mathop {=}\limits ^\mathrm{{def}}2^{\delta _{k_1,0}}\cdots 2^{\delta _{k_d,0}}\),

where \(\delta _{i,j}\) is the Kronecker delta function. From now on we always assume that \(a_{{\varvec{k}}} = a_{|{\varvec{k}}|}\) holds for all \({\varvec{k}}\in \mathbb {Z}^d\) from the cosine symmetry and use the notation of the component-wise absolute value, that is \(|{\varvec{k}}|\mathop {=}\limits ^\mathrm{{def}}(|k_1|,\dots ,|k_d|)\).

Remark 2.1

Boldface characters will be used throughout this paper to denote multi-indices, with \({\varvec{k}}=(k_1,\dots ,k_d)\) representing (non-negative) integers \(k_j\) for \(j=1,\dots ,d\). Furthermore, component-wise inequalities are used such that \({\varvec{k}}< {\varvec{n}}\) denotes \(k_j<n_j\) for all \(j=1,\dots ,d\). Similarly, the notations \({\varvec{k}}\le {\varvec{n}}\), \({\varvec{k}}>{\varvec{n}}\), and \({\varvec{k}}\ge {\varvec{n}}\) mean component-wise inequalities.

Plugging the Fourier expansion (2.2) in the nonlinearity \(\Delta ^p N(u)\) of the general model (1.1) results in

where \(q=2p \in \{0,2\}\) and

\(a = (a_{{\varvec{k}}})_{{\varvec{k}}\ge 0}\), and \(\mathcal {N}_{{\varvec{k}}}(a)\) is a nonlinear term involving discrete convolutions of a such that \(\mathcal {N}_{{\varvec{k}}}(a) = \mathcal {N}_{|{\varvec{k}}|}(a)\), \(\mathcal {N}_{{\varvec{k}}}(0)=0\) and \(D \mathcal {N}_{{\varvec{k}}}(0)=0\).

Plugging the Fourier series (2.2) into (1.1) and assuming (2.3) holds, we obtain the following general infinite-dimensional system of ordinary differential equations (ODEs)

where

is determined by the linear term \(\lambda _0 + \lambda _1 \Delta + \lambda _2 \Delta ^2\) of (1.1).

Our goal is to determine the time-dependent Fourier coefficients \(a(t)\mathop {=}\limits ^\mathrm{{def}}(a_{{\varvec{k}}}(t))_{{\varvec{k}}\ge 0}\) that solve (2.4) with an initial condition \(a(0)=(\varphi _{\varvec{k}})_{{\varvec{k}}\ge 0}\). The function a(t) will then provide, via the Fourier expansion (2.2), the solution of the IVP (2.1) with \(u_0(x) = \sum _{{\varvec{k}}\ge 0} \alpha _{{\varvec{k}}} \varphi _{{\varvec{k}}} \cos (k_1 L_1 x_1)\cdots \cos (k_d L_d x_d)\).

Before considering the infinite-dimensional problem, we explain how to obtain a numerical approximation of a(t) via a finite-dimensional truncation of the ODEs (2.4). To this end, we define a finite set of indices of “size \({\varvec{N}}=(N_1,\dots ,N_d)\)” as \({\varvec{F_N}}\mathop {=}\limits ^\mathrm{{def}}\{{\varvec{k}}\ge 0:{\varvec{k}}<{\varvec{N}}\}\). We note that \({\varvec{F_N}}=F_{N_1}\times \dots \times F_{N_d}\), where \(F_{N_j}\mathop {=}\limits ^\mathrm{{def}}\{k_j\ge 0:k_j<N_j\} = \{0,\dots ,N_j - 1\}\).

Truncating (2.4) to the size \({\varvec{N}}\), we fix a step size \(\tau >0\) and numerically solve the following truncated \(N_1N_2 \cdots N_d\)-dimensional IVP

The numerical solution of (2.5) is denoted by \((\bar{a}^{({\varvec{N}})}_{{\varvec{k}}}(t))_{{\varvec{k}}\in {\varvec{F_N}}}\) and has Fourier expansion

Moreover, we assume that \(\bar{a}^{(\varvec{N})}_{{\varvec{k}}}(t)\) is a polynomial given by

where \(T_\ell \) is the \(\ell \)-th order Chebyshev polynomial of the first kind defined on \([0,\tau ]\) and n denotes the size of Chebyshev coefficients. To get the coefficients \((\bar{a}_{\ell ,{\varvec{k}}})_{\begin{array}{c} 0\le \ell <n,\\ {\varvec{k}}\in {\varvec{F_N}} \end{array}}\) we use a standard Chebyshev interpolation technique (see Appendix A for more details).

Using the numerical solution (2.6), we define an approximate solution \(\bar{a}_{{\varvec{k}}}(t)\) of (2.5) as a natural extension of \((\bar{a}_{{\varvec{k}}}^{(\varvec{N})}(t))_{{\varvec{k}}\in {\varvec{F_N}}}\), that is

Remark 2.2

(Choice of basis functions) It is worth noting that by focusing our analysis on PDEs defined over hyper-rectangular domains with periodic boundary conditions, the use of a Fourier expansion in space is a natural choice. Extending this to alternative spatial expansions would require a non-trivial generalization, which lies beyond the scope of this work. On the other hand, we chose to solve the finite-dimensional IVP (2.5) using a Chebyshev expansion in time due to the favorable approximation properties of Chebyshev polynomials and our positive experience with them. However, other expansions, such as Taylor series, could also be employed to solve (2.5).

2.1 Function Space and Banach Algebra

Having introduced the finite dimensional IVP (2.5) together with the representation of its numerical approximation (2.6) expressed as a finite linear combination of Chebyshev polynomials, we are now ready to introduce several function spaces to validate the numerical approximation (2.7).

Let us define a Banach space of multi-index sequences as

where the weight \(\omega _{{\varvec{k}}}\) for \({\varvec{k}}\ge 0\) is defined by

The choice of the weight in (2.9) is to ensure that \(\ell _{\omega }^1\) is Banach algebra under the discrete convolution, that is

holds for all \(a, b\in \ell _{\omega }^1\), where the discrete convolution of a and b is defined by

Denote the time step \(J=(0,\tau ]\) for the fixed step size \(\tau >0\) and define a function space of the time-dependent sequences as \(C(J;\ell _{\omega }^1)\), which is the Banach space with the norm

The operator norm for a bounded linear operator \(\mathcal {B}\) from \(\ell _{\omega }^1\) to \(\ell _{\omega }^1\) is defined by

We also introduce subspaces of all multi-index sequences denoted by

Note that these subspaces are also Banach spaces, but are not Banach algebras under the discrete convolution. In addition, \(\ell _{\omega }^1=\ell _{\omega ,0}^1\) (\(q=0\)).

2.2 The Zero-Finding Problem for the Initial Value Problem

To solve rigorously the infinite-dimensional set of ODEs (2.4) on the sequence space (2.8), let us define three operators as follows. First, define \(\mathcal {L}:D(\mathcal {L})\subset \ell _{\omega }^1\rightarrow \ell _{\omega }^1\) to be a densely defined closed operator acting on a sequence \(\phi =(\phi _{{\varvec{k}}})_{{\varvec{k}}\ge 0}\) as

where the domain of the operator \(\mathcal {L}\) is defined by \(D(\mathcal {L})\mathop {=}\limits ^\mathrm{{def}}\left\{ \phi =(\phi _{{\varvec{k}}})_{{\varvec{k}}\ge 0}:\mathcal {L}\phi \in \ell _{\omega }^1\right\} \). Second, define \(\mathcal {Q}:\ell _{\omega }^1\rightarrow \ell _{\omega ,-q}^1\) to be a multiplication operator acting on a sequence \(\phi =(\phi _{{\varvec{k}}})_{{\varvec{k}}\ge 0}\) as

Third, define \(\mathcal {N}:\ell _{\omega }^1\rightarrow \ell _{\omega }^1\) to be a Fréchet differentiable nonlinear operator acting on a sequence \(\phi =(\phi _{{\varvec{k}}})_{{\varvec{k}}\ge 0}\) as

which depends on the nonlinear term of the target PDE, and satisfies \(\mathcal {N}(0)\) and \(D\mathcal {N}(0)=0\). Hence, for any fixed \(\psi \in \ell _{\omega }^1\), there exists a non-decreasing function \(g:(0,\infty )\rightarrow (0,\infty )\) such that

Using the above operators, we define a zero-finding problem whose root provide a solution of the ODEs (2.4) with initial condition \(\varphi \in \ell _{\omega }^1\). Let \(X\mathop {=}\limits ^\mathrm{{def}}C(J;{\ell _{\omega }^1})\), \(Y\mathop {=}\limits ^\mathrm{{def}}C(J;{\ell _{\omega ,-q}^1})\) and \(\mathcal {D}\mathop {=}\limits ^\mathrm{{def}}C\left( J; D(\mathcal {L})\right) \cap C^1\left( {(0,\tau )};\ell _{\omega }^1\right) {\cap \,C(J;\ell _{\omega }^1)}\). Note that X is the Banach space with the norm defined in (2.11). The map F acting on \(a\in {\mathcal {D}}\) is defined by

We find it convenient to denote the first and second components of (2.13) respectively by \((F(a))_1={\dot{a}}-\mathcal {L}a-\mathcal {Q}\,\mathcal {N}(a)\) and \((F(a))_{2}=a(0)-\varphi \). From the above construction, solving the initial value problem (2.1) reduces looking for the solution of the zero-finding problem for the map F given in (2.13).

As it is often the case in nonlinear analysis, the strategy to prove existence (constructively) of a solution of \(F=0\) is to turn the problem into an equivalent fixed-point problem, which we solve using a Newton–Kantorovich type theorem. The idea is to look for the fixed point of a Newton-like operator of the form \(a \mapsto a - DF(\bar{a})^{-1} F(a)\), where \(DF(\bar{a})\) is the Fréchet derivative of the map F at the numerical approximation \(\bar{a}\) given component-wise by (2.7). This construction is done next.

2.3 Exact Inverse of the Linearization and Definition of the Newton-Like Operator

To use a Newton-like operator to validate the root of F, we consider first the Fréchet derivative of the map F at \(\bar{a}\) defined in (2.7). For any \(h \in {\mathcal {D}}\), the Fréchet derivative is given by

where \(D\mathcal {N}(\bar{a})\) is the Fréchet derivative of \(\mathcal {N}\) at \(\bar{a}\).

To construct the inverse of \(DF(\bar{a})\), consider \((p,\phi )\in Y\times \ell _{\omega }^1\) and note that the unique solution of \(DF(\bar{a}) h = (p,\phi )\), that is of the linearized problem

is given by the variation of constants formula

In the previous formula, \(\{U(t, s)\}_{0\le s\le t\le \tau }\) denotes the evolution operator (cf. [37]) defined by a solution map of the first equation with \(p=0\) in (2.15). More precisely, we consider an initial value problem which is represented by

where \(\phi \in \ell _{\omega }^1\) is any initial data. Here s is also a parameter since the evolution operator U(t, s) is a two parameter family of bounded operators. We define the evolution operator by \(U(t,s)\phi \mathop {=}\limits ^\mathrm{{def}}b(t)\) for \(0\le s\le t\), where b is the solution of (2.17). Furthermore, it is known that the evolution operator satisfies the following properties:

1. (i)

   \(U(s,s)=\textrm{Id}\) (identity operator) and \(U(t,r)U(r,s)=U(t,s)\) for \(0\le s\le r\le t\le \tau \)

2. (ii)

   \((t,s)\mapsto U(t,s)\) is strongly continuous for \(0\le s\le t\le \tau \).

We show how to construct rigorously the evolution operator in Sect. 3.1. Assuming that this is done, the formula (2.16) gives us the inverse of \(DF(\bar{a})\) via the formula

Using the explicit expression for \(DF(\bar{a})^{-1}\), we define the Newton-like operator by

Remark 2.3

The Newton-like operator defined in (2.19) seems to be defined only on \({\mathcal {D}\subset X}\) by the formula \(a \mapsto a-DF(\bar{a})^{-1}F(a)\), but it can be defined on X using the following fact:

where the linear part including derivatives is canceled out. From the property of \(DF(\bar{a})^{-1}:Y\times \ell _{\omega }^1\rightarrow {\mathcal {D}}\), a fixed-point, say \(\tilde{a}\), of T satisfies \(\tilde{a}\in {\mathcal {D}}\). Combining this fact with the continuous property of the evolution operator, the resulting zero of the map defined in (2.13) also satisfies \(\tilde{a}\in C([0,\tau ];\ell _{\omega }^1)\).

Remark 2.4

The inverse map \(DF(\bar{a})^{-1}\) restores the smoothness of the solution. To demonstrate the recovery of smoothness from Y to X, letting \(\phi =0\) and \(p=\mathcal {Q}\psi \in Y\) (\(\psi \in X\)) in (2.18), one can show the inverse map denoted by

is bounded in X. By using the uniform control of the evolution operator given in (3.2), which is constructed in Sect. 3.1, it follows that

Moreover, it is obvious that the function b satisfying (2.15) belongs to both \(C\left( J; D(\mathcal {L})\right) \) and \(C^1\left( (0,\tau );\ell _{\omega }^1\right) \). This fact yields the smoothing recovery of the inverse map \(DF(\bar{a})^{-1}:Y\times \ell _{\omega }^1\rightarrow {\mathcal {D}}\).

Remark 2.5

It is also worth mentioning that there is another natural way of obtaining the operator T directly considering the IVP written by

Using the evolution operator, it follows that

and this formulation directly makes sense on X.

Now that the Newton-like operator T is defined, the remaining task in validating the local existence of a solution of \(F=0\) given in (2.13) is to show that \(T:B_{J}(\bar{a}, \varrho ) \rightarrow B_{J}(\bar{a}, \varrho )\) is a contraction mapping, for some \(\varrho >0\), where

Note that \(B_{J}(\bar{a}, \varrho ) \subset X\) is the ball of radius \(\varrho \) in X centered at the numerical approximation \(\bar{a}(t) = (\bar{a}_{{\varvec{k}}}(t))_{{\varvec{k}}\ge 0}\in X\) defined in (2.7).

In this section, we provide a numerical method to validate the solution of the zero-finding problem defined in (2.13). The main theorems for this method are based on Theorem 3.2 in Sect. 3.1 and Theorem 3.15 in Sect. 3.2. The techniques presented here are closely related to the work in [30].

First, in Sect. 3.1, we obtain a uniform bound for the evolution operator over the simplex \(\mathcal {S}_{J} \mathop {=}\limits ^\mathrm{{def}}\{(t, s): 0< s\,{<} \,t\le \tau \}\), which is presented in Theorem 3.2 (and Corollary 3.4). More precisely, by setting two initial data for (2.17) as \(b(s) = \phi \) and \({\mathcal {Q}} \psi \in \ell _{\omega ,{-q}}^1\), we obtain two computable constants \(\varvec{W^{\mathcal {S}_{J}}},~\varvec{W_q^{\mathcal {S}_{J}}}>0\) such that

where \(\gamma \in (0,1)\) is a parameter to be determined in Sect. 3.1.2. If \(q=0\) in the target ODEs (2.4), one can take \(\gamma =0\), resulting in \(\varvec{W^{\mathcal {S}_{J}}}=\varvec{W_{0}^{\mathcal {S}_{J}}}\). Note that the bound (3.2) simply shows that by introducing the weight \((t-s)^\gamma \) in time, the loss of regularity in space is compensated by the evolution operator. See also Remark 2.4.

Second, in Sect. 3.2, the above computable constants \(\varvec{W^{\mathcal {S}_{J}}}\) and \(\varvec{{W_q^{\mathcal {S}_{J}}}}\) are used to bound the action of the linear operator \(DF(\bar{a})^{-1}\) (see Lemma 3.13). This control is then used to verify that the Newton-like operator (2.19) has a unique fixed point in \(B_{J}(\bar{a}, \varrho )\), as defined in (2.20). This yields a rigorous validation of the local existence of the solution to the IVP (2.1). This result is formalized in our second main theorem, Theorem 3.15, which provides a hypothesis guaranteeing the local existence of the solution to the zero-finding problem in (2.13). This hypothesis can be numerically verified using interval arithmetic. Therefore, in the present method, we numerically confirm that the hypothesis holds, thereby validating the solution of the zero-finding problem.

3.1 Explicit Uniform Bounds of the Evolution Operator

In this section, we introduce explicit uniform bounds for the evolution operator. We present an explicit computable formula (see Equation (3.18)) for the bound \(\varvec{{W_q^{\mathcal {S}_{J}}}}\) satisfying (3.2) (see Theorem 3.2). Then, in Corollary 3.4, we also present an explicit formula (see Equation (3.28)) the bound \(\varvec{W^{\mathcal {S}_{J}}}\) satisfying (3.1).

Remark 3.1

While the theoretical existence of the evolution operator U(t, s) is well-established (cf. [37, Section 5.6, Theorem 6.1]), obtaining a rigorous enclosure of its action for a given \(\tau \) is difficult in practice. In practical computations, for a specified \(\bar{a}\) and \(\tau \), an explicit and rigorous bound on U(t, s) is achieved first by employing computer-assisted proofs to solve a finite dimensional projection of the linearized problem (2.17) and then by using theoretical estimates on the tail. This methodology yields the uniform bounds (3.1) and (3.2), used to solve the Cauchy problem.

Let \({\varvec{m}}= (m_1,\dots ,m_d)\ge 0\) be a size of the finite set \({\varvec{F_m}}\), where we suppose that each component \(m_i\) is small such that \({\varvec{m}}<{\varvec{N}}\). We define a finite-dimensional (Fourier) projection \(\varPi ^{(\varvec{m})}: \ell _{\omega }^1 \rightarrow \ell _{\omega }^1\) acting on \(\phi \in \ell _{\omega }^1\) as

We denote \(\phi ^{(\varvec{m})}\mathop {=}\limits ^\mathrm{{def}}\varPi ^{(\varvec{m})} \phi \in \ell _{\omega }^1\) and \(\phi ^{(\infty )}\mathop {=}\limits ^\mathrm{{def}}(\textrm{Id}-\varPi ^{(\varvec{m})}) \phi \in \ell _{\omega }^1\). Thus, any \(\phi \in \ell _{\omega }^1\) is uniquely decomposed by \(\phi =\phi ^{(\varvec{m})}+\phi ^{(\infty )}\).

We decompose \(b\in X\) in (2.17) into a finite part \(b^{({\varvec{m}})}\) and an infinite part \(b^{(\infty )}\) using the Fourier projection defined in (3.3), we construct the solution b. In the followings, two initial conditions are considered, that is \(b(s)=\phi \) and \(b(s)=\mathcal {Q}\psi \), but only the case of \(b(s)=\mathcal {Q}\psi \) is presented, because the other case is derived as a corollary to the presented case. For this purpose we consider yet another Cauchy problem with respect to the sequence \(c(t)=\left( c_{{\varvec{k}}}(t)\right) _{{\varvec{k}}\ge 0}\).

with the initial data \(c^{({\varvec{m}})}(s)={{\mathcal {Q}}}\psi ^{({\varvec{m}})}\) and \(c^{(\infty )}(s)={{\mathcal {Q}}}\psi ^{(\infty )}\). We remark that one can independently separate these system into the finite part and the infinite part. So we define the solution operator of (3.4) by extending the finite dimensional fundamental solution matrix \(C^{({\varvec{m}})}(t,s)\mathop {=}\limits ^\mathrm{{def}}\left( C^{({\varvec{m}})}_{{\varvec{k}},{\varvec{j}}}(t,s)\right) _{{\varvec{k}},{\varvec{j}}\in {\varvec{F_m}}}\) satisfying \(\sum _{{\varvec{j}}\in {\varvec{F_m}}} C^{({\varvec{m}})}_{{\varvec{k}},{\varvec{j}}}(t,s){\left( {\mathcal {Q}}\psi \right) }_{{\varvec{j}}}= c_{{\varvec{k}}}(t)\). Similarly, the solution operator of (3.5) is defined by using the solution map \(C^{(\infty )}(t,s){{\mathcal {Q}}\psi ^{(\infty )}}\mathop {=}\limits ^\mathrm{{def}}(c_{{\varvec{k}}}(t))_{{\varvec{k}}\notin {\varvec{F_m}}}\) for \((t,s)\in {\mathcal {S}_{J}}\). We extend the action of the operator \(C^{({\varvec{m}})}(t,s)\) (resp. \(C^{(\infty )}(t,s)\)) on \(\ell _{\omega }^1\) by introducing the operator \(\bar{U}^{({\varvec{m}})}(t,s)\) (resp. \(\bar{U}^{(\infty )}(t,s)\)) given by

At this point, the existence of the operator \(C^{(\infty )}(t,s)\) is not guaranteed, but will be verified via a method of rigorous numerics for the fundamental solution matrix introduced in Sec 3.1.1 and via an analytic approach for the evolution operator in Sec 3.1.2, respectively. Using the solution operators (3.6) and (3.7), we present the solutions of (3.4) and (3.5) as

Returning to consider the solution b(t), the solutions of (2.17) is represented by using the Fourier projection

Here, from the variation of constants, we have

By using this formula (3.8), the following theorem provides a sufficient condition (see Eq. 3.17) under which we have an explicit formula for the uniform bound of \(\varvec{{W_q^{\mathcal {S}_{J}}}}\) satisfying (3.2).

Theorem 3.2

Let \((t, s) \in {\mathcal {S}_{J}}\) and \(\bar{a}\) be fixed. Let \({\varvec{m}}\ge 0\) be the size of the Fourier projection such that \(\mu _{{\varvec{k}}}<0\) holds for \({\varvec{k}}\not \in {\varvec{F_m}}\). For such \({\varvec{m}}\ge 0\), define \(\mu _*\) as

and assume that there exists a constant \({W_{m,q}^{\mathcal {S}_{J}}}>0\) such that

Assume also that \(\bar{U}^{(\infty )}(t,s)\) satisfies

Moreover, let us assume the existence of constants \({W_{\infty ,q}^{\mathcal {S}_{J}}}>0\), \({\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}\ge 0\), \(\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}\ge 0\) satisfying

respectively. Assume also the existence of two constants \(\mathcal {E}_{m,\infty }^{J}\) and \({\mathcal {E}_{\infty ,m}^{J}}\) such that

hold for each \(t\in J\). Then, if

a uniform bound of the evolution operator U(t, s), that is \(\varvec{{W_q^{\mathcal {S}_{J}}}}>0\) in (3.2), is obtained by

where \(\Vert \cdot \Vert _1\) denotes the matrix 1-norm.

Note that the bounds \({{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\) and \({\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\) become smaller if an appropriate size of the Fourier projection \({\varvec{m}}\ge 0\) is chosen. In particular, the fact that \({\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\) can be made small by increasing \({\varvec{m}}\) helps insuring that condition (3.17) of Theorem 3.2 is satisfied. This can be intuitively understood from the following simplified case: if we assume \(W^{(\infty )}_{{q}}(t,s)\sim e^{\mu _*(t-s)}\) (where \(\mu _*\) is defined in (3.9)), then \({{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\sim 1/|\mu _*|\) and \({\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\sim 1/\mu _*^2\) hold. Therefore, as the absolute value of \(\mu _*\) become large, it is more likely that condition (3.17) will be satisfied.

Remark 3.3

Theorem 3.2 asserts that there exists the evolution operator \(\{U(t, s)\}_{0\le s\le t\le \tau }\) on \(\ell _{\omega }^1\) and it is uniformly bounded in the sense of (3.2) over the simplex \({\mathcal {S}_{J}}\). The specific bounds, including the \({W_{m,q}^{\mathcal {S}_{J}}}\), \(W^{(\infty )}_{{q}}(t,s)\), \({W_{\infty ,q}^{\mathcal {S}_{J}}}\), \({{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\), and \({\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\) bounds, are presented in the following sections of the paper. The \({W_{m,q}^{\mathcal {S}_{J}}}\) bound is discussed in Sect. 3.1.1 via rigorous numerics for the fundamental solution matrix of (3.4). The \(W^{(\infty )}_{{q}}(t,s)\) bound is presented in Sect. 3.1.2 considering the solution map of (3.5). The remaining bounds are presented in Sect. 3.1.3. Additionally, the \({\mathcal {E}_{m,\infty }^{J}} \) and \({\mathcal {E}_{\infty ,m}^{J}} \) bounds are determined for each case. These bounds are further introduced in Sect. 5.3.

Proof

Taking the \(\ell _{\omega }^1\) norm of (3.8), it follows from (3.10), (3.11), (3.15), and (3.16) that

Plugging each estimate in the other one, we obtain

and

Since \(\kappa =1-{W_{m,q}^{\mathcal {S}_{J}}}{\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}} {\mathcal {E}_{m,\infty }^{J}} {\mathcal {E}_{\infty ,m}^{J}} >0\) holds from the sufficient condition of the theorem, \(\mathcal {X}\) norm bounds of \(b^{(\varvec{m})}\) and \(b^{(\infty )}\) are given by

Finally, we have

\(\square \)

Corollary 3.4

Under the same assumptions in Theorem 3.2, assume that \(\bar{U}^{(\infty )}(t,s)\) satisfies

Assume also that there exists \({W_{\infty }^{\mathcal {S}_{J}}}>0\), \({{\overline{W}}_{\infty }^{\mathcal {S}_{J}}}\ge 0\), \({{{\overline{W}}}_{\infty ,q}^{\prime \mathcal {S}_{J}}}\ge 0\), \({{\overline{\overline{W}}}_{\infty ,q}^{\prime \mathcal {S}_{J}}}\ge 0\) satisfying

respectively. Then, taking \({\tilde{\kappa }}\mathop {=}\limits ^\mathrm{{def}}1-{W_{m,q}^{\mathcal {S}_{J}}}{{\overline{\overline{W}}}_{\infty ,q}^{\prime \mathcal {S}_{J}}} {\mathcal {E}_{m,\infty }^{J}} {\mathcal {E}_{\infty ,m}^{J}} \), the uniform bound \(\varvec{{W^{\mathcal {S}_{J}}}}>0\) in (3.1) is given by

The proof of this corollary and each bound are presented in Appendix C. We only note here that \({\tilde{\kappa }}>0\) holds under the same assumptions as in Theorem 3.2. This result follows from the inequality \({\tilde{\kappa }}\ge \kappa \,(>0)\), which is confirmed by the relation \({{\overline{\overline{W}}}_{\infty ,q}^{\prime \mathcal {S}_{J}}} \le \overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}\) in general (see (3.46) in Lemma 3.12 and (C.3) in Lemma C.2 for their explicit definitions).

3.1.1 Computing \({W_{m,q}^{\mathcal {S}_{J}}}\): Rigorous Numerics for the Fundamental Solution Matrix

Here, we start explaining how we get the \({W_{m,q}^{\mathcal {S}_{J}}}\) bound in (3.10). The method presented here is mostly standard for computer-assisted proofs and is closely connected to that dicussed in [30].

From (3.6), the uniformly bounded evolution operator \(\{U(t, s)\}_{0\le s\le t\le \tau }\) corresponds to the fundamental solution matrix \( C^{({\varvec{m}})}(t,s)\) of the finite system (3.4). To construct \(C^{({\varvec{m}})}(t,s)\), we use the relation \(C^{({\varvec{m}})}(t,s) = \Phi (t)\Phi (s)^{-1}\), where \(\Phi (t)\) is the \(m_1\dots m_d\) dimensional matrix-valued function such that

Furthermore, letting \(\Psi (s)\mathop {=}\limits ^\mathrm{{def}}\Phi (s)^{-1}\), \(\Psi (s)\) solves the adjoint problem of the fundamental system (3.29), that is

Using these formulas, the fundamental solution matrix \(C^{({\varvec{m}})}(t,s)\) is explicitly constructed via rigorous numerics, and one can obtain the \({W_{m,q}^{\mathcal {S}_{J}}}\) bound from the definition of \(\bar{U}^{({\varvec{m}})}(t,s)\) in (3.6) such that

Our main task here is to solve the finite differential systems (3.29) for \(t\in [0,\tau ]\) and (3.30) for \(s\in [0,\tau ]\). We present only the case of (3.29), as the other case can be handled by a straight forward extension. Let us define the fundamental solution matrix \(\Phi (t) \in M_{m_1\dots m_d}(\mathbb {R})\) as

From (3.29), each \(c_{{\varvec{k}}}^{({\varvec{j}})}\) (\({\varvec{k}},{\varvec{j}}\in {\varvec{F_m}}\)) solves the IVP of the linear system

In the following, we omit the superscript \(({\varvec{j}})\). By rescaling the time interval \([0,\tau ]\) to \([-1,1]\), and defining as \(G_{{\varvec{k}}}(c)\mathop {=}\limits ^\mathrm{{def}}\frac{\tau }{2}\left( \mu _{{\varvec{k}}}c_{{\varvec{k}}} + \textrm{i}^q ({\varvec{k}}{\varvec{L}})^q\left( D\mathcal {N}(\bar{a})c\right) _{{\varvec{k}}}\right) \), we can rewrite the IVP (3.32) as an integral equation given by

where \(b_{{\varvec{k}}}\) is the given initial data (determined as \(b_{{\varvec{k}}}=(\textrm{e}_{{\textbf {j}} })_{{\varvec{k}}}\) when \({\varvec{j}}\) is fixed). For each \({\varvec{k}}\in {\varvec{F_m}}\), we use the idea presented in [30, 39] and solve the finite number of integral equations (3.33) using Chebyshev series expansions, that is we expand \(c_{\varvec{k}}(t)\) as

where \(c_{-\ell ,{\varvec{k}}} = c_{\ell ,{\varvec{k}}}\) and \(t = \cos (\theta )\). For each \({\varvec{k}}\in {\varvec{F_m}}\), we expand \(G_{\varvec{k}}(c(t))\) using a Chebyshev series, that is

where

and \(N_{\ell ,{\varvec{k}}}(c)\mathop {=}\limits ^\mathrm{{def}}\left( (\tau /2)D\mathcal {N}\left( \bar{a}\right) c\right) _{\ell ,{\varvec{k}}}\) is the Chebyshev coefficients of the Fréchet derivative \(D\mathcal {N}(\bar{a})\) acting on c. Letting \(N_{\varvec{k}}(c) \mathop {=}\limits ^\mathrm{{def}}(N_{\ell ,{\varvec{k}}}(c))_{\ell \ge 0}\), \(\psi _{\varvec{k}}(c) \mathop {=}\limits ^\mathrm{{def}}(\psi _{\ell ,{\varvec{k}}}(c))_{\ell \ge 0}\) and noting that \((\lambda _{\varvec{k}}c_{\varvec{k}})_\ell = \lambda _{\varvec{k}}c_{\ell ,{\varvec{k}}}\), we get that

Combining expansions (3.34) and (3.35) leads to

and this results in solving \(f=0\), where \(f = \left( f_{\ell ,{\varvec{k}}} \right) _{{\ell \ge 0,\,{\varvec{k}}\in {\varvec{F_m}}}}\) is given component-wise by

Hence, for \(\ell >0\) and \({\varvec{k}}\in {\varvec{F_m}}\), we aim at solving

Finally, the problem that we solve is \(f=0\), where \(f = \left( f_{\ell ,{\varvec{k}}} \right) _{\ell ,{\varvec{k}}}\) is given component-wise by

Define the operators (acting on Chebyshev sequences, that is for a fixed \({\varvec{k}}\)) by

Using the operators \(\mathcal {T}\) and \(\Lambda \), we may write more densely for the cases \(\ell >0\) and \({\varvec{k}}\in {\varvec{F_m}}\)

Hence,

Denoting the set of indices \(\mathcal {I}= \{ (\ell ,{\varvec{k}}) \in \mathbb {N}^{d+1}: \ell \ge 0 \text { and } {\varvec{k}}\in {\varvec{F_m}}\}\) and denote \(f=(f_{\varvec{j}})_{\varvec{j}\in \mathcal {I}}\). Assume that using Newton’s method, we computed \({\bar{c}} = ({\bar{c}}_{\ell ,{\varvec{k}}})_{\begin{array}{c} \ell =0,\dots ,n-1 \\ {\varvec{k}}\in {\varvec{F_m}} \end{array}}\) such that \(f({\bar{c}}) \approx 0\). Fix \(\nu _{{C}} \ge 1\) the Chebyshev decay rate and define the weights \(\omega _{\ell ,{\varvec{k}}} =\alpha _{\ell ,{\varvec{k}}} \nu _{{{C}}}^\ell \), where \(\alpha _{\ell ,{\varvec{k}}}\mathop {=}\limits ^\mathrm{{def}}2^{\delta _{\ell ,0}}\alpha _{{\varvec{k}}}\). Given a sequence \(c = (c_{\varvec{j}})_{\varvec{j}\in \mathcal {I}} = (c_{\ell ,{\varvec{k}}})_{\begin{array}{c} {\varvec{k}}\in {\varvec{F_m}}\\ \ell \ge 0 \end{array}}\), denote the Chebyshev-weighed \(\ell ^1\) norm by

The Banach space in which we prove the existence of the solutions of \(f=0\) is given by

The following result is useful to perform the nonlinear analysis when solving \(f=0\) in \(\mathcal {X}^{\varvec{m}}_{{C}}\).

Lemma 3.5

For all \(a,b \in \mathcal {X}^{\varvec{m}}_{{C}}\), \(\Vert a * b\Vert _{\mathcal {X}^{\varvec{m}}_{{C}}} \le \Vert a\Vert _{\mathcal {X}^{\varvec{m}}_{{C}}} \Vert b\Vert _{\mathcal {X}^{\varvec{m}}_{{C}}}\).

The computer-assisted proof of existence of a solution of \(f=0\) relies on showing that a certain Newton-like operator \(c \mapsto c-Af(c)\) has a unique fixed point in the closed ball \(B_r({\bar{c}}) \subset \mathcal {X}^{\varvec{m}}_{{C}}\), where r is a radius to be determined. Let us now define the operator A. Given n, a finite number of Chebyshev coefficients used for the computation of c, denote by \(f^{(n,\varvec{m})}:\mathbb {R}^{nm_1\cdots m_d} \rightarrow \mathbb {R}^{nm_1\cdots m_d}\) the finite dimensional projection used to compute \({\bar{c}} \in \mathbb {R}^{nm_1\cdots m_d}\), that is

First consider \(A^\dagger \) an approximation for the Fréchet derivative \(Df({\bar{c}})\):

where \(Df^{(n,\varvec{m})}({\bar{c}}) \in M_{nm_1\cdots m_d}(\mathbb {R})\) denotes the Jacobian matrix. Consider now a numerical inverse \(A^{(n,\varvec{m})}\) of \(Df^{(n,\varvec{m})}({\bar{c}})\), that is \(\Vert I - A^{(n,\varvec{m})} Df^{(n,\varvec{m})}({\bar{c}}) \Vert \ll 1\). We define the action of A on a vector \(c \in \mathcal {X}^{\varvec{m}}_{{C}}\) as

The following Newton–Kantorovich type theorem (for linear problems posed on Banach spaces) is useful to show the existence of zeros of f (similar formulations can be found for instance in [15, 30, 40]).

Theorem 3.6

Assume that there are constants \(Y_0, Z_0, Z_1 \ge 0\) having that

If

then for all

there exists a unique \({\tilde{c}} \in B_r({\bar{c}})\) such that \(f({\tilde{c}}) = 0\).

Proof

We omit the details of this standard proof. Denote \(\kappa \mathop {=}\limits ^\mathrm{{def}}Z_0 + Z_1 <1\). The idea is show that \(T(c)\mathop {=}\limits ^\mathrm{{def}}c - A f(c)\) satisfies \(T(B_r({\bar{c}})) \subset B_r({\bar{c}})\) and then that \(\Vert T(c_1)-T(c_2) \Vert _{\mathcal {X}^{\varvec{m}}_{{C}}} \le \kappa \Vert c_1-c_2 \Vert _{\mathcal {X}^{\varvec{m}}_{{C}}}\) for all \(c_1,c_2 \in B_r({\bar{c}})\). From the Banach fixed point theorem, there exists a unique \({\tilde{c}} \in B_r({\bar{c}})\) such that \(T({\tilde{c}}) = {\tilde{c}}\). The condition (3.39) implies that \(\Vert \textrm{Id} - A A^\dagger \Vert _{B(\mathcal {X}^{\varvec{m}}_{{C}})} < 1\), and by construction of the operators A and \(A^\dagger \), it can be shown that A is an injective operator. By injectivity of A, we conclude that there exists a unique \({\tilde{c}} \in B_r({\bar{c}})\) such that \(f({\tilde{c}}) = 0\). \(\square \)

Summing up the above, for a given Fourier projection dimension \(\varvec{m}=(m_1, \ldots , m_d)\), we apply Theorem 3.6 to validate the solutions of \(m_1 \cdots m_d\) problems of the form \(f=0\), as given in (3.36), which is equivalent to the IVP of linear system (3.32). This approach yields a sequence of solutions \(c^{({\varvec{j}})}:[0,\tau ] \rightarrow \mathbb {R}^{m_1\cdots m_d}\) (\({\varvec{j}}\in {\varvec{F_m}}\)) with the following Chebyshev series representation:

Thus, we obtain the explicit form of the fundamental solution matrix \(\Phi (t)\) that satisfies (3.29). Using this representation and \(\Psi (s)\), which is the solution to the adjoint problem (3.30), we can derive the \({W_{m,q}^{\mathcal {S}_{J}}}\) bound defined in (3.31) via rigorous evaluation of Chebyshev polynomials with interval arithmetic.

The rest of this section is dedicated to the explicit construction of the bounds \(Y_0\), \(Z_0\) and \(Z_1\) of Theorem 3.6.

The bound \(Y_0\). The first bound will be given by

Using that the term \(N_{\varvec{k}}\) is a finite polynomial and that the numerical solution \({\bar{c}} = ({\bar{c}}_{\ell ,{\varvec{k}}})_{\begin{array}{c} \ell =0,\dots ,n-1 \\ {\varvec{k}}\in {\varvec{F_m}} \end{array}}\), there exists a \(J \in \mathbb {N}\) such that \(f_{\ell ,{\varvec{k}}}({\bar{c}})=0\) for all \(\ell \ge J\). In other words, the term \(f({\bar{c}})=(f_{\varvec{j}}({\bar{c}}))_{\varvec{j}\in \mathcal {I}}\) has only finitely many nonzero terms. Hence, the computation of \(Y_0\) satisfying \(\Vert Af({\bar{c}})\Vert _{\mathcal {X}^{\varvec{m}}_{{C}}} \le Y_0\) is finite and can be rigorously computed with interval arithmetic.

The bound \(Z_0\). The computation of the bound \(Z_0\) satisfying (3.37) requires defining the operator

which action is given by

Using interval arithmetic, compute \(Z_0\) such that

The bound \(Z_1\). For any \(c \in B_1(0)\), let

which is given component-wise by

To simplify the notation of the bound \(Z_1\), lets us first define component-wise uniform bounds \({\hat{z}}_{\ell ,{\varvec{k}}}\) for \( 0 \le \ell < n\) and \({\varvec{k}}\in {\varvec{F_m}}\) such that \(|z_{\ell ,{\varvec{k}}}(\bar{a},c)| \le {\hat{z}}_{\ell ,{\varvec{k}}}\) for all \(c \in \mathcal {X}^{\varvec{m}}_{{C}}\) with \(\Vert c\Vert _{\mathcal {X}^{\varvec{m}}_{{C}}} \le 1\). The first case \(\ell = 0\) can be bounded by

For the case \(\ell = 1,\dots , n-1\), the challenging part to bound is the non-linear term \(N_{\varvec{k}}(c^{c^{(\infty ,\varvec{m} )}})\). First we need to recall that the term \(N_k\) can be express as a finite polynomial in \(\bar{a}\) such that \(N_{\varvec{k}}(c) = \sum _{i = 1}^p \beta _{i} (\bar{a}^i* c)\) for some \(p \in \mathbb {N}\) and \(\beta _i \in \mathbb {R}\) with \(i = 1, \dots , p\). Before bounding this non-linear term, we first look at how to bound the discrete convolution of \(\bar{a}^p\) with \(c^{(\infty , \varvec{m})}\) for some \(p >0\). This convolution can be written as

Using this notation, we see that

where

which can easily be computed using interval arithmetic. For the cases \(1\le \ell <n\) and \({\varvec{k}}\in {\varvec{F_m}}\), this leads to the bound

for all \(c \in B_1(0)\), where \(|\mathcal {T}|\) denotes the operator with component-wise absolute values. Given \(c \in B_1(0)\), we use Lemma 3.5 to conclude that the bound \(Z_1\) is given by

Hence, by construction

satisfies (3.38).

3.1.2 Generation of the Evolution Operator and \(W^{(\infty )}_{{q}}(t,s)\) Bounds

In the followings, unless otherwise noted the index \({\varvec{m}}\) is the one set in Theorem 3.2 and \(\mu _*\) is that in (3.9). Let \(\ell ^1_{\infty }\mathop {=}\limits ^\mathrm{{def}}\left( \textrm{Id}-\varPi ^{({\varvec{m}})}\right) \ell ^1_{\omega } = \left\{ \left( a_{{\varvec{k}}}\right) _{{\varvec{k}}\ge 0}\in \ell ^1_{\omega }:a_{{\varvec{k}}}=0~({\varvec{k}}\in {\varvec{F_m}})\right\} \) endowed with the norm \(\Vert a\Vert _{\ell ^1_{\infty }}\mathop {=}\limits ^\mathrm{{def}}\sum _{{\varvec{k}}\not \in {\varvec{F_m}}}|a_{{\varvec{k}}}|\omega _{{\varvec{k}}}\). Consider a linear operator \(\mathcal {L}_\infty \) whose action is restricted on \(\ell ^1_{\infty }\), that is

where the domain of the operator \(\mathcal {L}_\infty \) is \(D(\mathcal {L}_\infty ) = \left\{ a\in \ell ^1_{\infty }: \mathcal {L}_\infty a\in \ell ^1_{\infty }\right\} \). More explicitly, \(\mathcal {L}_\infty = \left( \textrm{Id}-\varPi ^{({\varvec{m}})}\right) \mathcal {L}\) holds. The operator \(\mathcal {L}_\infty \) generates the semigroup on \(\ell ^1_{\infty }\), which is denoted by \(\left\{ e^{\mathcal {L}_\infty t}\right\} _{t\ge 0}\). Furthermore, the action of such a semigroup \(\left\{ e^{\mathcal {L}_\infty t}\right\} _{t\ge 0}\) can be naturally extended to \(\ell ^1_{\omega }\) by

To get a norm bound of the evolution operator \(\bar{U}^{(\infty )}(t,s)\), let us prepare two lemmas.

Lemma 3.7

Let \(\left\{ e^{\mathcal {L}_\infty t}\right\} _{t\ge 0}\) be the semigroup on \(\ell ^1_{\omega }\) defined in (3.40). For \(\gamma \in (0,1)\), \(\xi \in (0,1]\), and \(\mu _* < 0 \) defined by (3.9), the following estimates:

and

hold for \(\psi \in {\ell _{\omega }^1}\) and \(t\in J\), where

Remark 3.8

The parameter \(\gamma \in (0,1)\) is taken to make the constant \(C_\infty \) in (3.43) bounded. For example, let us set \(\mu _{{\varvec{k}}}=({\varvec{k}}{\varvec{L}})^2(-\varepsilon ^2({\varvec{k}}{\varvec{L}})^2+1)-\sigma \) and \(q=2\), which is the case of our example in Sect. 6. In this case, \(\mu _{{\varvec{k}}}\) is a fourth-order term with respect to \({\varvec{k}}\), whereas \(({\varvec{k}}{\varvec{L}})^q\) is a second-order term. Consequently, it is necessary for \(\gamma \) to be asymptotically greater than or equal to 0.5 to ensure the boundedness of \(C_\infty \).

Proof

Let \(\psi \in {\ell _{\omega }^1}\). It follows for \({\varvec{k}}\notin {\varvec{F_m}}\) from (3.9) that

where we used the inequality

Therefore, it yields from the definition (3.40) that

This directly gives the semigroup estimate (3.42). In the case of \(q=0\), one can directly get the estimate (3.41) by setting \(\gamma =0\), \(\xi =0\), and \(C_\infty =1\). \(\square \)

From Lemma 3.7, we show the existence of the solution of (3.5) with any initial sequence \(c^{(\infty )}(s)=\mathcal {Q}\psi ^{(\infty )}\) in the following theorem. Consequently, the existence of the evolution operator \(\bar{U}^{(\infty )}(t,s)\) is obtained. To state the result, we recall \(\bar{a}(t)\) defined in (2.7) and recall (2.12) that, for each \(t\in J\), there exists a non-decreasing function \(g:(0,\infty )\rightarrow (0,\infty )\) such that \(\left\| D\mathcal {N}(\bar{a}(t))\phi \right\| _{\omega }\le g(\Vert \bar{a}(t) \Vert _{\omega })\Vert \phi \Vert _{\omega }\) for all \(\phi \in \ell _{\omega }^1\).

Theorem 3.9

Let \(\left\{ e^{\mathcal {L}_\infty t}\right\} _{t\ge 0}\) be the semigroup on \(\ell ^1_{\omega }\) defined in (3.40) and assume that Lemma 3.7 holds for \(\gamma \in (0,1)\) and \(\xi \in (0,1]\). There exists a unique solution of the integral equation:

Hence, such a solution \(c^{(\infty )}\) solves the infinite-dimensional system of differential equations (3.5). It yields that the evolution operator \(\bar{U}^{(\infty )}(t,s)\) exists and the following estimate holds

where \(W_{{q}}^{(\infty )}(t,s)\) is defined by

Here, \(\textrm{B}(x,y)\) denotes the Beta function.

Proof

For a fixed \(s> 0\), let us define a map \(\mathcal {P}\) acting on the \(c^{(\infty )}(t)\) as

and let a function space \(\mathcal {X}_{\infty }\) be defined by

Consider \(\beta >1\) and define a distance on \(\mathcal {X}_{\infty }\) as

where H(r) is defined by

Then, \((\mathcal {X}_{\infty },{\textbf{d}})\) is a complete metric space. Let us denote \(g(r)\equiv g(\Vert \bar{a}(r)\Vert _{\omega })\) for the simplicity. We prove that the map \(\mathcal {P}\) becomes a contraction mapping under the distance \({\textbf{d}}\) on \(\mathcal {X}_{\infty }\). For \(c_1^{(\infty )},~c_2^{(\infty )}\in \mathcal {X}_{\infty }\), it follows using (3.42) that

Since \(\beta >1\), the map \(\mathcal {P}\) becomes a contraction mapping on \(\mathcal {X}_{\infty }\). This yields that the solution of (3.5) uniquely exists in \(\mathcal {X}_{\infty }\), which satisfies (3.44).

Moreover, letting

it follows from (3.44) using (3.42) that

where \(G(t)\mathop {=}\limits ^\mathrm{{def}}C_\infty (t-s)^\gamma \). From the Gronwall lemma [41, Chap. 7, p. 356], it follows that

Since G(t) is a non-decreasing function, one gets that

Furthermore,

holds. Here, \(\textrm{B}(x,y)\) is the Beta function.

Then, we conclude that the following inequality holds for any \(\psi \in {\ell _{\omega }^1}\):

where \(W_{{q}}^{(\infty )}(t,s)\) is defined in (3.45). When \(q=0\), an analogous discussion holds by setting \(\gamma =0\), \(\xi =0\), and \(C_\infty =1\). It directly follows (3.45) in the case of \(q=0\).

\(\square \)

3.1.3 The Other Bounds

Here, we will give \({W_{\infty ,q}^{\mathcal {S}_{J}}}\), \({{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\), and \({\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\) bounds satisfying (3.12), (3.13), and (3.14), respectively. Firstly, let us rewrite \(W^{(\infty )}_{0}(t,s)\) in (3.45) as \(e^{\vartheta (t-s)}\), where \(\vartheta =-|\mu _*|+g(\Vert \bar{a}\Vert )\).

Lemma 3.10

Set \(q=0\) and define the constants \({W_{\infty ,0}^{\mathcal {S}_{J}}}>0\), \({{\overline{W}}_{\infty ,0}^{\mathcal {S}_{J}}}\ge 0\), \({\overline{\overline{W}}_{\infty ,0}^{\mathcal {S}_{J}}}\ge 0\) as

respectively. Then \(W^{(\infty )}_{{0}}(t,s)\), defined in (3.45), obeys the inequalities (3.12), (3.13), and (3.14) with \(\gamma =0\).

Proof

First, we note that from (3.45) in the case of \(q=0\)

The inequality (3.12) holds. Second, we note that

and that

These yield that (3.13) holds. Third, note that

Finally, we have

Hence, (3.14) holds. \(\square \)

Remark 3.11

The inequalities (3.13) and (3.14) hold regardless of the positivity or negativity of the variable \(\vartheta \). Furthermore, in the above proof, we used the monotonicity of the functions

with respect to \(t\in J\).

Next, we consider the case of \(q>0\). Let us rewrite \(W^{(\infty )}_{q}(t,s)\) in (3.45) as \(C_\infty (t-s)^{-\gamma }e^{-\iota (t-s) + \vartheta (t-s)^{1-\gamma }}\), where

Similar to Lemma 3.10, we have the following lemma:

Lemma 3.12

Define the constants \({W_{\infty ,q}^{\mathcal {S}_{J}}}>0\), \({{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\ge 0\), \({\overline{\overline{W}}_{\infty ,q}^{\mathcal {S}_{J}}}\ge 0\) as

respectively. Then \(W^{(\infty )}_{q}(t,s)\), defined in (3.45), obeys the inequalities (3.12), (3.13), and (3.14).

Proof

First, it follows from (3.45) that

The inequality (3.12) holds. Second, we note that

and that

These yield that (3.13) holds. Third, note that

Finally, we have

Here we used the formula \(\textrm{B}(x,x) = 2\textrm{B}(x+1,x)\) for \(x>0\) in the last line. Hence, (3.14) holds. \(\square \)

3.2 Validation Theorem for Local Inclusion of the Solution

Lemma 3.13

Recall the definition of \(DF(\bar{a})^{-1}\) in (2.18) and assume that the uniform bounds satisfying (3.1) and (3.2) hold. Let p, \(\psi \in X\) and \(\phi \in \ell _{\omega }^1\). Then \(\mathcal {Q}\psi + p \in Y\) and

Proof

Recalling (2.18), and applying (3.1) and (3.2), it follows that

\(\square \)

Remark 3.14

The parameter \(\gamma \) is taken as \(\gamma =0\) in the case of \(q=0\) in Lemma 3.13.

Theorem 3.15

(Local existence of the solution to IVP) Given the approximate solution \(\bar{a}\in {\mathcal {D}}\) of (2.5), assume that \(\left\| (F(\bar{a}))_1\right\| =\Vert {\dot{\bar{a}}}-\mathcal {L}\bar{a}-\mathcal {Q}\mathcal {N}(\bar{a})\Vert \le \delta \) and \(\Vert (F(\bar{a}))_{2}\Vert _{\omega }=\Vert \varphi -\bar{a}(0)\Vert _{\omega }\le \varepsilon \) hold. Assume also that for \(a_1,a_2\in B_J\left( \bar{a},\varrho \right) \) there exists a non-decreasing function \(L_{\bar{a}}:(0,\infty )\rightarrow (0,\infty )\) such that

Let us assume that there exists \(\varvec{{W^{\mathcal {S}_{J}}}}\), \(\varvec{{W_q^{\mathcal {S}_{J}}}}>0\) satisfying (3.1), (3.2), respectively. Define

If there exists \(\varrho _0>0\) such that

then there exists a unique \(\tilde{a}\in B_J\left( \bar{a},\varrho \right) \) satisfying \(F(\tilde{a})=0\) defined in (2.13). Hence, the solution of the IVP (2.1) exists locally in J.

Remark 3.16

It is worth noting that the non-decreasing function \(L_{\bar{a}}(\varrho )\) depends on the nonlinear term \(\mathcal {N}\). For example, if the nonlinearity is cubic, i.e.,

\(\mathcal {N}_{{\varvec{k}}}(a) = -(a^3)_{{\varvec{k}}}\),

it follows that for \(a_1,a_2\in B_J\left( \bar{a},\varrho \right) \)

Therefore, \(L_{\bar{a}}(\varrho )=3\varrho \left( 2\Vert \bar{a}\Vert +\varrho \right) \) holds in this case.

Proof

We prove that the operator T defined in (2.19) is the contraction mapping on \(B_J(\bar{a},\varrho _0)\) defined in (2.20). Firstly, for any \(a\in B_J(\bar{a},\varrho _0)\), we have using (2.14) and (2.19)

It follows from (3.47) and the assumption of the theorem that

From the assumption of the theorem \(p_{\varepsilon }\left( \varrho _0\right) \le \varrho _0\), \(T(a)\in B_J(\bar{a},\varrho _0)\) holds.

Secondary, we show the contraction property of T on \(B_J(\bar{a},\varrho _0)\). From (2.19) we have for \(a_1, a_2\in B_J(\bar{a},\varrho _0)\)

It follows from (3.47) and the assumption of the theorem that

Hence, \({\frac{\tau ^{1-\gamma }\varvec{{W_q^{\mathcal {S}_{J}}}}}{1-\gamma }L_{\bar{a}}(\varrho _0)}<p_{\varepsilon }(\varrho _0)/\varrho _0\le 1\) holds. Therefore, it is proved that the operator T is the contraction mapping on \(B_J(\bar{a},\varrho _0)\). \(\square \)

The rest of this section is dedicated to the construction of bounds \(\varepsilon \) and \(\delta \) of Theorem 3.15.

The bound \(\varvec{\varepsilon }\). We obtain the \(\varepsilon \) bound of Theorem 3.15. It follows from (2.6) and (2.7) that

Then the \(\varepsilon \) bound is given by

The bound \(\varvec{\delta }\).

We next introduce the \(\delta \) bound of Theorem 3.15. From the definition of the operator F in (2.13) and that of the approximate solution (2.7), we recall

We suppose that the derivative of \(\bar{a}\) with respect to t is denoted by

where the coefficients \(\bar{a}_{\ell ,{\varvec{k}}}^{(1)}\) can be computed by the explicit formula, see, e.g., [42, Sec. 2.4.5],

On the one hand, for \({\varvec{k}}\in {\varvec{F_N}}\), it follows from (3.49) and the fact \(|T_\ell (t)|\le 1\) for all \(\ell \ge 0\) that

where \(\mathfrak {m}_\ell \) denotes the multiplicity for the Chebyshev coefficients

and \(\mathcal {N}_{\ell ,{\varvec{k}}}(\bar{a})\) is the Chebyshev-Fourier coefficients of the nonlinear term \(\mathcal {N}_{{\varvec{k}}}(\bar{a}(t))\) such that

On the other hand, in the case of \({\varvec{k}}\notin {\varvec{F_N}}\), we have from (3.49)

It is worth noting that if we consider the cubic nonlinearity, i.e., \(\mathcal {N}_{{\varvec{k}}}(a) = -(a^3)_{{\varvec{k}}}\), for instance, the above coefficients are given by the discrete convolution for \(d+1\) dimensional tensor.

Furthermore, thanks to the finiteness of the nonzero elements in \(\bar{a}\), these coefficients also have a finite number of nonzero elements. The indices of these nonzero elements are \({\varvec{k}}\in \varvec{F_{3N-2}}\) in the Fourier dimension and \(\ell =0,\dots ,3(n-1)\) in the Chebyshev dimension. This implies that the last term of (3.50) and that of (3.51) are finite sums in this case.

Finally, the defect bound \(\delta \) is given from (3.50) and (3.51) by

In this section, a multi-step scheme is introduced to extend the existence time of the solution to the initial value problem (2.1). This technique is similar to the domain decomposition of the time interval, which is standard for computer-assisted proofs of integrating differential equations presented in [35, 43], etc. By defining the F map within discrete time intervals and numerically verifying the zero-finding problem for the coupled system, it is possible to rigorously include the solution in multiple time intervals. In this context, we consider a sequence of time steps denoted by \(J_i = (t_{i-1}, t_i]\) (\(i=1,2,\dots ,K\)), where \(0 = t_0< t_1< \dots <t_K\) and the step size \(\tau _i\mathop {=}\limits ^\mathrm{{def}}t_i-t_{i-1}\) can be changed. On each time step \(J_i\), we represent the function spaces X, Y, and \(\mathcal {D}\) as \(X_i\), \(Y_i\), and \(\mathcal {D}_i\), respectively. We denote the Banach space

The norm on the space \(\varvec{X}\) is defined for given \(\varvec{a}=(a^{J_1},a^{J_2},\dots ,a^{J_K})\in \varvec{X}\) by

where the superscript is used to indicate the dependence on each time step. Similarly, we use the bold style notations in the same manner for denoting \(\varvec{\mathcal {D}}\mathop {=}\limits ^\mathrm{{def}}\prod _{i=1}^{K} \mathcal {D}_i\), \(\varvec{Y}\mathop {=}\limits ^\mathrm{{def}}\prod _{i=1}^{K} Y_i\), \(\varvec{\ell }_{\omega }^1=\prod _{i=1}^{K} \ell _{\omega }^1\), and \(\varvec{Y}\times \varvec{\ell }_{\omega }^1\mathop {=}\limits ^\mathrm{{def}}\prod _{i=1}^{K} Y_i\times \ell _{\omega }^1\).

The mapping \(F:{\varvec{\mathcal {D}}}\rightarrow \varvec{Y}\times \varvec{\ell }_{\omega }^1\) for the initial value problem over multiple time steps is defined by

where \(F_1:{\mathcal {D}_1}\rightarrow Y_1\times \ell _{\omega }^1\) and \(F_i:{\mathcal {D}_{i-1}}\times {\mathcal {D}_i}\rightarrow Y_i\times \ell _{\omega }^1\) (\(i=2,3,\dots ,K\)) are given by

As we denoted in Sect. 2.2, the first and second components of \(F_i\) (\(i=1,2,\dots , K\)) are also denoted by \((F_i)_1\) and \((F_i)_{2}\), respectively. From the above construction, solving the initial value problem (2.1) over multiple time steps reduces into the zero-finding problem for the map F given in (4.1).

The approximate solution is set by \(\varvec{\bar{a}}=\left( \bar{a}^{J_1},\bar{a}^{J_2},\dots ,\bar{a}^{J_K}\right) \), which is numerically computed on each time step \(J_i\). The Fréchet derivative at the approximate solution \(DF(\varvec{\bar{a}})\) is denoted by

for \(\varvec{h}=(h^{J_1},h^{J_2},\dots ,h^{J_K})\).

To define the Newton-like operator for validating the solution over multiple time steps, we obtain a formula for the inverse of the Fréchet derivative, that is, we solve \(DF(\varvec{\bar{a}})\varvec{h}=(\varvec{p},\varvec{\phi })\) for given \(\varvec{p}=(p^{J_1},p^{J_2},\dots ,p^{J_K})\in \varvec{Y}\) and \(\varvec{\phi }=(\phi ^{J_1},\phi ^{J_2},\dots ,\phi ^{J_K})\). In order to perform this task, we first construct the evolution operators \(\{U_{J_i}(t,s)\}_{t_{i-1}\le s\le t\le t_i}\) (\(i=1,2,\dots ,K\)) validated on each time step \(J_i\) via Theorem 3.2 in Sect. 3.1 and then, denoting \({\mathcal {S}_{J_i}} \mathop {=}\limits ^\mathrm{{def}}\{(t, s): t_{i-1}< s < t\le t_{i}\}\) \((i=1,2,\dots ,K)\), we obtain positive constants \(\varvec{{W^{\mathcal {S}_{J_i}}}}\), \(\varvec{{W_q^{\mathcal {S}_{J_i}}}}\), \(\varvec{{W^{(t_i,J_i)}}}\), \(\varvec{{W^{(t_i,J_i)}_q}}\), \(\varvec{{W^{(t_i,t_{i-1})}}}\) such that

respectively. Note that these constants fall into three types: First, \(\varvec{{W^{\mathcal {S}_{J_i}}}}\) and \(\varvec{{W_q^{\mathcal {S}_{J_i}}}}\) are bounds for the two parameters (t, s), with the superscript \(\mathcal {S}_{J_i}\) indicating the dependency \((t,s)\in \mathcal {S}_{J_i}\). Second, \(\varvec{{W^{(t_i,J_i)}}}\) and \(\varvec{{W^{(t_i,J_i)}_q}}\) are bounds for the one parameter s, where the superscript \((t_i,J_i)\) indicates that \(t=t_i\) is fixed and \(s\in J_i\). Third, \(\varvec{{W^{(t_i,t_{i-1})}}}\) denotes the operator norm at the specific times \(t=t_i\) and \(s=t_{i-1}\).

Remark 4.1

Note that computing each evolution operator \(U_{J_i}(t,s)\) by solving the linearized problems and then obtaining the bounds in (4.3) at each step \(i=1,\dots ,K\) can be done independently. Hence each computation can be made in parallel, which implies that the computational cost is additive in time rather than multiplicative.

Having obtained the evolution operators \(U_{J_i}(t,s)\) on each \(J_i\), the solution \(\varvec{h}\) of \(DF(\varvec{\bar{a}})\varvec{h}=(\varvec{p},\varvec{\phi })\) is derived iteratively as follows: The first step is the same as that presented in Sect. 2.3. The variation of constants formula yields

Next, the \(h^{J_2}\) is given by

Here, we denote \(U_{J_1\cup J_2}(t,t_0)\mathop {=}\limits ^\mathrm{{def}}U_{J_2}(t,t_1)U_{J_1}(t_1,t_0)\) (\(t\in J_2\)). We iterate this process for \(i=1,2,\dots ,K\). The resulting \(h^{J_K}\) is obtained for \(t\in J_K\) by

where

Finally, we get an unique expression of the inverse of \(DF(\varvec{\bar{a}})\) as

The Newton-like operator is defined by

Therefore, we validate the solution of \(F(\varvec{a})=0\) given in (4.1) to show that \(T:\varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }) \rightarrow \varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho })\) is a contraction mapping, for some \(\varvec{\varrho }>0\), where

Lemma 4.2

For each \(i=1,2,\ldots ,K\), consider the positive constants \(\varvec{{W^{\mathcal {S}_{J_i}}}}\), \(\varvec{{W_q^{\mathcal {S}_{J_i}}}}\), \(\varvec{{W^{(t_i,J_i)}}}\), \(\varvec{{W^{(t_i,J_i)}_q}}\), \(\varvec{{W^{(t_i,t_{i-1})}}}\) satisfying (4.3), and let \(\varvec{W_\ell }>0\), \(\varvec{W_{X}}>0\) and \(\varvec{W_{Y}}>0\) be defined by the matrix infinity norm of the following lower triangular matrices:

where \({\hat{\tau }}_i\mathop {=}\limits ^\mathrm{{def}}\frac{\tau _i^{1-\gamma }}{1-\gamma }\) (\(i=1,2,\dots ,K\)), and

respectively. Here, we denote \(\varvec{W^{(t_i,t_j)}}=\varvec{W^{(t_i,t_{i-1})}}\cdots \varvec{W^{(t_{j+1},t_j)}}\) with \(0\le j<i\le K-1\). Then, these constants obey the following inequality:

where \(\Vert \varvec{\phi }\Vert _{\varvec{\omega }}\mathop {=}\limits ^\mathrm{{def}}\max \{\Vert \phi ^{J_i}\Vert _\omega \}\).

Proof

For given \(\varvec{\phi }=(\phi ^{J_1},\phi ^{J_2},\dots ,\phi ^{J_K})\in \varvec{\ell _{\omega }^1}\), \(\varvec{\psi }=(\psi ^{J_1},\psi ^{J_2},\ldots ,\psi ^{J_K})\in \varvec{X}\) and \(\varvec{p}=(p^{J_1},p^{J_2},\dots ,p^{J_K})\in \varvec{X}\), taking \(X_i\) norm in each component of (4.6), we have from (3.47)

where \({\hat{\tau }}_i=\frac{\tau _i^{1-\gamma }}{1-\gamma }\) (\(i=1,2,\dots ,K\)). We also have

where we used

Similarly, we obtain the K-th component

where \(\varvec{W^{(t_i,t_j)}}=\varvec{W^{(t_i,t_{i-1})}}\cdots \varvec{W^{(t_{j+1},t_j)}}\) (\(0\le j<i\le K-1\)). Summing up the above inequalities and taking the maximum norm of each components, the \(\varvec{X}\) norm of (4.6) is given by

\(\square \)

Theorem 4.3

(Local existence over multiple time steps) Given the approximate solution \(\varvec{\bar{a}}=\left( \bar{a}^{J_1},\bar{a}^{J_2},\dots ,\bar{a}^{J_K}\right) \in {\varvec{\mathcal {D}}}\), assume that

and

Assume also that for \(a_1^{J_i},a_2^{J_i}\in B_{J_i}\left( \bar{a}^{J_i},\varvec{\varrho }\right) \) \((i=1,2,\dots ,K)\) there exists a non-decreasing function \(L_{\bar{a}^{J_i}}:(0,\infty )\rightarrow (0,\infty )\) such that

Set \(\varvec{\delta }=\max \{\delta ^{J_i}\}\), \(\varvec{\varepsilon }=\max \{\varepsilon ^{J_i}\}\) and define

where positive constants \(\varvec{W_\ell }\), \(\varvec{W_X}\), and \(\varvec{W_Y}\) are defined in Lemma 4.2. If there exists \(\varvec{\varrho }_0>0\) such that

then there exists a unique \(\varvec{\tilde{a}}\in \varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\) satisfying \(F(\varvec{\tilde{a}})=0\) defined in (4.1). Hence, the solution of the IVP (2.1) exists locally over the multiple time steps.

Proof

We prove that the operator T defined in (4.7) is the contraction mapping on \(\varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\) in (4.8). Firstly, for any \(\varvec{a}\in \varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\), we have using (4.2) and (4.7)

It follows using Lemma 4.2 and the assumption of the theorem that

From the assumption of the theorem \(p_{\varvec{\varepsilon }}(\varvec{\varrho }_0)\le \varvec{\varrho }_0\), then \(T(\varvec{a})\in \varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\) holds.

Secondary, we show the contraction property of T on \(\varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\). From (4.7) we have for \(\varvec{a}_1, \varvec{a}_2\in \varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\)

Using Lemma 4.2 and the assumption of the theorem, we obtain

Hence, \(\varvec{W_X}\max _{i}\left\{ L_{\bar{a}^{J_i}}(\varvec{\varrho }_0)\right\} <p_{\varvec{\varepsilon }}(\varvec{\varrho }_0)/\varvec{\varrho }_0\le 1\) holds. Therefore, it is proved that the operator T is the contraction mapping on \(\varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\). \(\square \)

4.1 Bounds on Each Time Step

The remaining bounds to complete the proof of the local existence over the multiple time steps are positive constants \(\varvec{{W^{(t_i,J_i)}}}\), \(\varvec{{W^{(t_i,J_i)}_q}}\), and \(\varvec{{W^{(t_i,t_{i-1})}}}\) such that \(\sup _{s\in J_i}\Vert U_{J_i}(t_i,s)\Vert _{B(\ell _{\omega }^1)}\le \varvec{{W^{(t_i,J_i)}}}\), \(\sup _{s\in [t_{i-1},t_i)}(t_i-s)^\gamma \Vert U_{J_i}(t_i,s)\mathcal {Q}\Vert _{B(\ell _{\omega }^1)}\le \varvec{{W^{(t_i,J_i)}_q}}\), and \(\Vert U_{J_i}(t_i,t_{i-1})\Vert _{B(\ell _{\omega }^1)}\le \varvec{{W^{(t_i,t_{i-1})}}}\) (\(i=1,2\dots ,K\)), respectively.

To get these bounds, we consider the linearized problem (2.17) and represent the solution of (2.17) as \(b_s(t)\equiv b(t)\) (\((t,s)\in {\mathcal {S}_{J_i}}\)). Since we have a rigorous inclusion of the fundamental solution in Sect. 3.1.1, one can obtain two bounds \({W_{m,q}^{(t_i,J_i)}}>0\) and \( {W_{m}^{(t_i,t_{i-1})}}>0\) from the definition of \(\bar{U}^{({\varvec{m}})}(t,s)\) in (3.6) such that

and

respectively.

Returning to consider \(b_s(t)\) in Sect. 3.1, we recall the bounds (3.21) and (3.22). We define a matrix whose entries are the coefficients of the bounds as

For any \(s\in J_i\), we have

where \({\varvec{U}_{kj}^{J_i}}\) denotes the (k, j) entry of \({\varvec{U}^{J_i}}\). Similarly, we also define another matrix to obtain the \(\varvec{{W^{\mathcal {S}_{J_i}}}}\) bound in (3.28) as

We also represent the norm bound \(b_s(t)\) as

with the same manner.

Using these bounds, since the solution \(b_s(t)\) of (2.17) satisfies the integral equations in (3.8), the finite part by the Fourier projection is bounded using the \({W_{m,q}^{(t_i,J_i)}}\) bound in (4.10), (3.13), (3.14), (3.15), and (3.20) by

On the other hand, the infinite part is bounded using (3.12), (3.13), (3.14), (3.16), and (3.19) by

Consequently, the \(\varvec{{W^{(t_i,J_i)}_q}}\) bound is given by

Moreover, an analogous discussion using the bounds in Corollary 3.4 directly yields the \(\varvec{{W^{(t_i,J_i)}}}\) bound given by

Next, setting \(s=t_{i-1}\) in (3.8) and the initial data \(\phi \in \ell _{\omega }^1\) instead of \(\mathcal {Q}\psi \), the finite part \(b^{({\varvec{m}})}_{t_{i-1}}(t_i)\) is bounded using the bounds \( {W_{m}^{(t_i,t_{i-1})}}\) and \({W_{m,q}^{(t_i,J_i)}}\), given in (4.11) and (4.10) respectively, (3.15), (3.20), (3.25), and (3.27) by

The \(\ell _{\omega }^1\) norm bound of the infinite part \(b^{(\infty )}_{t_{i-1}}(t_i)\) is derived using \(W^{(\infty )}(t,s)\) in (3.23), (3.16), (3.19), (3.26), and (3.27) by

The \(\varvec{{W^{(t_i,t_{i-1})}}}>0\) bound satisfying \(\Vert U_{J_i}(t_i,t_{i-1})\Vert _{B(\ell _{\omega }^1)}\le \varvec{{W^{(t_i,t_{i-1})}}}\) is finally given by

4.2 Error Bounds at the End Point

When the local existence of the solution is proved, we obtain error bounds at the end point of \(J_i\) as follows: Firstly, setting \(t=t_1\) in (3.48), it follows from Lemma 3.13

where \(\varvec{\varrho }_0>0\) is the validated radius in Theorem 4.3 of the ball defined in (4.8). Secondly, setting \(t=t_2\) in (4.9), the exact form of the inverse (4.4) yields

Similarly, we have an error bound at \(t=t_K\) from (4.5) and (4.9)

Finally, the \(\varepsilon ^{J_i}\) bound satisfying

is a numerical error between two different approximate solutions, which is caused by the Chebyshev interpolation. This is tiny error in general and easily computed by the form

where \(n_i\) denotes the size of Chebyshev coefficients on \(J_i\).

In the following sections, we demonstrate some applications of our rigorous integrator. Firstly, in this section, we provide a computational approach of proving the existence of global in time solutions. We then apply the provided method for computer-assisted proofs of global existence of solutions to initial value problems of the 3D/2D Swift–Hohenberg equation. We combine the rigorous forward integration with explicit constructions of trapping regions to prove global existence converging to asymptotically stable nontrivial equilibria. Secondly, in the next section, the 2D Ohta–Kawasaki equation is considered to deal with derivatives of the nonlinear term. All computations in this study were conducted on Windows 10, Intel(R) Core(TM) i9-10900K CPU @ 3.70GHz using MATLAB R2022a with INTLAB - INTerval LABoratory [44] version 11, which supports interval arithmetic. The code used to generate the results presented in the following sections is freely available from [45].

5.1 Constructing Trapping Region

We propose a strategy to demonstrate global existence of (1.1) via the mechanism of convergence towards an asymptotically stable equilibrium solution. We are specifically looking at cases where the nonlinearity of (1.1) is of the form (2.3) with \(q = 0\). We denote the vector field of the infinite-dimensional system of ODEs (2.4) by

Assume the existence of an equilibrium solution \(\tilde{a}\in \ell _{\omega }^1\), that is

and assume that the Fréchet derivative \(Df(\tilde{a})\) is an invertible linear operator with real eigenvalues. In addition, let \(\mathcal {G}:\ell _{\omega }^1 \rightarrow \ell _{\omega }^1\) be defined by

so that \(\mathcal {G}(0)=f(\tilde{a})=0\) and \(D\mathcal {G}(0) = Df(\tilde{a}) - Df(\tilde{a})=0\) hold.

The following theorem asserts existence of the trapping region associated with the asymptotically stable equilibrium solution. Note that a similar argument was consider in [46] to construct an attracting neighborhood in Fisher’s equation, which was combined with a computer-assisted proof to prove global existence.

Theorem 5.1

(Global existence via asymptotic convergence) Let \(\tilde{a}\in \ell _{\omega }^1\) be an equilibrium solution of (5.1). Assume that there exist \(C\ge 1\) and \(\lambda >0\) such that

for all \(t>0\). Consider \(\epsilon >0\) small enough so that

Since \(\mathcal {G}(0)=0\) and \(D\mathcal {G}(0)=0\), pick \(\rho >0\) (as large as possible) such that

If

then for all \(t > 0\)

More explicitly, if \(a(0) \in B_{\frac{\rho }{C}}(\tilde{a}) \subset \ell _{\omega }^1\), then the solution a(t) is bounded (it stays in \(B_{\frac{\rho }{C}}(\tilde{a})\)), it exists globally in time and it satisfies

Proof

While the proof is standard in the theory of differential equations, we present it here since it provides a computational procedure to show global existence and asymptotic convergence.

Consider \(t>0\) and let \(h(t) \mathop {=}\limits ^\mathrm{{def}}a(t) - \tilde{a}\). From (5.2), one gets that

with \(\mathcal {G}(0)=0\) and \(D\mathcal {G}(0)=0\). The variation of constants formula in \(h \in \ell _{\omega }^1\) for (5.6) yields

and therefore

Combining (5.3) and (5.7), one obtains

Now, as long as \(\Vert h(s)\Vert _{\omega } \le \rho \) for \(s \in [0,t]\), then combining assumption (5.5) and (5.8) yields

Letting \(y(t) \mathop {=}\limits ^\mathrm{{def}}e^{\lambda t} \Vert h(t)\Vert _{\omega } \ge 0\), one obtains from the last inequality that

Applying Gronwall’s inequality to (5.9)

which implies that

Observe that if \(\Vert h(0)\Vert _{\omega } \le \rho /C\), then (5.10) guarantees that \(\Vert h(t)\Vert _{\omega } \le \rho e^{-\epsilon t} < \rho \) for all \(t > 0\). Going back to the original coordinates \(a(t) = h(t) + \tilde{a}\), then we conclude that

\(\square \)

Next, we provide an explicit procedure which will be used to conclude, via a successful application of Theorem 5.1, that a solution exists globally in time and that it converges to an asymptotically stable equilibrium solution.

Procedure 5.2

To study the asymptotic behaviour of solutions close to an asymptotically stable equilibrium solution, perform the following steps rigorously:

(P1):

Compute \(\tilde{a}\in \ell _{\omega }^1\) such that \(f(\tilde{a})=0\);

(P2):

Define \(\mathcal {G}\) from (5.2);

(P3):

Compute \(C>0\) and \(\lambda >0\) such that (5.3) holds;

(P4):

Consider \(\epsilon >0\) small enough (in fact the smaller the better) so that (5.4) holds;

(P5):

Let \(\delta \mathop {=}\limits ^\mathrm{{def}}\lambda - \epsilon > 0\);

(P6):

Compute \(\rho =\rho (\delta ,C)>0\) (as large as possible) such that (5.5) holds;

(P7):

As we integrate equation (5.1), verify that \(\Vert a(\tau ) - \tilde{a}\Vert _{\omega } \le \frac{\rho }{C}\), for some \(\tau \ge 0\);

(P8):

From Theorem 5.1, conclude that a(t) stays in \(B_{\rho }(\tilde{a})\) for all \(t \ge \tau \) (it is bounded), it exists globally in time and it satisfies \(a(t) \rightarrow \tilde{a}\) as \(t \rightarrow +\infty \).

If Procedure 5.2 is successfully applied, then as we perform the rigorous integration, we verify if \(a(\tau ) \in B_{\frac{\rho }{C}}(\tilde{a})\). If so, then the solution is trapped in the ball \(B_{\rho }(\tilde{a})\) for all time \(t \ge \tau \). For this reason, we say that the set \(B_{\rho }(\tilde{a})\) is a trapping region.

5.2 Semigroup Estimates

We will show how one get the constants \(C\ge 1\) and \(\lambda >0\) satisfying

This corresponds to (P3) in Procedure 5.2. Since the equilibrium solution is asymptotically stable, such \(\lambda >0\) is expected to always exit. Before presenting the computation of the constants, we first need to construct an operator P, for which the column will be an approximation of the eigenvectors of \(Df(\tilde{a})\). Let \({\bar{\lambda }}_{{\varvec{k}}}\) (\({\varvec{k}}\in {\varvec{F_M}}\)) be approximate eigenvalues of the truncated matrix of \(Df(\tilde{a})\) with the size \({\varvec{M}}=(M_1,\dots ,M_d)\).^{Footnote 1} These eigenvalues can be sorted in increasing order with respect to \({\varvec{k}}\). To present a decomposition of an infinite dimensional operator, say \(\mathcal {A}:\ell _{\omega }^1\rightarrow \ell _{\omega }^1\) for example, into finite and tail parts, we denote a finite truncation of the size \({\varvec{M}}\) by \(\mathcal {A}_1:\ell _{\omega }^1\rightarrow \mathbb {C}^{N}\), where \(N=\prod _{i=1}^dM_i\), and the tail part by \(\mathcal {A}_\infty :\ell _{\omega }^1\rightarrow \ell _{\infty }^1=(\textrm{Id}-\varPi ^{({\varvec{M}})})\ell _{\omega }^1\), respectively. Let \(P_1\) be an invertible matrix whose columns are numerical approximation of the eigenvectors corresponding to \({\bar{\lambda }}_{{\varvec{k}}}\). Then we set an operator P and its inverse \(P^{-1}\) as

Therefore, using these operators, we define the operator \(M\mathop {=}\limits ^\mathrm{{def}}P^{-1}Df(\tilde{a})P\). Since P is an approximation of the eigenvectors of \(Df(\tilde{a})\), and \(P^{-1}\) is the exact inverse of P, we can see M as being an approximate standard diagonalization of \(Df(\tilde{a})\). For this reason we will refer to M as the pseudo-diagonalization of \(Df(\tilde{a})\). The motivation for the pseudo-diagonalization is presented later in this section. We can rewrite M with the form

where \(\Lambda _1=\textrm{diag}({\bar{\lambda }}_{{\varvec{k}}})\) (\({\varvec{k}}\in {\varvec{F_M}}\)) and \(\Lambda _\infty =\textrm{diag}(\mu _{{\varvec{k}}})\) (\({\varvec{k}}\notin {\varvec{F_M}}\)) are the finite/infinite diagonal matrices, \({\mathscr {E}}_1^1\in \mathbb {C}^{N\times N}\) contains numerical errors of diagonalization, \({\mathscr {E}}_1^\infty :\ell _{\infty }^1\rightarrow \mathbb {C}^{N}\), \({\mathscr {E}}_\infty ^1:\mathbb {C}^{N}\rightarrow \ell _{\infty }^1\), and \({\mathscr {E}}_\infty ^\infty :\ell _{\infty }^1\rightarrow \ell _{\infty }^1\). More precisely,

Using this notation, we get the bound

By construction of the operator P and \(P^{-1}\), the computation of \(\Vert P \Vert _{B(\ell _\omega ^1)}\) and \(\Vert P^{-1} \Vert _{B(\ell _\omega ^1)}\) are finite and given by

The next step is to define the necessary conditions to use the following theorem.

Theorem 5.3

[47, Theorem B.1] Consider the linear operators M, \(\Lambda \), \({\mathscr {E}}: \mathbb {C}^{N} \times \ell _{\infty }^1\rightarrow \mathbb {C}^{N} \times \ell _{\infty }^1\) defined in (5.11). We require \(\Lambda \) to be densely defined and \({\mathscr {E}}\) to be bounded. Suppose that \(\Lambda _1\) is diagonal and that \(\Lambda _\infty :\ell _{\infty }^1\rightarrow \ell _{\infty }^1\) has a bounded inverse.

Let \(\mu _1\), \(\mu _\infty > 0 \) such that for \(t>0\) and \((\phi _1,\phi _\infty )\in \mathbb {C}^{N}\times \ell _{\infty }^1\) we have

Fix constants \(\delta _a\), \(\delta _b\), \(\delta _c\), \(\delta _d>0\) satisfying

and set

Assume the following inequalities hold

Then, we have the following estimates of semigroup generated by M:

where

Computing \(C_s\) and \(\lambda _s\) using interval arithmetic and combining them with (5.12) give us

where

It is worth mentioning that if we applied Theorem 5.3 directly to the operator \( \Vert e^{Df(\tilde{a}) t} \Vert _{B(\ell _\omega ^1)}\) without pseudo-diagonalizing first, we would have \(\lambda < 0\), and thus we would not have the hypothesis needed for Theorem 5.1 to prove global existence. Also, the truncation size \({\varvec{M}}\) of the finite part of the pseudo-diagonalization is chosen such that the conditions (5.13) are satisfied and such that \(\lambda > 0\).

5.3 Computer-Assisted Proofs

Let us apply Procedure 5.2 to study global existence and asymptotic convergence of solutions in the Swift–Hohenberg equation, which is introduced in [48] to describe the onset of Rayleigh–Bénard convection

where \(\sigma \in \mathbb {R}\) is a parameter. It is worth mentioning that several computer-assisted proofs have been presented in the last twenty years to study the dynamics in (5.14), including constructive proofs of steady states [40, 49,50,51,52], global dynamics [53], chaos [54, 55], stable manifold of equilibria [47] and solutions to initial value problems in the one-dimensional case [34, 35].

Here, we impose the homogeneous Neumann boundary condition on 3D prism domain \([0,\pi /L_1] \times [0,\pi /L_2]\times [0,\pi /L_3]\) with \((L_1,L_2,L_3) = (1,1.1,1.2)\) (resp. 2D rectangle domain \([0,\pi /L_1] \times [0,\pi /L_2]\) with \((L_1,L_2)=(1,1.1)\)). Using the general notation (1.1), \(\lambda _0=\sigma -1\), \(\lambda _1=-2\), \(\lambda _2=-1\) and \({\Delta ^p N(u)} = -u^3\) (\(p=0\)). In this case, the corresponding differential equation (5.1) in the space of Fourier coefficients in \(\ell _{\omega }^1\) is given by

where \((\mathcal {L}a)_{\varvec{k}}= (\sigma - ( 1- ({\varvec{k}}{\varvec{L}})^2)^2)a_{{\varvec{k}}} (=\mu _{{\varvec{k}}}a_{{\varvec{k}}})\) and the cubic term is the usual convolution defined by

For the part (P1) of Procedure 5.2, we use the tools of computer-assisted proofs for equilibria of PDEs defined on 3D/2D domains with periodic boundary conditions (e.g. see [40, 51, 56]) to compute \(\tilde{a}\in \ell _{\omega }^1\) such that \(f(\tilde{a}) = \mathcal {L}\tilde{a}- \tilde{a}^3\). A brief introduction of such tools is presented in Appendix D. For the part (P2), we note that for the Swift–Hohenberg equation,

For the part (P3), the explicit construction of the constants \(C>0\) and \(\lambda >0\) satisfying (5.3) is performed rigorously via Theorem 5.3 in Sect. 5.2. For the part (P4), fix \(\epsilon = 10^{-16}\), which is small enough so that \(\delta \mathop {=}\limits ^\mathrm{{def}}\lambda -\epsilon >0\) (the part (P5)). For the part (P6), one must pick \(\rho =\rho (\delta ,C)>0\) such that (5.5) hold. In our context, we use the Banach algebra structure of \(\ell _{\omega }^1\) to reduce the verification of (P6) to

for all \(\phi \in \ell _{\omega }^1\) with \(\Vert \phi \Vert _{\omega } \le \rho \). Assuming \(\phi \ne 0\), the previous inequality boils down to verify that

for all \(\phi \in \ell _{\omega }^1\) \(\Vert \phi \Vert _{\omega } \le \rho \). A sufficient condition for this to hold is that

The largest positive \(\rho \) which satisfies this inequality is given by

For the part (P7), we integrate (5.1) via the multi-step scheme introduced in Sect. 4 to prove that \(\Vert a(t_K) - \tilde{a}\Vert _{\omega } \le \frac{\rho }{C}\) holds for some \(t_K \ge 0\). To perform this, we denote \(\mathcal {N}_{{\varvec{k}}}(a)=-(a^3)_{{\varvec{k}}}\) in (2.4). The Fréchet derivative of \(\mathcal {N}\) at \(\bar{a}\) is given by \(D\mathcal {N}(\bar{a})\phi =-3\left( \bar{a}^2*\phi \right) =-3(\bar{a}*\bar{a}*\phi )\) for \(\phi \in \ell _{\omega }^1\) because of the fact

Then the function g satisfying (2.12) is given by \(g(\Vert \psi \Vert _{\omega }) = 3\Vert \psi \Vert _{\omega }^2\) for \(\psi \in \ell _{\omega }^1\). If we take \(\psi =\bar{a}(t)\) for a fixed t, this g can also be expressed as \(g(\Vert \bar{a}(t)\Vert _{\omega })=3\Vert \bar{a}^2(t)\Vert _{\omega }\). This is because that \(\bar{a}^2(t)\) is finite-dimensional sequence in the Fourier dimensions. In addition, \(\bar{a}(t)\in \ell _{\omega ^{-1}}^\infty \) holds from its finiteness. From Lemma B.2, if we set \(c=c^{(\infty )}\in \ell _{\omega }^1\) and \(a=\bar{a}^2(t)\in \ell _{\omega ^{-1}}^\infty \) in (B.2), it follows that

where \({\varvec{N}}\) is the size of the Fourier dimension of \(\bar{a}\). Let

Using this notation, it follows from (3.3) that

Therefore, the constant \({\mathcal {E}_{m,\infty }^{J}}\) satisfying (3.15) is given by

where \(\Psi _{\ell ,{\varvec{k}}}\left( \bar{a}^2\right) \) denotes the Chebyshev coefficients of \(\Psi _{{\varvec{k}}}\left( \bar{a}^2(t)\right) \). Similarly, we have for \({\varvec{k}}\notin {\varvec{F_m}}\)

To get the \({\mathcal {E}_{\infty ,m}^{J}}\) bound satisfying (3.16), we have from the Banach algebra property

Hence, we obtain \({\mathcal {E}_{\infty ,m}^{J}}\) in (3.16) as

where \(\left( \bar{a}^2\right) _{\ell ,{\varvec{k}}}\) also denotes the Chebyshev coefficients of \(\left( \bar{a}^2(t)\right) _{{\varvec{k}}}\). Finally, recalling Remark 3.16, we have \(L_{\bar{a}}(\varrho )=3\varrho \left( 2\Vert \bar{a}\Vert _X+\varrho \right) \) in Theorem 3.15.

Our integrator based on Theorem 4.3 proves that there exists a solution of Swift–Hohenberg equation (5.14) in \(\varvec{B}_K(\varvec{\bar{a}}, \varvec{\varrho }_0)\) defined in (4.8). Then we rigorously compute the error bound at the end point, say \(\varvec{\varrho }_K\), by (4.12). We verify that the solution is in the trapping region via the following inequality:

The second term in the last inequality is rigorously computable via the result of computer-assisted proofs for the equiliubrium based on interval arithmetic. Consequently, if the last inequality holds, then we have a computer-assisted proof of global existence of the solution to Swift–Hohenberg equation (5.14) for the part (P8). We note that, while the global existence of solutions of the Swift–Hohenberg equation is already established, our numerical integrator enables us to begin with any initial condition far from equilibrium, without prior knowledge of the equilibrium to which it will converge. Moreover, our method is not restricted to gradient systems and could potentially be applied to more complex problems.

5.3.1 3D Swift–Hohenberg Equation

From now on, we fix \(\sigma =0.04\) in the 3D case of (5.14). We set \({\varvec{N}}=(4,4,4)\) for the size of Fourier and \({\varvec{m}}=(2,1,1)\) for the Fourier projection to obtain the solution map shown in Sect. 3.1. We also set \(\nu _{{F}}=1\) for the \(\ell _{\omega }^1\) norm. As shown in Fig. 1f, we have a stable (stripe pattern) equilibrium of Swift–Hohenberg equation (5.14) via the tools of computer-assisted proofs. The explicit construction of the constants C and \(\lambda \) introduced in Sect. 5.2 yields that \(C=1.0216\) and \(\lambda =0.0799\) for the part (P3). The largest positive \(\rho \) defined in Sect. 5.3 is obtained as \(\rho = 0.0988\) (the part (P6)). Therefore, our target neighborhood of the equilibrium is \(\rho /C = 0.0967\).

The profiles of solution to the 3D Swift–Hohenberg: From the initial data (a), our integrator proves the existence of solution locally in time. After i (\(i=25\), 100, 150, \(\dots \)) steps, the time evolution of solution profile is fully changed in the 3D prism domain. Consequently, after 287 steps, we verified that the exact solution enclosure is in the trapping region in the part (P7) of Procedure 5.2. Therefore, we proved that there exists a global in time solution from the initial data (a) to the stripe pattern equilibrium (f)

Full size image

The initial data is set as

with

Setting \(\tau _i=0.25\), our rigorous integrator proves that the solution of the IVP of (5.14) exists at least to \(t_K=71.75\) (\(K=287\)). The resulting error bound \(\varvec{\varrho }_0\) via Theorem 4.3 is given by

In this case we also obtain \(\varvec{\varrho }_K=4.7241\cdot 10^{-6}\). The value of \(\varvec{\varrho }_K + \Vert \bar{a}^{J_K}(t_K) - \tilde{a}\Vert _{\omega }\) is bounded by 0.0962, which is less than \(\rho /C\). Hence, for the part (P8), the proof of global existence of the solution to the Swift–Hohenberg equation (5.14) is completed.

Several profiles of the solution are presented in Fig. 1, demonstrating that the evolutionary behavior of the solution undergoes significant changes within the 3D prism domain. We believe that such a capability to visualize the solution’s evolution in a higher domain, underpinned by mathematical proof, provides a new perspective of pattern dynamics.

5.3.2 2D Swift–Hohenberg Equation

Next let us consider the 2D case of Swift–Hohenberg equation (5.14). The parameter is fixed as \(\sigma = 3\). Computational parameters are set by \({\varvec{N}}=(12,12)\), \({\varvec{m}}=(3,3)\), \(\nu _{{F}}=1.0\). Our tool of computer-assisted proofs obtains two stable equilibria of Swift–Hohenberg equation (5.14) shown in Fig. 2. Both equilibria are asymptotically stable. We apply the provided approach of computer-assisted proofs for existence of global in time solutions.

Profiles of two stable equilibria of the 2D Swift–Hohenberg equation: a Stripe pattern equilibrium solution \(\tilde{u}_1(x)\) and b spot pattern equilibrium solution \(\tilde{u}_2(x)\)

Full size image

Stripe pattern equilibrium. The constants C and \(\lambda \) are determined by \(C=2.7295\) and \(\lambda =1.9406\), respectively. The largest positive value of \(\rho \) introduced in Sect. 5.3 is identified as \(\rho = 0.1138\). This leads the radius of our target neighborhood near the equilibrium solution, given as \(\rho /C = 0.0416\). The initial data \(u_0(x)\) we used here is plotted in Fig. 3a. Using the provided integrator, we achieve a rigorous inclusion of the solution for the IVP up to \(t_K=1.4141\) (where \(K=181\)), by setting \(\tau _i=7.8125\cdot 10^{-3}\) (\(i=1,2,\ldots ,K\)) uniformly for each step.

The resulting error bound \(\varvec{\varrho }_0\) is given by

The \(\varvec{\varrho }_K\) bound is obtained by \(\varvec{\varrho }_K=2.6919\cdot 10^{-6}\). Recalling the part (P7) of Procedure 5.2, the value of \(\varvec{\varrho }_K + \Vert \bar{a}^{J_K}(t_K) - \tilde{a}\Vert _{\omega }\) is bounded by 0.0414. This is definitely less than \(\rho /C\). Therefore, we have a computer-assisted proof of global existence of the solution asymptotically converging to the stripe pettern equilibrium solution plotted in Fig. 2a.

The initial data \(u_0(x)\) and the difference between the stable stripe pattern equilibrium solution \(\tilde{u}_1(x)\) and \(u_0(x)\): to demonstrate computer-assisted proofs for global existence beyond the neighborhood of the equilibrium solution, our integrator is essential for capturing variations in the solution profile

Full size image

As shown in Fig. 3b, there is a substantial difference between the equilibrium solution \(\tilde{u}_1(x)\) and the initial data \(u_0(x)\) in this case. Consequently, our rigorous integrator plays a crucial role in providing computer-assisted proofs for the global existence of solutions originating from initial data, particularly those distant from the neighborhood of equilibria.

Spot pattern equilibrium. Let us consider another solution converging to the spot pattern equilibrium solution in Fig. 2b. In this case, the values of C and \(\lambda \) are determined to be \(C=6.7813\) and \(\lambda =0.3505\), respectively. The largest positive \(\rho \) is found to be \(\rho = 0.0046\). Consequently, these values lead to a radius of the target neighborhood around the equilibrium, calculated as \(\rho /C = 6.746\cdot 10^{-4}\). The initial data, denoted by \(u_0(x)\), is plotted in Fig. 4a, which is slightly close to the target equilibrium solution. Nevertheless, the integrator helps us to propagate the rigorous inclusion of the solution into the trapping region. In other words, without the integrator, the global existence of the solution cannot be proved. The integration is carried out until \(t_K=0.375\) (with \(K=12\)), setting \(\tau _i=3.125\cdot 10^{-2}\) (\(i=1,2,\ldots ,K\)).

The error bound \(\varvec{\varrho }_0\) in Theorem 4.3 is obtained by

The \(\varvec{\varrho }_K\) bound is \(\varvec{\varrho }_K=9.6731\cdot 10^{-7}\). Then the value of \(\varvec{\varrho }_K + \Vert \bar{a}^{J_K}(t_K) - \tilde{a}\Vert _{\omega }\) is bounded by \(6.669\cdot 10^{-4}\). This value is less than \(\rho /C\). From Theorem 5.1, there exists a solution globally in time, which converges to the spot pattern equilibrium solution.

The initial data \(u_0(x)\) and the difference between the stable spot pattern equilibrium solution \(\tilde{u}_2(x)\) and \(u_0(x)\): the initial data is close to the equilibrium solution but it is still out of the trapping region. Our integrator is necessary to bring the rigorous inclusion of the solution into the trapping region

Full size image

A notable aspect of this result is that the global existence of the solution has only been proved for initial values close to the equilibrium solution. The reason is the exceptionally small trapping region required for proving global existence, which is related to the stability of this equilibrium solution. In other words, this equilibrium solution has a slow stable manifold, and on such a slow manifold the rigorous integrator is difficult to succeed for a long time period. To successfully demonstrate computer-assisted proofs for global existence from initial values further from the equilibrium, it is necessary to progressively adjust the step size \(\tau _i\) of each time step. Or the accuracy of the approximate solution needs to be improved, which leads to smaller bounds \(\delta ^{J_i}\) and \(\varepsilon ^{J_i}\) in Theorem 4.3. These improvements should be positioned as future works.

In this section, as the second application of the provided integrator, we consider the 2D Ohta–Kawasaki equation

which is a nonlinear PDE used in the study of microphase separation in diblock copolymers, and which is pivotal in material science for predicting and understanding the self-organizing structures in soft condensed matter systems. Developed by Ohta and Kawasaki [57], it models the complex pattern formation in polymer blends, driven by the balance between entropic effects and chemical incompatibility. Note that the Ohta–Kawasaki equation (6.1) has been studied recently with the tools of rigorous numerics, including constructive computer-assisted proofs of existence of steady states [14, 56, 58, 59] and connecting orbits [60].

In this paper, we set \(\epsilon =0.4\), \(\sigma =1\), \(m=0\) (the average of solution, i.e., \(\frac{1}{|\Omega |}\int _\Omega u(x)dx\)) and consider the equation (6.1) under the homogeneous Neumann boundary condition on the 2D rectangle domain \([0,\pi /L_1] \times [0,\pi /L_2]\) with \((L_1,L_2)=(1,1.1)\). We note that the average of solution denoted by m is conserved for any time. In other words, the zero mode of the Fourier coefficient is fixed as \(m=0\) in this example. Using the general notation (1.1), \(\lambda _0=-\sigma \), \(\lambda _1=-1\), \(\lambda _2=-\epsilon ^2\) and \({\Delta ^p N(u)} = \Delta (u^3)\) (that is \(p=1\) and \(N(u)=u^3\)). The corresponding ODEs (2.4) for the time-dependent Fourier coefficients is given by

Setting \(\mu _{{\varvec{k}}}=-\sigma +({\varvec{k}}{\varvec{L}})^2-\epsilon ^2({\varvec{k}}{\varvec{L}})^4\), \(q=2p=2\), and \(\mathcal {N}_{\varvec{k}}(a)=(a^3)_{\varvec{k}}\) in the framework, we integrate (2.4) via the multi-step scheme in Sect. 4 until some \(t_K \ge 0\).

Unlike the earlier Swift–Hohenberg cases, achieving rigorous integration of the Ohta–Kawasaki equation for long time \(t_K\) is challenging. This difficulty primarily arises from the wrapping effect encountered in rigorous forward-time integration, which stems from the overestimation of several bounds presented in Sect. 3. Therefore, for successful integration up to a large \(t_K\), it is essential to carefully choose the time steps \(J_i\) as defined in Sect. 4. In other words, a meticulous partitioning of the time interval \([0,t_K]\) allows for long-time integration. Thus, we employed an “adaptive” procedure based on the error bound at the end point of each time step, as introduced in Sect. 4.2. The step sizes \(\tau _i\) of time steps are automatically adjusted to ensure that the ratio of the error bound at \(t=t_{i-1}\) to the one at \(t=t_i\) (\(i=1,2,\dots ,K\)) does not exceed a predetermined threshold, set in this case to 1.2.

Note that while this “adaptive” procedure is similar to the time-stepping method, the resulting validated error bound is obtained all at once over the entire time interval \([0,t_K]\). This approach corresponds to the multi-step method introduced in Sect. 4, rather than the standard step-by-step procedure.

The Fréchet derivative of \(\mathcal {N}\) at \(\bar{a}\) is the same as that of Swift–Hohenberg, given by \(D\mathcal {N}(\bar{a})\phi =3\left( \bar{a}^2*\phi \right) \) for \(\phi \in \ell _{\omega }^1\). Accordingly, the function g satisfying (2.12) is also the same. As a result, \({\mathcal {E}_{m,\infty }^{J}}\) and \({\mathcal {E}_{\infty ,m}^{J}}\) bounds and \(L_{\bar{a}}(\varrho )=3\varrho \left( 2\Vert \bar{a}\Vert _X+\varrho \right) \) in Theorem 3.15 are equal to those previously determined. The main difference from the case of the Swift–Hohenberg equation in the previous section is the derivative at the nonlinear term. We set \(\nu _{{F}}=1.001\), \(\gamma =0.5\), and \(\xi =0.8\) in Sect. 3.1. The other computational parameters are determined by \({\varvec{N}}=(12,12)\) and \({\varvec{m}}=(3,3)\).

6.1 Computer-Assisted Proofs

We consider two initial data, whose time evolutionary solution profiles go to different stationary states as shown in Fig. 5. Such two initial datum are given by the form

with

Two stationary states of the 2D Ohta–Kawasaki equation (6.1): From the different initial data \(\varphi ^i\) (\(i=1,2\)) in (6.2), the solution profile changes to each stationary state. Such an asymptotic dynamics of (6.1) is determined by the parameters \(\epsilon \) and m (see, e.g., [58, 59] for more details)

Full size image

Computational results of rigorous integration for the Ohta–Kawasaki equation (6.1): Results of the adjusted step size are shown in a, b. In both cases, the initial step size was set to \(\tau _1=0.0625\). After about 100 time steps, smaller step sizes were chosen to control the wrapping effect of the rigorous enclosure. This wrapping effect arises from an overestimation of the uniform bounds for the evolution operator presented in Sect. 3.1. By using Theorem 4.3, the rigorous error bound displayed in c was obtained all at once over the entire interval [0, t]. Note that these results were not obtained using a step-by-step procedure

Full size image

Stripe pattern state. Set the initial data as \(\varphi ^1\) in (6.2). We chose the step size \(\tau _i\) at each time step \(J_i\), as illustrated in Fig. 6a. Resulting computer-assisted proof assures that the solution of the IVP of (6.1) exists at least to \(t_K=8.4666\) (\(K=194\)). The rigorous error bound \(\varvec{\varrho }_0\) via Theorem 4.3 satisfies

We also have \(\varvec{\varrho }_K=1.2838\).

The results shown in Fig. 6a and c indicate that after about 100 steps, while using a smaller step size helps to control the increasing error bound to some extent, such error bound becomes larger and more difficult to handle. Eventually, this error escalation makes the hypothesis of Theorem 4.3 unverifiable. However, our result marks the first implementation of rigorous integration method for time-dependent PDEs in higher spatial dimensions, including nonlinear terms with derivatives. We believe that this method constitutes a significant advance in the methodology of rigorous forward integration for PDEs and demonstrates satisfactory performance as a rigorous integrator.

Spot pattern state. Setting the initial data as \(\varphi ^2\) in (6.2), the choice of the step size is displayed in Fig. 6b. Computer-assisted proofs verifying the existence of a solution for the IVP of (6.1) succeeded up to \(t_K=9.0577\) (with \(K=216\)). The rigorous error bound \(\varvec{\varrho }_0\), as in Theorem 4.3, is confirmed to satisfy the following inequality:

Notably, \(\varvec{\varrho }_K\) at the final step is obtained as 2.0647.

Figure 6c also shows the escalation of the error bounds, due to the wrapping effect. In this case, the error bound is slightly smaller than that of the previous (stripe pattern) case, indicating the success of the integrator over more multiple time steps. To achieve successful long-time rigorous integration, it is necessary either to make the step size more reasonable, as discussed at the end of Sect. 5.3.2, or to reduce the defects at each time step.

We conclude this section by mentioning that while it would be natural to consider a computer-assisted proof for global existence in the Ohta–Kawasaki equation, as was done in Sect. 5 for Swift–Hohenberg, this would require introducing a nontrivial construction for the trapping region (as the one used in the paper [30, 31]), and therefore we do not follow this path here.

In this paper, we have developed a general constructive method to compute solutions for IVPs of semilinear parabolic PDEs. Our approach combines the principles of semigroup theory with the practicality of computer-assisted proofs. A feature of our approach over the ordinary mathematical-analytic approach is the Fourier–Chebyshev series expansion, which derives a numerical candidate for the solution, thereby enhancing computational accessibility and feasibility. We have constructed a two-component zero-finding problem as a direct derivation from the original PDEs. To bring in the fixed-point argument, we have introduced a Newton-like operator at the numerically approximated solution. Central to our approach is the use of the evolution operator, which is the solution map for linearized equations at the numerical approximation. By inverting the linearized operator “by hand”, the Newton-like operator is effectively derived in a more direct way. The existence and local uniqueness of the fixed point of the Newton-like operator, together with its rigorously determined bounds, provide the rigorous inclusion of the solution to the IVPs. Computer-assisted proof is done by numerically verifying the hypothesis of the contraction mapping via interval arithmetic.

We have further generalized our approach to a multi-step scheme that extends the existence time of solutions to IVPs. This involves considering a coupled system of the zero-finding problem over multiple time steps, where the Newton-like operator at numerically approximated solution is derived by constructing the inverse of the linearized operator. By numerically validating the hypothesis that the Newton-like operator becomes the contraction map, we obtain a rigorous inclusion of the solution over these multiple time steps. This multi-step approach should enhance the applicability of our approach not only for IVPs but also for boundary value problems in time, which is the subject of future studies.

Constructing the trapping region based on the mechanism of convergence towards asymptotically stable equilibria, we have also presented proofs of global existence of solutions to three-dimensional PDEs (Swift–Hohenberg) converging to a nontrivial equilibrium. This is the first instance of CAPs being applied to global existence of solutions to 3D PDEs, offering a new perspective in the field. Moreover, our integrator has been applied to the 2D Ohta–Kawasaki equation, dealing with derivatives in the nonlinear terms. The multi-step scheme, together with the smoothing property of the evolution operator, has allowed for some effective control of the wrapping effect, and long-time rigorous forward integration has been successfully achieved. Our study not only contributes to the theoretical aspects of these equations but also provides practical tools and methods for computer-assisted proofs for their resolutions, opening new avenues for study and application in the field of computer-assisted proofs.

We conclude this paper by highlighting several potential extensions of our rigorous integrator. First, our approach provides a cost-effective way to implement computer-assisted proofs for multi-steps. This aspect is particularly noteworthy because it allows the linearized problems to be solved independently and a priori at each time step, making the process naturally parallelizable and keeping the computational cost additive rather than multiplicative. In other words, the manual inversion of the linearized operator, a critical part of our approach, significantly reduces the computational complexity by avoiding the inversion of large finite-dimensional matrices. We believe this feature could be beneficial in studying boundary value problems over time (such as periodic orbits and connecting orbits), which is a topic for exploration in future research.

Second, exploring beyond the spatial domain assumption (rectangular domain) is an interesting direction of this study. For example, it is worth considering the potential of extending computer-assisted proofs to the time-dependent PDEs on more complicated domains, such as disks/spheres, etc. This would require suitable basis functions to approximate solutions satisfying boundary conditions. For effective computer-assisted proofs, the property of geometric decay for the coefficients of the series expansion (such as Fourier/Chebyshev) is desirable, and the property of Banach algebra for discrete convolutions is useful for handling function products.

Finally, a challenging application of our approach lies in addressing the incompressible 3D Navier–Stokes equation on the domain \(\mathbb {T}^3\). Here, we could use the Fourier setting of the vector fields as described in [15] to construct the associated infinite-dimensional ODEs. One could also obtain a trapping region of the solution asymptotic to the stable equilibrium (zero) solution. Achieving computer-assisted proofs for such complex equations is one of the ultimate goals in this line of research.