Rigorous computation of non-uniform patterns for the 2-dimensional Gray-Scott reaction-diffusion equation

In this paper a method to rigorously compute several non trivial solutions of the Gray-Scott reaction-diffusion system defined on a 2-dimensional bounded domain is presented. It is proved existence, within rigorous bounds, of non uniform patterns significantly far from being a perturbation of the homogenous states. As a result, a non local diagram of families that bifurcate from the homogenous states is depicted, also showing coexistence of multiple solutions at the same parameter values. Combining analytical estimates and rigorous computations, the solutions are sought as fixed points of a operator in a suitable Banach space. To address the curse of dimensionality, a variation of the existing technique is presented, necessary to enable successful computations in reasonable time.

a pair of reactions involving cubic autocatalysis [3]. It consists in the following reaction-diffusion system where U and V are the concentrations of chemical reactants U and V, d 1 , d 2 are the diffusion coefficients, F is the feed rate and (F + κ) the removal rate of V. Experiments and numerical investigation reported by Pearson [4] and in [5,6] reveal a rich and complex structure in the solutions for the Gray-Scott equation, including self-replicating patterns, oscillating patterns, spots annihilation, irregular patterns, spatio-temporal chaos. Since then, the behaviour and the dynamics of the patterns in the Gray-Scott system has been extensively studied through experimental observations, numerical techniques and analytic approach.
In two dimensional domains, particular solutions like single spot solutions, ring-like solitons and stripe patterns have been found and analysed in [9][10][11]15]. Also, an hybrid analytic-numerical approach is adopted in [12] to detect multi-spot quasi equilibrium patterns.
In this paper analytical theory and rigorous numerics are combined in a computer assisted technique, that allows to successfully treat non-local problems while preserving the mathematical rigorousness. We prove existence, together with precise bounds, of several non homogenous stationary patterns, arranged on continuous branches, for the 2-dimensional Gray-Scott system defined on a bounded domain Ω and with no-flux boundary condition. Because of the extreme difficulty in studying non linear differential equations by purely analytical approach, any new methodology that demonstrates qualitative and quantitative properties of solutions for systems of nonlinear PDEs is a valuable achievement by its own. Moreover, the results of this work can contribute to the interesting discussion concerning the genesis of localised patterns and the underlying mechanisms that drive the transition from one dynamics to another. For instance, they can be useful to understand the self-replicating pattern dynamics, that is, the itinerary process that moves a localised trigger to a stable stationary or oscillating pattern by multiple splitting of the pulses. Attempts to answer these questions are presented in [13,14]. Looking at the global bifurcation diagram of stationary solutions, in [14] it has been suggested that a hierarchy structure of saddle-node points is responsible for the occurrence of self-replicating pattern dynamics. In [13] the relationship between different bifurcation branches of steady states is considered as a criteria for the onset of spatio-temporal chaos. The spatio-temporal chaotic dynamics emerges by unfolding an itinerant heteroclinic cycle. Also, self-destruction and spots annihilation can be explained as an orbit starting close to an unstable pattern whose unstable manifold is connected to the homogenous state. From these considerations, it is evident that the bifurcation diagram of steady states together with the stability analysis plays a key role in the comprehension of the mechanism behind the complex pattern dynamics.
In order to introduce the results, let us present more in details the system we are dealing with. Following [7], the Gray-Scott system is equivalent to the following For any choice of λ and γ, problem (2) admits a stable uniform stationary solution (u, v) = (1, 0) called p 1 and, for λ ≥ 4, two branches of uniform stationary solutions denoted by p 2 (λ) and p 3 (λ). Local bifurcation analysis predicts the values of γ where families of non-constant stationary solutions bifurcate from the homogeneous states. We investigate various bifurcation branches of non-constant solutions emanating from the homogenous states and the goal is to prove existence of several non constant patterns, providing explicitly the numerical data together with the enclosing bounds.   Looking at Figure 1, we prove the following: -Any bullet in Figure 1(a) represents a non-uniform stationary solution for system (2) belonging to branches that bifurcate from the homogenous state p 2 (λ). -Any bullet in Figure 1(b) represents a non-uniform stationary solution for system (2) belonging to branches that bifurcate from the homogenous state p 3 (λ).
Clearly, by changing the parameter λ and the size of the rectangular domain Ω, many non-constant solutions on different branches can be computed and rigorously validated. Also, most of the families depicted in Figure 1 can be easily extended. The computational methods used to prove Theorem 1 has roots in the radiipolynomial technique, introduced in [17] with the aim at solving nonlinear differential problems. In this paper a variant of the computational method is presented. It consists in the introduction of a sharpening parameter that allows to obtain more precise estimates of the nonlinear terms without enlarging the dimension of the finite dimensional approximation. Indeed, using the computational technique in the form it has been developed so far, the computational cost and the time required for the enclosure of a single solution stated in Theorem 1 would be so large as to be practically unfeasible.
The paper is organised as follows. In section 2 we discuss the stationary solutions for system (2) and, for fix λ, we find the bifurcation values γ k for branches of non uniform solutions. Section 3 concerns the radii polynomial technique. First we present the computational method in the general form together with the notion of radii polynomials. Then, in section 3.2 some fundamental formulas for the analytic estimates of the convolution sums are recalled. In section 4 we adapt the computational method to the Gray-Scott system. We rephrase the differential system into an equivalent fix point problem and we introduce the Banach space within which we look for the solutions. Then we write the radii polynomials as they would be without the modifications. In section 4.3 we explain in more details the reason why the radii polynomials technique, so far presented, fails in proving the existence of the solutions we are interested in. In section 4.4 the sharpening parameter is introduced and the new scheme to bound the nonlinear terms is discussed. Finally, in section 5 we are concerned with the computational results.

