Optimisation of the total population size with respect to the initial condition for semilinear parabolic equations: Two-scale expansions and symmetrisations

In this article, we propose in-depth analysis and characterisation of the optimisers of the following optimisation problem: how to choose the initial condition $u_0$ in order to maximise the spatial integral at a given time of the solution of the semilinear equation $u_t-\Delta u=f(u)$, under $L^\infty$ and $L^1$ constraints on $u_0$? Our contribution in the present paper is to give a characterisation of the behaviour of the optimiser $\overline{u}_0$ when it does not saturate the $L^\infty$ constraints, which is a key step in implementing efficient numerical algorithms. We give such a characterisation under mild regularity assumptions by proving that in that case $\overline{u}_0$ can only take values in the"zone of concavity"of $f$. This is done using two-scale asymptotic expansions. We then show how well-known isoperimetric inequalities yield a full characterisation of maximisers when $f$ is convex. Finally, we provide several numerical simulations in one and two dimensions that illustrate and exemplify the fact that such characterisations significantly improves the computational time. All our theoretical results are in the one-dimensional case and we offer several comments about possible generalisations to other contexts, or obstructions that may prohibit doing so.


Scope of the article
In this article, we propose to establish several results concerning an optimal control problem for a class of semilinear parabolic equations. Some aspects of this problem have been initially addressed by two of the authors in [44]. In the setting under consideration, the control variable to optimise is the initial condition. As we will see throughout the statement of the results, the form (e.g. convex or concave) of the semilinearity plays a crucial role in the analysis and calls for a detailed study of second order optimality conditions, which is our main result, Theorem 1. In the case of convex semilinearities, using rearrangement arguments, we can give a full characterisation of maximisers, see Theorem 2. Using Theorem 1, we can improve an algorithm initially developed in [44], and we display numerical results in Section 4.
Initial motivation of the paper The origin of this paper is the study of an optimal control problem that arises naturally in mathematical biology and that deals with bistable reactiondiffusion equations. Namely, for a semilinear equation, what is the best possible initial condition ("best" being understood as maximising the integral of the solution at a certain time horizon T )? The complicated behaviour of bistable non-linearities, which are neither convex nor concave, makes the analysis of this query very intricate. The two aforementioned results, Theorems 1 and 2, enable us to show how complicated the behaviour of maximisers can be for such non-linearities. Bistable equations are of central importance in mathematical biology [43] and, very broadly speaking, model the evolution of a subgroup of a population. Among their many applications, one may mention chemical reactions [47], neurosciences [19], phase transition [35], linguistic dynamics [48] or the evolution of diseases [43]. The last interpretation is of particular relevance to us, given that this model is used to design optimal strategies in order to control the spread of several mosquito borne diseases such as the dengue; this was the main motivation in [3]. The strategy is to release a certain amount of Wolbachia carrying mosquitoes (Wolbachia is a bacterium that inhibits the transmission of mosquito borne diseases that individuals inherit from their mother) in a population of wild mosquitoes that can potentially transmit the diseases, in order to maximise the proportion of Wolbachia carrying mosquitoes at the final time. In mathematical terms: given a time horizon T , How should we arrange the initial population in order to maximise the population size at T ?
Even without having stated it formally, we can make two observations on this problem: the first one is that , since the variable of the equation is the proportion of a subgroup, we need to enforce pointwise (L ∞ ) constraints. The second one is that we naturally have to add an L 1 constraint for modelling reasons. Both of these constraints can in practice be very complicated to handle.
Optimisation problems in mathematical biology Let us briefly sketch how this problem fits in the literature devoted to such optimisation and control problems for mathematical biology. Optimisation problems for reaction-diffusion equations have by now gathered a lot of attention from the mathematical community. Most of these optimisation problems are set in a stationary setting, that is, assuming that the population has already reached an equilibrium, and the main problems that have been considered often deal with the optimisation of the spatial heterogeneity [18,22,23,29,30,38,39,41,45] (we also refer to the recent surveys [26,40]); most of these works deal with monostable nonlinearities. We also point to the recent [14] for the study of an optimal control problem for parabolic monostable equations. On the other hand, optimisation problems for bistable equations, which are the other paradigmatic class of equations in mathematical biology [43], have received a less complete mathematical treatment, but are now the topic of an intense research activity from the control point of view, see [3,44] and the references therein. Related optimal control problems are not yet fully understood. More generally, less attention has been devoted to optimisation problem with respect to the initial condition for such semilinear evolution equations.

Statement of the problem
We work in Ω = (0; π). We consider a C 2 function f : [0; 1] → IR, and the associated parabolic equation where u 0 is an initial condition satisfying the constraint 0 u 0 1.
Since our initial motivation, as explained in the first paragraph of this introduction, is to maximise the proportion of a subgroup of a population, such an L ∞ constraint is natural. At the mathematical level, it should be noted that we could carry out the same analysis with any L ∞ constraint of the form 0 u 0 κ by a simple change of variable. We define, for any T > 0, the functional The goal is to maximise J T with respect to u 0 . Since we are then again wondering how to maximise the proportion of a subgroup by controlling its distribution at the initial time, it is natural to introduce a L 1 constraint on u 0 . This constraint is encoded by a parameter m ∈ (0; |Ω|) which is henceforth fixed.
These considerations lead us to defining our admissible class as and the variational problem under scrutiny throughout this paper is This problem (P f ) was directly addressed by two of the authors in [44], where expressions for the first and second order optimality conditions were provided. We need to recall them, to motivate and contextualise our results: if we consider u 0 ∈ A and an admissible perturbation h 0 at u 0 (by "admissible perturbation" we refer to the fact that h 0 belongs to the tangent cone to the set A at u 0 . This tangent cone is the set of functions h ∈ L ∞ (Ω) such that, for any sequence of positive real numbers (ε n ) n∈IN decreasing to 0, there exists a sequence of functions (h n ) n∈IN ∈ L ∞ (Ω) IN converging to h as n → +∞, and u 0 + ε n h n ∈ A for every n ∈ IN) then the first order Gâteauxderivative of J T at u 0 in the direction h 0 is where p solves the adjoint equation Here u is the solution of (1) with initial condition u 0 . The main result in [44] states the following: Theorem. [44] There exist a solution u 0 ∈ A of (P f ). Moreover, setting u as the solution of (1) associated with this optimal initial data and p as the unique solution of (5) for u = u, there exists a non-negative real value c such that Finally, for almost every x ∈ {p(0, ·) = c}, one has and the left-hand side belongs to L p loc (Ω). The characterisation of u 0 with the help of p is almost complete here, except on the singular arc ω = {0 < u 0 < 1}. Note first that this singular arc might have a positive measure. It was even proved in [44] that, if f is concave, then ω ≡ Ω. If f ′ is monotonic, equation (6) admits a unique solution and thus fully characterizes u 0 . But for a bistable nonlinearity f θ (u) = u(1 − u)(u − θ), equation (6) might have two solutions, one belonging to [0; η] and the other to (η; 1], where η = η(θ) ∈ (0; 1) is the unique real number such that f θ is convex in [0; η) and concave in (η; 1]. It is then necessary to distinguish between these two possible roots in order to completely characterize u 0 with p. From the numerical point of view, the characterisation given by this result naturally leads to a gradient descent algorithm, which is not well-posed if we are not able to characterise u 0 on ω. Let us briefly describe this algorithm, which we detail further in section 4 of the present paper, to explain the core difficulty and how our theoretical results enable us to bypass it: starting from an initial configuration u 0 0 , we seek to improve it to obtain a better admissible candidate u 1 0 . We first compute the adjoint state p 0 0 associated with u 0 0 . The problem arises if p 0 0 has "flat zones", in other words if there exists c 0 (necessarily unique) such that On the set {p 0 0 > c 0 }, we replace u 0 0 by 1, while on {p 0 0 < c 0 } we replace u 0 0 by zero. On {p 0 0 = c 0 } we must replace u 0 0 with a root of (6). (6) can have two roots µ 0 1 and µ 0 2 . These roots can be distinguished by the convexity of f : up to a relabelling, f ′′ (µ 0 1 ) > 0 and f ′′ (µ 0 2 ) 0. In [44], the two possibilities were explored successively, which led to high computational costs. This was the main limitation of the numerical approach of [44]. Theorem 1 of the present paper shows that one should choose µ 0 2 . This significantly improves the running time of our algorithm and we refer to section 4 for examples.

Related works
A related problem has first been addressed by Garnier, Hamel and Roques in [21], where the authors consider a bistable reaction term f (u) := u(1 − u)(u − θ), with θ ∈ (0, 1), over the full line Ω = IR. In this earlier paper, the authors did not investigate (P f ), but they tried to optimize the initial datum in order to ensure the convergence to u ≡ 1 when t → +∞. They investigated numerically the particular case u 0 := 1 (−α− m 2 ,−α) + 1 (α,α+ m 2 ) , and proved that in some situations the initial datum associated with α = 0 might lead to extinction (that is, u(t, x) → 0 as t → +∞), while a positive α > 0 might lead to persistence (that is, u(t, x) → 1 as t → +∞). Also, numerics for more general classes of initial datum indicate that fragmentation might favor species persistence. Hence, even if the problem we consider here is a bit different, we expect the maximiser to be fragmented, that is, non-smooth, for bistable nonlinearities.
More recently, this problem was also addressed in [32], in a slightly more general form, and for the criterion´Ω |1 − u(T, x)| 2 dx. These authors investigated in particular various conditions ensuring that the maximiser u 0 is constant with respect to x, and, reciprocally, that the constant initial datum is a local maximiser.