Some analytical preliminaries: bifurcation values
Stationary patterns of (2) consist in solutions of the elliptic problem Throughout this paper we denote by h k the eigenvalues of and by φ k the associated eigenfunctions. We restrict our analysis to the case of rectangular domain, Ω = [0, . Denoting x = (z, y) ∈ Ω, the eigenvalues and eigenfunctions are explicitly given as for any k = (k (1) , k (2) ) ∈ N 2 . For any value of λ and γ, system (3) admits the uniform state (u, v) = (1, 0), denoted by p 1 . At λ = 4 a bifurcation occurs and two branches of uniform solutions are created. The set of uniform solutions is the following: Remark 1 Note that for any uniform solution (u * , v * ), except p 1 , it holds u * v * = 1.
For each λ > 4 we now detect the values of γ where families of not uniform solutions bifurcate from the homogenous state. Denote by (u * , v * ) any of the constant stationary solutions p 2 (λ), p 3 (λ). For any k ∈ N 2 \ {(0, 0)}, define a perturbation of the steady state in the direction of the eigenfunction φ k . Following the standard bifurcation theory, a bifurcation occurs at value γ k if for any ε small enough it exists a solution for (3) of the form (u k (ε), v k (ε); γ(ε)), being γ(ε) = γ k + εγ + o(ε). Inserting (u k (ε), v k (ε)) and γ(ε) into the elliptic problem (3) and using the fact that (u * , v * ) is a solution for any γ, we obtain the following Neglecting the infinitesimal terms, for Remark 1, it remains to solve the system 0 whose solution is given by and We summarise the result in the next Lemma.

The Computational method
The method here presented is one of the so called verification methods [18]. Suppose we are looking for solutions of f (x) = 0 where f is defined in some Banach space (X, · X ). The goal of a verification methods is to prove the existence of a genuine solution x * of f (x) = 0 in the neighbourhood of a numerical approximate solutionx, providing explicit and rigorous bounds for x * −x X . In this paper we adopt the radii polynomial technique that relies on the Banach fixed point theorem. In the last ten years the radii polynomial technique has been exploited to study a variety of problems in areas ranging from, but not restricted to, ODEs, PDEs, delay differential equations and dynamical systems. See for instance [16, 17, 25-27, 29, 30] and reference therein. In [8] the radii polynomial technique has been applied to the 1-dimensional Gray-Scott system to establish the existence of symmetric homoclinic orbits to the steady state p 1 = (1, 0). More recently, in [25] bifurcation branches of solutions for a system of one dimensional reaction diffusion equations have been rigorously validated using a similar technique.
Of course there are many other methods for computer assisted and rigorous computational study of differential equations, based for instance on topological arguments like the Conley index theory, covering relations with cone conditions, the Shauder fixed point theorem. We refer to [20][21][22][23][24] and the references therein for a more detailed discussion of topological tools applied to rigorous computations.
Before presenting the technique, let us fix some notation. Boldface is used to denote multi-indeces as k = (k (1) , . . . , k (d) ) ∈ Z d . When applied to multi-indeces, the absolute value | · | acts component-wise, |k| := (|k (1) |, . . . , |k (d) |). Inequality in between multi-indeces is intended component-wise, that is for k, n ∈ Z d , k < n means k (j) < n (j) for all j = 1, . . . , d. The same for ≥, ≤, >. In the following we will introduce the multi-indeces m, M and M all in N d , respectively called finite dimensional parameter, FFT dimension parameter and computational parameter, satisfying m ≤ M ≤ M . Associated to them, we define the set of indices F m , F M , F M as F m := {k ∈ Z d : |k| < m} and similarly for the others. Note that the difference F M \ F m = {k ∈ Z d | ∃j : m (j) ≤ |k (j) | < M (j) } contains all the indices having at least one ( at not necessarily all) component

Overview of the radii polynomial method
We describe the technique in the context of solving the PDE defined on a bounded domain Ω ⊂ R d , subjected to prescribed boundary conditions. Here L : D(L) ⊂ H → H is a linear operator densely defined on a Hilbert space H and p is the degree of nonlinearity, c p ∈ R. The first step of the method consists in writing an equivalent problem defined on a suitable Banach space, so that the zeros of f (x) correspond univocally to solutions of (8). Assume that the Hilbert space H admits an orthogonal basis formed by the eigenfunctions of L(u), say {φ k } k with eigenvalues {h k } k and assume that the domain of L is given by By expanding (8) on the basis {φ k } k , it follows that the solutions of (8) correspond to the zeros of f (x), x = {x k } k , given by If the eigenfunctions basis {φ k } k is of the form of a Fourier basis, as it will be in our case, the inner product < u p , φ k > H is equal to convolution products of the Fourier coefficients {x k } k , that is For a vector s = (s (1) , . . . , s (d) ) ∈ N d , let the weight ω s k be defined as Define the s-norm of a sequence x = {x k } k as x s := sup k∈Z d {ω s k |x k |} and introduce the Banach space The space X s is the space of sequences algebraically decaying, with decay rate s, and it is the space where we look for solutions of f (x) = 0. This choice is motivated by the fact that any solution u of the PDE (8) is expected to be smooth and the sequence of Fourier coefficients of a smooth function decays exponentially fast to zero. Therefore, it is reasonable to look for the solution {x k } k in the space of sequences that decay at least as fast as { 1 The proof of the existence of solutions for f (x) = 0 relays on a contraction mapping argument applied to an operator T , defined below, whose fixed points are in one to one correspondence with the zeros of f . The operator T depends on the Frechet derivative of f evaluated at a numerical approximate solutionsx. For Letx F m be a numerical solution for f (m) (x) = 0 and denote byx the embedding of By construction, J −1 is an approximate inverse of Df (x). Note that for the tail part, that is for k / ∈ F m , the operator J −1 does not rely on numerical computations, rather it is defined taking into account the linear contribution of f k only. The fixed point operator is defined as T : The goal of the technique is to detect a ball B(x, r) ⊂ X s within which T is a contraction. That is done through the notion of the radii polynomials, a finite set of inequalities in the variable r, the verification of which implies that T is a contraction on some ball B(x, r). They are defined as combination of the following bounds.
Let {Y k } k and {Z k (r)} k be two sequences such that, for any k ∈ Z d , Furthermore, for a choice of a computational parameter Proof The result follows as an application of the Banach fixed point theorem on the operator T . See for instance [17,19].
Before discussing some analytical estimates used in the construction of the bounds Y k and Z k let us close this section with a remark.

Remark 2
-The method requires that the eigenvalues and eigenfunctions of the linear operator L(u) are explicitly known. That is the reason why we consider only rectangular domain Ω with Neumann boundary conditions. However, only few and trivial modifications are required if periodic boundary conditions are posed. -The fact that the polynomial nonlinearity is transformed into convolution sums is a peculiarity of the exponential function basis. When a different complete function basis is chosen the polynomial u p may be translated in something similar to the convolution sum, as in the case of the Chebyshev function basis ( that are nothing more that Fourier series in disguise [28]), or in something completely different as in the case of Hermite functions.

Analytical estimates of the convolution sums
One of the fundamental steps of the technique is the definition of the bounds Y and Z satisfying the inequalities (10) and (11). In the construction of these bounds it is required the computation of several convolution sums of the form and sharp estimates for convolutions series A deep and systematic analysis of the latter series is done in [31], where analytical estimates have been proved for any convolutions terms. The former one (12) is a finite sum and the Fast Fourier Transform (FFT) can be used to faster the computation. Nevertheless, the computational time dramatically increases as far as M , and in particular the dimension d, increase. Moreover, aiming at a rigorous enclosure, these operations must be performed in interval arithmetics regime, hence requiring an even larger computational resources. For instance, if d = 3 and M = 200, the rigorous computation of (12) involves computing the FFT of sequences of size 10 9 , that is a serious computational task.
These complications and in particular the curse of dimensionality have been addressed in [32]. In that work the authors improved the analytical estimates of [31] and introduced the FFT parameter M with the aim of reducing the computational cost required to enclose the solution of the finite sum (12). For convenience, we report the results that will be used in the sequel of this paper. The detailed construction of the various parameters is given in the Appendix.
The first estimate concerns the upper bound of the convolution series (13) when the indices are free to move in the full set Z d .
be defined as in (33). Then, for any k ∈ Z d , Proof See Lemma 2.1 in [31]. Now, we consider the case when some of the indices are restricted to take value outside a certain set. Suppose M = (M (1) , . . . , M (d) ) has be defined so that M (j) ≥ 6 and M ≥ M . The notation {k 1 , . . . , k +1 } ⊂ F M means that at least one of the k j is not in F M . Such a set corresponds to the complement of be defined as in (35). For any k ∈ F M and = 0, . . . , n − 1 Proof See Lemma 2.2 in [31].
ii) The boundary conditions in (3) are automatically satisfied by any u, v given in (14).
Again from (5), we remind that Since to solve system (3) it corresponds to solve the infinite dimensional algebraic system In each g k we replace the second line with the sum of the first and the second and we obtain the equivalent system From the second equation it follows thus, plugging v k into the first equation, we end up with a system of equations in the variables {u k } k only. However, to easy the notation, we continue to use the variable v k , having in mind that any v k is function of u k . Let us define and, given a sequence x = {x k } k , we adopt the notation qx for the sequence qx = {q k x k } k . In particular v = qu. Also, introducê Remark 3 The values of the parameters λ and γ of interest in this work are such that λγ < 1. In this situation it holds |q k | ≤q for any k ∈ Z 2 .
The system we need to solve is in the unknown u = {u k } k∈Z 2 . We look for solutions in the Banach space where the s-norm is the same as in section 3.1. Introduce to denote the coefficients of the linear terms in f k .

Construction of the polynomial bounds
To avoid unnecessary verbosity, most of the technicalities are skipped, rather we report the explicit definition of the bounds Y and Z. For more details we refer to [32]. Choose the finite dimensional parameter m, the computational parameter M so that M ≥ 3(m − 1) − 1 and the FFT dimension parameter M so that m ≤ M < M .
Suppose that a numerical solutionū F m is computed so that f (m) (ū F m ) ≈ 0. Denote byū the embedding ofū F m into X, i.e.ū = (ū F m , 0 I m ), and definev = qū. Compute the finite dimensional matrix J −1 m and the infinite dimensional operator J −1 as in (9). For this, it is convenient to report the jacobian matrix Df (u). Having in mind that the v k is function of u k , it holds The construction of Y k readily follows from the definition (10): The last line is due to the choice of M , for which (ū * v * v) k = 0 for any k / ∈ F M . The polynomial Z k (r) is defined as and it is rigorously computed, since only a finite number of computations are involved.
Combining rigorous computation and the analytical estimates, all the other bounds can be defined as: The tail coefficient (11) is given bỹ M r +Z

A priori analysis of the radii polynomials
According to Lemma 2, the existence of a stationary solution for the Gray-Scott equation is assured if the radii polynomials are all negative for some positive value of r. A careful analysis of the coefficients of the polynomials clarifies under which conditions the computation will be successful and feasible.
For k ∈ F M the polynomial p k (r) is of the form p k (r) = a 0 + a 1 r + a 2 r 2 + a 3 r 3 . The constant term a 0 is always positive and it is related to the accuracy of the numerical solution. It can be lowered if a better approximate solution is produced, for instance by increasing the finite dimensional projection m.
Also a 2 and a 3 are positive and they are related to the rate of expansion of the operator T due to superlinear terms in the function f (x). A necessary condition for the polynomial p k (r) to be negative somewhere is that a 1 is negative. Looking to the case k ∈ F M \ F m , the coefficient a 1 is of the form a 1 = 1 Since Z (1) k decays like Cω −s k (because of the properties of the convolutions), it is enough that |µ k | is increasing to ensure, for m large enough, that a 1 is negative. The coefficients µ k s are related to the eigenvalues of the linear operator in the PDE. In particular, high order linear operators yield larger µ k and enhance the contractivity of T , that translates into a smaller value for a 1 and a possible choice of smaller m. Besides that preliminary theoretical argument, a crucial issue is the feasibility of the computation. Indeed, it is necessary to choose the various parameters m, M and M relatively small, so to abet practical computation. In particular, since the convolution product is the most costly computation, and it appears several times in (20), it is fundamental to keep M as small as possible.
Let us now inspect how the growth rate of µ k influences the choice of M in the context of the PDE under investigation. Looking at the polynomial p k (r) with k ∈ F M \ F M , we realize that The parameter M must be chosen large enough so that µ k remains reasonable larger than α (3) k for any k / ∈ F M . Again, as a general statement, a larger growth for µ k allows a smaller choice of M . For two-dimensional PDE problems the coefficient α (3) k rapidly grows to large values, hence the parameter M needs to be selected quite large. Our situation is even worst due to the presence ofq = γ −1 that multiply the effect of α (3) k . To have a quantitative idea, consider the family of solutions labeled with (0, 3) in Fig. 1(a). It bifurcates at the value γ ≈ 0.0087. Taking into account all the constant factors and the s-norm ofū and v ( s = (2, 2)) it is necessary to choose M = (M 1 , M 2 ) with both M 1 , M 2 larger than 250. To perform cubic convolutions of sequences with around (500) 2 terms is extremely computationally expensive. We underline that it is the combination of the high dimensionality of the problem and the presence of a small constant γ in front of the linear term that makes the computation of the radii polynomials not affordable. Indeed the method was successfully applied in [31] to prove equilibria for the 2D and 3D Cahn-Hilliard equation with constant factor in the range γ ≈ [0.05, 0.1] and γ ≈ [0. 3,1] respectively. Rather, in [32] it is proven existence of solutions for the Allen-Cahn equation with much smaller factor γ ≈ 0.001 but only for the one dimensional case.
In conclusion, we are not claiming that the technique fails in proving the solutions we are interested in, but that, just for one tentative enclosure, one needs to perform a number of computations so large that can not be addressed by a standard computer in reasonable time. Therefore, to the goal of proving Theorem 1, the method so far described is not feasible. To overcome these hurdles, in the next section we introduce some modifications. k . Although it was not said explicitly, the role of these terms is to provide a bound for the series

New definition of the bound Z
and k 1 +k 2 +k 3 =k |k 1 |,|k 2 |∈Fm ,k 3 ∈Z 2 The first is obtained applying Lemma 4 and the uniform estimates |ū k | ≤ ū s ω −s k , |v k | ≤ v s ω −s k and |q k 3 | ≤q. Similarly, the second bound has been defined applying Lemma 3.
Rather than uniformly estimate each coefficient |ū k |, |v k | as done above, we can separate in the series (22) the significative contributions from the smaller one. Indeed, it could be the case that few of theū k 's andv k 's are large, while the others are extremely small. Such a behaviour is expectable, for instance, when the aimed solution belongs to a branch bifurcating from the constant state. Typically, in this situation, only one mode is dominating.
For a given , let us define We refer to as the sharpening parameter and we assume that < min{ ū s , v s } so that both F m (ū, ), F m (v, ) are not empty. They are clearly finite, being subsets of F m .

The series (22) is decomposed into
(25) The effect of the sharpening parameter is twofold. First, it allows to isolate the main contribution of the series, given by the first sum on the right hand side of (25) that will be rigorously computed. All the others sums are uniformly bounded by means of the bound ε (3) k , see Lemma 4. We note that anytime the index |k 1 | ranges outside F m (ū, ) (resp. |k 2 | ranges outside F m (v, )) it holds |ū k 1 | ≤ ω −s The last estimate is a much sharper than |ū k | ≤ ū s ω −s k previously adopted. Explicitly, we obtain ω s k .
The advantage of this approach is to enfeeble the effect of γ −1 hidden inq. On the other side we have to calculate one more sum at any stage, requiring further computational time. However, if the mass ofū is concentrated in few Fourier coefficients, such a computation is pretty fast. Inserting these newly defined bounds, we redefine Z  (20) 2q ū s v s ε

Concerning the bound Z
(2) k , we aim at replacing the terms 6 v sq ε and similarly the terms 2 ū sq 2 ε The former are defined so that Arguing as before, i.e. collecting the largest contributions in the series, and using the bound q k ≤q and Lemma 4, we write Similarly, the same procedure provides bounds for the second of (26) and for the We omit the details. Collecting all the terms, we reformulate bounds as reported before.
It remains to define the tail boundZ M . By definition, we have to find a constant Z For k ∈ F M , Z (1) k is now given by (28) According to [31] the last two terms in (28) can be easily bounded in terms ofα given in (34). The first one is less than 2qS k where Having in mind that k = (k (1) , k (2) ), k 1 = (k 2 ), M = (M (1) , M (2) ), s = (s (1) , s (2) ), let us introduce the following quantities.
The next Lemma provides an uniform bound for ω s k S k for any k ∈ F M .

Results
In this section we report the results obtained by applying the radii polynomial technique, in the form discussed in the previous section, for the enclosure of non constant solutions for system (3). All the rigorous computations and the visualisation have been performed in Matlab with the interval arithmetic package Intlab [33]. We remind that the domain Ω is given as the rectangular In any computation the decay rate s is fixed equal to s = (2, 2). We do not investigate the performance of the enclosure as s varies. At this regard, an interesting analysis has been performed in [25].
For convenience, let us recall the homogenous states First we discuss how to compute the numerical solutions, then we report some details about the rigorous enclosure.

Computing the numerical solution
For a choice of the domain sizes L 1 , L 2 and of the parameter λ we numerically compute several patterns lying on the same family bifurcating from one of the homogenous state p 2 (λ), p 3 (λ). Each solutions is intended in Fourier space, that is, we numerically find a zero of a finite m-dimensional projection of system (16) in the unknowns X = {γ, {u k } k∈F m , {v k } k∈F m }.
According to the analysis performed in section 2, for a choice of k we set that value γ = γ k and we perturb the homogenous state in the direction of the eigenfunction φ k . This perturbed state, together with γ k , is now considered as initial data for a Newton iteration scheme. Once a first non homogenous state is computed, we adopt a numerical continuation scheme to compute other solutions on the same branch. We slightly move γ in the direction of the branch, we perturb the previous computed state along the tangent direction and we look for a new solution in the unknowns {u k } k∈F m . There are different parameters one can play with that influence the outcome of the process. For example the length of the increment of γ, the size m of the finite dimensional projection and also the accuracy requested for the numerical solution. In our computations we assume that the numerical solutionx is achieved when the sup-norm of f (m) (x) is less than 10 −10 . The parameter m = (m, m) ranges in between 12 ≤ m ≤ 23 and different ∆γ-step are considered. Here we discuss the results stated in Theorem 1 and depicted in Figure 1. A magnification of Figure 1(a) is given in Figure 2. Any bullet in the figures represents a validated non-constant pattern. The graphs depict γ against the squared L 2 norm of the v − v 0 , where v 0 denotes the v component of the homogenous state. The different families are labelled according to the bifurcation point, that is the (k (1) , k (2) )-branch is the one that bifurcates from the homogenous state at value γ (k (1) ,k (2) ) .   Table 1. Table 1 reports the parameters m, M , M and the sharpening coefficient used in the computation of the solutions marked in Figure 2 together with the enclosure radius. Since s = (2, 2), the value of r provides also a bound for the C 2 distance of the exact solution from the numerical approximation.

Rigorous enclosure results
From the data in the Table 1, we realise that larger finite dimensional parameter and larger computational parameter are required when γ decreases and the size of k increases. The role of γ is already discussed in section 4.3. Although the introduction of the sharpening coefficient enables to weaken the aftereffect of γ, a dependence between the parameters is still present and a smaller γ requires a larger M and M . On the other side, k denotes the Fourier coefficient of u and v that has been initially perturbed at the bifurcation point. Due to the nonlinear effects, such a perturbation is spread on the neighbourhood coefficients. Hence a larger k requires larger m and M to take into account a large number of contributing    Table 1: Values of the parameters used in the computation of the solutions labeled in Figure 2. The last column is the radius of the ball in X s , s = (2, 2) around the numerical solution within which the genuine solution is guaranteed to exists.
coefficients. Also, further the solution is from the bifurcation point, more important is the effect of the nonlinearity. A similar analysis, here skipped, can be performed for the enclosure of the patterns bifurcating from the homogenous state p 3 (λ) and reported in Fig. 1(b). In Figure 3 some of such patterns are plotted.

Example 2. Effect of the sharpening parameter.
We now investigate the role of the newly introduced sharpening parameter .
Remember that the small is the higher is its effect, in the sense that the sets F m (ū, ), F m (v, ) given in (24) are larger. Consider the solution in correspondence to γ = 0.0108 on the (1, 2) bifurcating branch in Figure 1(b). In Figure 4 both the components u and v of the pattern are depicted. For different choices of , in Table  2 we list the computational parameters M and M needed to obtain a successful validation of the numerical solution. The finite dimensional projection parameter is m = 13. We see that for large extremely large values for the computational parameters are requested. The last row in the table concerns a value close to the largest possible choice of , according to the condition ≤ min{ ū s , v s }. Therefore, if the sharpening parameter is not considered at all, that is, if we follow the approach of [32], even larger computational parameters must be selected and a rigorous enclosure of the solution would be practically impossible.   The solution belongs to the family (1,2) reported in Figure 1(b) at γ = 0.0108.

Conclusion
In this paper a new development of the radii polynomial technique is proposed with the aim of rigorously computing non-uniform patterns for the two dimensional Gray-Scott system. It is demonstrated that the new approach succeeds in enclosing several complex patterns, solutions that could not be proved following the previous approach present in the literature. We plan to develop further the technique and to investigate in different directions. For instance, the three dimensional case may be considered. We propose to combine the method with the theory introduced in [25] to compute parameter dependent smooth branches of solutions. Then a systematic analysis of the bifurcation branches of non-uniform patterns may be performed, to prove whether some branches undergo secondary bifurcations or reconnect to the uniform solution. Also, a more challenging project is to rigorously prove existence of heteroclinic connections and homoclinic cycles.