Main results of the paper
The main contributions of this paper are the following: • When u 0 does not saturate the L ∞ constraints (i.e. when the set ω := {0 < u 0 < 1} has positive measure), we prove in theorem 1 that any maximiser u 0 must necessarily be in a zone of concavity of f : in ω, f ′′ (u 0 ) 0.
• When f is convex, theorem 2 characterizes explicitely a global maximiser, using rearrangement techniques. Here, the presence of Neumann boundary conditions prohibits using in a straightforward manner the results of [11], and we need to adapt some points of the proof to this case.
• When f is a bistable non-linearity, we improve the algorithm initially introduced in [44] and display several numerical simulations. One-dimensional simulations are displayed that exemplify the fact that theorem 1 significantly improves the computational time of optimisation algorithms. We also provide two-dimensional simulations.

Characterisation of the singular arc
Let us first recall the expression of the second order derivative [44]: where p solves (5) and h solves We can now state our main result: Assume Ω = (0; π). Let u 0 be a solution of (P f ). If the set Ω c := {x ∈ Ω : 0 < u 0 (x) < 1} has a positive measure then, for almost every interior point x of Ω c , there holds f ′′ (u 0 (x)) 0.
Remark 1. The method we put forth is reminiscent of one that was used in [44] to study the case of a constant initial condition and to prove that such a constant u 0 was always a maximiser in the case of monostable non-linearities. Here, working with interior point greatly complexifies the situation and calls for two-scale asymptotic expansions.
The main drawback to our approach is that it can not cover the case of singular (e.g. Cantorlike) singular arcs, and it is a very interesting question to prove that such a property holds for any point of the singular arc Ω c . We comment on the main difficulties of this approach in the conclusion and only state, for the moment, that the main problem is related to the ubiquitous problem of separation of phases in homogenisation [2].

Convex non-linearities and rearrangements
We present, in this section, a characterisation of maximisers when the non-linearity f is convex, using rearrangement and symmetrisation techniques. It should be noted that, since we are working with Neumann boundary conditions, it is not possible to use directly well-known parabolic isoperimetric inequalities [11,42,8]. We refer to [36,11,49] and the references therein for an introduction to parabolic isoperimetric inequalities, and only underline here that the most precise results available in the literature only encompass the case of Dirichlet boundary conditions. For Neumann boundary conditions, a large literature [16,20,33] is devoted to such questions. Usually, it involves a comparison of the solution with the solution of a Dirichlet, or of a mixed Dirichlet-Neumann problem, and makes use of constants appearing in relative isoperimetric inequalities. It is unclear whether these comparison results could be used in our case. Since we are working in the one-dimensional case, a direct adaptation of the proof of [11] yields the required results.
In order to state our result, let us introduce the following notation: letũ 0 be defined as Theorem 2. Assume f is a convex, C 1 function such that f (0) = 0. Thenũ 0 is a solution of (P f ).

Remark 2.
• It should be noted thatũ 0 appears as the solution of many other optimisation problems in population dynamics, most specifically for the monostable case [13,23,38] with Neumann boundary conditions.
This provides an example of non-uniqueness of the maximiser.

Notations, plan of the proof and first simplification
Second order optimality conditions We recall the expression of the second order derivative of J : for an admissible perturbation h 0 , we have where p solves (5) and h solves (8).
Let u 0 be a solution of (P f ). We assume that the set Ω c := {0 < u 0 < 1} has a non-empty interior and we want to prove that f ′′ (u 0 ) 0 almost everywhere on the interior of this set. To do this, we need the following expression of second order optimality conditions: where h is the solution of (8) associated with the initial condition h 0 .
Proof of Lemma 1. We first notice the following thing: let, for any n ∈ IN * , the set F n be defined as Then, for any This is a consequence of the fact that, for any τ such that |τ | is small enough, u 0 + τ h 0 is an admissible initial condition.
Let us now consider h 0 ∈ L ∞ (Ωc) satisfying´Ω h 0 = 0. We define, for any n ∈ IN, For every n ∈ IN, h 0,n is supported in F n and verifies´Ω h 0,n = 0. Hence, defining, for any n ∈ IN, h n as the solution of (8) associated with the initial condition h 0,n we havë However, there holds which, by standard parabolic estimates, entails where h is the solution of (8) associated with the initial condition h 0 . Passing to the limit n → ∞ in (12) yields the conclusion.
Since we want to retrieve, from the second order optimality conditions (11), an information of f ′′ (u 0 ), we need to find a perturbation h 0 such that the ensuing solution h satisfies, roughly speaking, for a certain function g. As h is a solution of a parabolic equation, one possibility to obtain such a behaviour is to choose a highly oscillating initial condition, say h 0 (x) = cos(kx) with a large integer k. This would give h 2 (t, x) ≈ e −tk 2 cos(kx) 2 , which, thanks to the Laplace method, does concentrate around t = 0 up to a proper rescaling. This however is not particularly convenient, as such a perturbation is not admissible: it is not supported in Ωc. To overcome this difficulty, we need to truncate such highly oscillating perturbations, thus choosing a perturbation of the form θ(x) cos(kx), with θ a cut-off function, leading to two-scale asymptotic expansions. The objective is to pick the correct function θ.
In order to summarise our approach, let us fix notations: we pick an optimiser u 0 , we define Ωc := {0 < u 0 < 1}, and we setΩc as the interior of Ωc. To prove Theorem 1 we argue by contradiction: assume that, for some δ > 0, SinceΩc is an open set, we can write it as a union of intervals By (13), there exists n 0 ∈ IN such that so that there exists ǫ > 0 such that, for the same n 0 , we have We fix such an ǫ > 0.
As p(0, ·) > 0 by the parabolic maximum principle, (15) yieldŝ We approximate in L 1 (a n0 ; b n0 ) the function 1 E by a sequence {ψ k } k∈IN of uniformly bounded, non-negative, C ∞ functions that are compactly supported in (a n0 ; b n0 ) ⊂Ωc. In particular, the sequence {ψ 2 k } k∈IN also converges to 1 E in L 1 (Ω), so that (16) implies that for K large enougĥ We fix K large enough so that (17) holds and we set, for this index K, The sequence of truncated, highly oscillating initial conditions is simply ensures that´Ω h k,0 = 0. This constant does not play a role in the upcoming analysis for the following reason: 1. First, by setting h k,0 := θ(x) cos(kx) and by defining h k as the solution of (8) associated with the initial condition h 0 k we shall show thaẗ for some constant C. This is the core of the proof, and will take up the remainder of this section of the paper.
It is also immediate by parabolic regularity to obtain that the sequence {h k } k∈IN is uniformly bounded in L 2 ((0; T ) × Ω).
2. Second we observe that, as θ ∈ C 4 , the Riemann-Lebesgue lemma in particular ensure that 3. If we now set z as the solution of (8) associated with the initial condition θ(·), the solution h k associated with the (admissible) initial condition h k,0 is given by h k = h k + α k z. Then the second order derivative in the admissible direction h k,0 is given by Taking into account (19) and the fact that 4. As a consequence, the theorem is proved, provided we can prove (19), and we henceforth focus on this point.

Asymptotic expansion of h k
Let θ be given as above.We consider the following sequence of equations: let, for any k ∈ IN, h k be the solution of In this context, it is natural [1] to look for a two-scale asymptotic expansion of h k of the form which, after a formal identification at the first and second order, gives the following equations on h 0 k and h 1 k : and Equation (22) can be solved explicitly, giving This, in turn, allows to solve equation (23) as Proposition 1. The asymptotic expansion (21) is valid in L 2 (Ω) in the following sense: there exists M > 0 that depends on the time horizon T such that, if we define then, for any t ∈ (0; T ), In particular,¨( Proof of Proposition 1. To prove this Proposition, we write down the equation satisfied by R k . Straightforward computations show that R k solves and all the functions on the right hand side are evaluated at (k 2 t, x, kx) (we dropped this for notational convenience). We now introduce the following notations: First of all, since 0 u 1 and f ∈ C 1 , there exists M 0 > 0 such that We gather the main estimates on source terms in the following Lemma: Proof of Lemma 2. We prove the three estimates separately. Let us recall the following consequence of the Laplace method: for any integer m ∈ N * , one haŝ Proof of (30) By the triangle inequality we get, for any t ∈ (0; T ), We first use the explicit expressions (24)-(25) for h 0 k and h 1 k to obtain, using θ L ∞ 1, and integrating this inequality between 0 and T giveŝ In the same way, we have, for any t ∈ (0; T ), for some constant C. Taking the square root and integrating in time we get, for a constant C ′ , Using (I m ) with m = 2 givesˆT for some constant M 1 . Summing these two contributions gives (30).
Proof of (31) This follows from the same arguments, by simply observing that Proof of (32) We once again split the expression and estimate separatelŷ dt and 1 kˆT 0 We first observe that for any t ∈ (0; T ), we have In particular, for any t ∈ (0; T ) so that the Laplace method (I m ) with λ = 2 gives the bound for some constant M 2 . The proof of the control of the second term follows along exactly the same lines.
Let us now prove estimate (26). The equation on R k rewrites Multiplying the equation by R k and integrating by parts in space gives In other words, bounding W 0 by M 0 and defining g(t) := R k (t, ·) 2 L 2 (Ω) we obtain Furthermore, R k (0, ·) = 0. We thus obtain, by the Gronwall Lemma, for any t ∈ (0; T ), Hence, by Lemma 2 we get for some constant N 0 and any t ∈ (0; T ),

Back to the proof
We turn back to the proof of Theorem 1 and, more precisely, to the proof of (19).
Proof of Theorem 1. We use the same θ as above and the same notation α k as in the introduction of the proof (equation (18)). Let us now consider the initial perturbation h k,0 := θ(x)(cos(kx) + α k ). We recall that z is the solution of (8) with initial condition θ and that h k is the solution of (8) with initial condition θ(·) cos(k·). Hence, since the equation is linear, Since θ ∈ C 4 , the Riemann-Lebesgue lemma implies We define so that the second order derivative of J T in u 0 rewrites as We focus on the first term: From the assumptions on f and the estimates on u and p, it easy to see that Gathering (40) and (27) it follows thaẗ and similarly, gathering (30) and (26) we obtain Let us now study the term Once again we split the expression. Applying the Cauchy-Swchartz inequality, and using estimates (33)- (34) it follows that the second term verifies The last step in the above expression follows directly from (I m ) with m = 2. We obtain in a similar way the following estimate on the third term: In this case we applied (I m ) for m = 3.
As cos(k·) 2 = 1 2 (1 + cos(2·)) ⇀ k→∞ 1 2 , it follows that for any k large enough G(0) > 0. Furthermore, from the Laplace method, when k → ∞, one has Gathering the estimates in (41), (42), (43), (44), (46) and plugging them into the second derivative of the functional J T given by (38) it follows thaẗ We go back to (38). By (36)-(37) and by (47) we have This means that for k suficiently large, which contradicts the fact that u 0 is a maximiser of J T . The proof of the Theorem is complete.

Proof of Theorem 2
The proof follows essentially from the same arguments as in [11]. We thus only present the main steps that are in order so as to apply the methods of [11]. We define g(u) := f (u) + cu, with c > f ′ L ∞ , so that g is increasing.
Reduction to a bang-bang maximiser. We recall that bang-bang functions are defined as characteristic functions of subsets of Ω, that is, functions only taking values 0 and 1. First, as f is convex, it follows from the same arguments as Proposition 6 of [44] that J T is convex. Hence, one can restrict to maximisers among the extremal points of A. These are exactly bang-bang function: u 0 satisfies u 0 = 0 or 1 almost everywhere on (0, π).
Reduction to a periodic problem. Next, for all t ∈ [0, T ], we extend u(t, ·) to (−π; π) by symmetrisation with respect to 0, and we then extend it to IR by 2π-periodicity. The Neumann boundary conditions at x = 0 and x = π ensure that the extended function is of class C 1 , and it thus satisfies the equation on the torus: Let us also recall some basic facts about rearrangements.

Periodic rearrangements
We recall the definition of the periodic rearrangement: for any periodic function u : T → IR + if we identify T with [−π; π] there exists a unique symmetric (with respect to 0) non-increasing function u ⋆ : T → IR that has the same distribution function as u. u ⋆ is called the periodic rearrangement of u. We recall that the distribution function of u is and that u ⋆ , the periodic rearrangement of u, is the left inverse of µ u .

Proposition 2.
Let v the solution of (49) associated with the initial datum u ⋆ 0 . Then In particular, taking t = T and r = π: As explained in the introduction, Proposition 2 follows simply by adapting minor points in the proofs of [11], and so we will simply indicate the main steps. The core idea is the following: the comparison results and Talenti-type inequalities one finds in the rearrangement literature rely on integrating the solution of the equation we are working with on its level sets and using the isoperimetric inequality. Thus, these proofs generally work only in the case of Dirichlet boundary conditions (for a recent, analogous result in the case of Robin boundary conditions we refer to [?]), as these conditions guarantee that, if the solution is non-negative, none of its level sets intersects the boundary of the domain. However, in our case, since the Neumann boundary conditions allow to symmetrise the solution u and to obtain a solution on the torus T, the boundary of the domain is empty and the isoperimetric inequality holds, so that the proofs are identical.
Proof of Proposition 2. For the sake of readability, we break down the main steps in establishing the inequality • Comparison result for elliptic equations: the first step is to compare the solutions of two elliptic problems. Let ε > 0, let ϕ ∈ L 2 (T) , ϕ 0, let ψ ∈ L 2 (T) , ψ 0 satisfying and let w ϕ , z ψ be the solutions to and Then there holds: To obtain (53), we may follow the standard steps of [50]: we assume that the level sets of (51) have measure zero (to cover the case of level sets of positive measure, one can argue as in [50], to which we refer for the sake of brevity). Let τ > 0 be a real number. Integrating where the last inequality comes from the Hardy-Littlewood inequality. We recall that from the co-area formula for a.e. τ , From the Cauchy-Schwarz inequality and the isoperimetric inequality we obtain Since w ⋆ ϕ is the left inverse of µ wϕ ,standard arguments [8] imply It should be noted that if we work with (52) instead of (51) every inequality becomes an equality since z ψ = z ⋆ ψ , and thus z ψ satisfies Furthermore, integrating (51) and (52) on the torus we obtain, by the equimeasurability of a function and its rearrangement, By the maximum principle, Z ϕ 0, which concludes the proof.
• Comparison result for parabolic equations: for this second step, we follow the strategy of [11], which relies on a Picard iteration scheme. Namely, we discretise the parabolic problem in time: let N ∈ IN * be a discretisation step. For u 0 ∈ A, we define the sequences {u 0,k , v k } k=0 ,... ,N as the solutions to and respectively. As g is convex and increasing, we have g(v) ⋆ = g(v ⋆ ). Thus, we can prove inductively that for every N and for every k ∈ {0 , . . . , N } there holds and it remains to pass to the limit N → ∞ to recover the result.
Conclusion. Assume that u 0 is a bang-bang maximiser of problem (P f ). Symmetrise it and extend it by periodicity. Consider the symmetric decreasing rearrangement u ⋆ 0 of its extension. Then´π −π v(T, x)dx ´π −π u(T, x)dx by Proposition 2, where v is the solution of the periodic equation (49) associated with the initial datum u ⋆ 0 . Clearly, v(t, ·) and u(t, ·) are symmetric with respect to x = 0 for all time t > 0. Hence,´π 0 v(T, x)dx ´π 0 u(T, x)dx. Also, one easily remarks that v restricted to (0, π) is the solution of the parabolic equation with Neumann boundary conditions (1), associated with the initial datum u ⋆ 0 restricted to (0, π). On the other hand, as u 0 is bang-bang, one has u ⋆ 0 = 1 (−m,m) . Hence, 1 (0,m) increases the criterion in (P f ). Thus, it is a solution of (P f ).

Numerical analysis in the bistable framework
As we explained in the introduction, the behaviour of optimisers vary wildly depending on the shape of the reaction term f . To exemplify this phenomenon, we use the bistable non-linearity that motivated [44], namely, f (u) := u(1 − u)(u − θ), with θ ∈ (0, 1).
When considering the optimisation problem (P f ), the fact that the set {p = c} may have a positive measure or, in other words, that a solution may not be the characteristic function of a set, leads to several difficulties in terms of numerical methods, because standard gradient methods or fixed-point algorithms fail to capture what this so-called "singular arc" should be replaced with.
Let us first recall the main principles of the numerical algorithm introduced in [44] and explain the difficulty related to {p = c} further. Given the initial condition at the n-step u n 0 , we construct u n+1 0 = u n 0 + h n 0 , where h n 0 maximises (4) and is an admissible perturbation. Since the adjoint at the n-th step p n 0 may have level sets of positive measure, one can not directly apply the bathtub principle and choose h n 0 as the difference of characteristic functions of two level sets of p n 0 ; we must thus describe what happens on the singular arc, that is, on the level set {p n 0 = c n } where c n is chosen so that |{p n 0 > c n }| < m , |{p n 0 c n }| > m , |{p n 0 = c n }| > 0.
We first define, in this case, u n+1 0 = 1 on {p n 0 > c n }, u n+1 0 = 0 on {p n 0 < c n }, and it remains to fix the value of u n+1 0 on ω n . Defining ω n := {p n 0 = c n } and discretising equation (5) on ω n we obtain, with an explicit finite difference scheme and the value on u n+1 0 on ω n must be a root of (61). However, for bistable non-linearities, this equation may have two roots, say µ n 1 and µ n 2 . In this case, these two roots can be distinguished through the convexity of f . In other words, if we have two roots, up to relabelling, In [44] this difficulty is overcome by examining the two different possibilities and choosing the best one, which significantly lessens the performance of the algorithm, but Theorem 1 allows to overcome this difficulty by choosing directly the root µ n 2 , which is in the "concavity" zone of f .

Comparison of different numerical methods in the one-dimensional case
In this section, we want to study an example in order to compare the performances of our numerical algorithm with other well known optimisation algorithms to solve general nonlinear problems under constraints. More precisely, we will consider the following numerical methods: • Method 1: The numerical algorithm introduced in [44], which we improve using Theorem 1, and that will be referred to as our algorithm.
• Method 2: The interior-point method, which is used to solve optimisation problems with linear equality and inequality constraints by applying the Newton method to a sequence of equality constrained problems. For a more detailed description of this method see for instance [15].
• Method 3: The sequential quadratic programming (SQP), which solves a sequence of optimisation sub-problems, each of which optimizes a quadratic model of the objective function subject to a linearisation of the constraints, see for instance [46].
• Method 4: The simulated annealing method, which is a probabilistic technique used to approximate global optimisation in a large search space. See for instance [28] for more details on this technique.
We used the MATLAB platform to perform the simulations. Methods 2 and 3 are already coded in the MATLAB function "fmincon" while methods 1 and 4 were coded for the experiment.
Setting the data Let us consider Ω = (−50; 50), and m = 13; thus the admissible set is defined as follows: Note that this set is defined by two inequalities and an equality constraint. We aim at maximising the quantity J T (u 0 ) :=´Ω u(T, x)dx for T = 25 and we use a bistable reaction term f (u) := u(1 − u)(u − 0.25).
In order to compare the performance of the four algorithms under the same conditions, we consider the same discretisation of Ω. Moreover, the solution of the equation is systematically computed by the Crank-Nicolson method, and, for the initialisation we consider the same u 0 0 given by a single block of mass 13. The value of the objective function at each iteration is numerically approximated by the rectangle rule. In particular, for the initialisation we have J 25 (u 0 0 ) = 29.42. The results of the simulations are shown in Fig. 1 and Table 1. For this example, our algorithm turns out to be faster than other well-known algorithms. Moreover, the evaluation of the objective function differs in less than 1% with respect to the best result obtained with the sequential quadratic programming method which takes more than twice the run-time of our algorithm. Though the solution given by the sequential quadratic programming method is clearly more regular than the others, the profile of the local optimisers found by simulated annealing and by our algorithm do not seem to be far from this profile. Indeed, the solutions obtained through Methods 1, 3 and 4 are qualitatively similar. On the other hand, the interior-point method gives a significantly different optimum, which seems to point out the good performance of our algorithm. It is important, however, to highlight that since uniqueness is not guaranteed in general, one can not ensure that the algorithms have converged to a global maximiser but only to a local one.

Numerical simulations in the two-dimensional case
We only considered in the present paper the one-dimensional case. We now display some numerical results obtained in dimension 2, for which new patterns might arise.
To solve the reaction-diffusion equation in the two-dimensional case, we consider the alternating direction implicit method (ADI) which is a classical method to solve parabolic problems in two or three dimensions. As in the one-dimensional case, the algorithm and routines were coded in MATLAB.
We consider a square domain Ω = (−10; 10) × (−10; 10), discretised uniformly by squares of side dx = 0.22 . We fix T = 30 for all subsequent simulations. We first tackle the case of a bistable reaction term In a second paragraph, we study the monostable case The justification for this second case is that this is a non-concave monostable non-linearity. It is hence not covered by the theoretical results of [44].

The bistable case
Example 1 The algorithm is initialised with a ball of full density located in the middle of the domain Ω. The mass is fixed to m = 5.8π, see Fig. 2(a). After 20 iterations, the algorithm converges to the local optimum showed in Fig. 2(b). The evolution of the objective function J 30 through iterations is showed in Fig. 3 (a) u 0 0 (b) u 0 = u 20 0 Figure 2: In the left-hand side, we show the input of the algorithm, given by the ball of radius r = √ 5.8 centered at the origin. In the right-hand side, we display the local optimum found by the numerical algorithm after 20 iterations, which looks radial, but is no longer a bang-bang distribution: it does not only take values 1 and 0. One might see that the local optimum found by the numerical algorithm is no longer a bangbang function but a circular ball with less mass in the middle and a slightly bigger ratio. Looking at the adjoint state defined as the solution of equation (5), associated to this initial data, one might see that the area in the middle of the circle corresponds to a set where the adjoint state remains constant, see Fig. 4. Figure 4: The figure shows the surface given by the solution p(0, x) of the adjoint problem defined by the eq. 5 associated to the initial data u 0 found by our algorithm. The plane colored in gray, is associated to the value c described in Theorem 1.2.1 and thus for every x ∈ Ω such that p 0 (x) = c, one has 0 < u 0 (x) < 1, see Fig. 2

(b).
Example 2 In this case, we keep the same discretization and initial mass m of the previous example, but we consider an initial data which is a stripe of full density dividing our domain into two equal regions of zero density, see Fig. 5 (a). The algorithm converges after 38 iterations and the local optimum is displayed in Fig. 5 (b). The corresponding variations of the objective function is showed in Fig. 6 (a) u 0 0 (b) u 0 = u 38 0 Figure 5: In the left-hand side is showed the input of the algorithm, given by the stripe of width r = 0.91 centred at the origin. In the right-hand side, the local optimum found by the numerical algorithm after 38 iterations. We observe that, in this case, the value of the objective function remains very low during the first 20 iterations. This fact, together with the radial geometry of the optimum found by the algorithm suggests that this stripe geometry is not optimal. It should also be pointed out that the geometry of the local optimum is interesting: indeed, it shows regions of zero density (i.e. the optimum u 0 found by the algorithm is equal to 0 in theses regions) in the middle of regions of full density (i.e. where u 0 = 1), which exemplifies the phenomenon described in the one-dimensional case in [21].
Another relevant feature is that the optima found in the first and second examples are different, which indicates that our algorithm converge to local optima, and thus that the choice of the initial distribution u 0 0 is crucial.

Example 3
For this example we keep the settings of the previous one, but we consider a higher initial mass m = 27. The geometry of the initial distribution is a stripe of full density dividing the domain into two regions of zero density, like in the example 2, see Fig. 7(a). The corresponding local optimum found by the numerical algorithm is showed in Fig. 7(b). As in the previous case, the optimiser reflects a low density zone ringed by a high density region. This gap is clearly filled by diffusion as time evolves.
This example suggests once again the non optimality of stripe-like initial distributions. Note from Fig. 8 that, despite the considerable increase of the initial mass with respect to example 2, the values of the objective function J 30 (u 0 0 ) associated to the stripe is of the order of 10 −5 , which is very low compared with the value associated to the final distribution.
Finally, let us mention that a possible approximation of the maximiser was discussed in the Appendix of [34]. Namely, in this thesis, the author replaced the maximiser u 0 by its mean on each of the connected components of {u 0 > 0}. This gives pretty good results in several cases. It would be good to manage to quantify analytically the difference of criterion between this approximated initial datum and the global maximiser.

The non-concave monostable case
We now present some simulations in the case f (u) = (u + 0.25)u(1 − u),  still working under the constraint that 0 u 0 1. The motivation behind this case is that this non-linearity is monostable on (0; 1), i.e. it only has one stable equilibrium, but it is not concave. As a consequence, the theoretical approach developed in [44] can not guarantee that the optimiser u 0 is the characteristic function of a subset of Ω.
The parameters of the simulations are still the same: T = 30, and the initial configuration is the same as in Example 1: the initialisation is a ball of full density, with mass m = 5.8π located in the middle of the domain Ω. We refer to Fig. 9 and 10. We should however point out that in this simulation, the value of the objective function remains almost constant, despite the fact that the final configuration is very different from the initial one. It is plausible that such monostable non-linearities converge too quickly to the equilibrium.

Conclusion, open problems and possible extensions
We make, in this conclusion, several concluding remarks and comments about possible generalisations and extensions of the results presented in this paper. For each of them, we try to present the arguments that have led us to the conclusion that other approaches were necessary in general.

Regarding the regularity of the singular arc
One of the main drawbacks of Theorem 1 is the regularity assumption on the singular arc. Namely, we obtain, for a maximiser u 0 of J T , the characterisation f ′′ (u 0 ) 0 only on the interior of the singular arc {0 < u 0 < 1}. The method presented in this paper (two-scale expansions) strongly relies on the smoothness properties of the cut-off function θ.
It may be tempting, in the one dimensional case, to overcome this difficulty arguing as in [37]: if we simply assume that the singular arc ω is measurable, but has positive measure, one can show that for every K ∈ IN there exists h K ∈ L 2 (Ω), supported in ω, that writes h K = k K α k,K cos(k·) , h K L 2 = 1.
Using the fact that h K only has high Fourier modes, one may hop for a two-scale expansion of the form This is however a priori prohibited by the problem of separation of phase: to obtain such a description, one needs welll-separated phases, in the sense of [2]. In the context of Fourier series, this would require, at the very least, that h K should write as a lacunary Fourier series (typically, h K = ∞ j=0 α j,K cos(K j x)). However, for such lacunary Fourier series, Zygmund's theorem (see [25] for instance) prohibits that they have compact support, so that admissible perturbations can not have this structure. This is a major drawback, and it is unclear whether or not one may be able to overcome this difficulty via a similar approach, or if an entirely new strategy needs to be devised.
Another approach would be to prove some regularity on ω ensuring that almost every of its point lie in its interior. This is satisfied for example if u 0 is Riemann integrable (since, due to Lebesgue's characterisation of Riemann integrable functions, almost every point is a continuity point of u 0 ). Riemann integrability is satisfied by BV functions. Unfortunately, we were not able to push the regularity further than L ∞ .

Monostable non-linearities
As seen in subsection 4.2.2 of this paper, the numerical approach we propose, based on Theorem 1, works for general monostable non-linearities. The theoretical tools are, however, not sufficient at this level to fully characterise optimisers. An interesting question would be to discuss whether or not optimisers in the monostable case are always bang-bang, are if some degeneracy zones can appear.

The singular arc in higher dimensions
It may be plausible to adapt the methods of Theorem 1 to obtain a characterisation of the singular arc analogous to that of Theorem 1 in the case Ω = N i=1 [0; a i ], a i > 0. To do so, the main difference with our proof would be to replace the initial perturbation θ(x) cos(kx) with N i=1 θ(x) cos(kx i ).

Rearrangement inequalities for other types of boundary conditions
In this work, we mostly dealt with the case of Neumann boundary conditions in the one-dimensional case. We ought to note two things: first, the proof of theorem 1 should hold in the case of Dirichlet or of Robin boundary conditions, provided the functions cos(k·), in the proof, are replaced with the Dirichlet or Robin eigenfunctions of the laplacian in the interval. Second, regarding theorem 2, the same type of results can be obtained in a straightforward manner for the case of Dirichlet boundary conditions, by applying directly [11]. The case of Robin boundary conditions may be encompassed by using the recent Talenti inequalities obtained in this case in [6]. Addressing the problem on the full line IR is more tricky, since in this case even the existence of a maximiser is unclear. We plan on investigating such matters in future works.