Numerical reconstruction based on Carleman estimates of a source term in a reaction-diffusion equation.

In this article, we consider a reaction-diﬀusion equation where the reaction term is given by a cubic function and we are interested in the numerical reconstruction of the time-independent part of the source term from measurements of the solution. For this identiﬁcation problem, we present an iterative algorithm based on Carleman estimates which consists of minimizing at each iteration cost functionals which are strongly convex on bounded sets. Despite the nonlinear nature of the problem, we prove that our method globally converges and the convergence speed evaluated in weighted norm is linear. In the last part of the paper, we illustrate the eﬀectiveness of our method with several numerical reconstructions in dimension one or two.


Introduction
Let Ω be a C 2 bounded domain of R d for d = 1, 2 or 3 and T > 0. We consider the following reaction-diffusion equation    ∂ t u(t, x) − ∆u(t, x) + u 3 (t, x) = σ(x)h(t, x), (t, x) ∈ (0, T ) × Ω, u(t, x) = g(t, x), (t, x) ∈ (0, T ) × ∂Ω, u(0, x) = u • (x), x ∈ Ω, (1.1) where g is the Dirichlet boundary data and u • is the initial condition. In the right hand side of the first equation, we assume that the time-dependent function h is known and we focus on the reconstruction of σ which is assumed to depend only on the spatial variable. To identify this unknown, we have two kinds of measurements, the flux of the solution on a part of a boundary and the solution in the whole domain at a given time: m(t, x) := ∇u(t, x) · n(x), (t, x) ∈ (0, T ) × Γ, r(x) := u(T 0 , x), x ∈ Ω, (1.2) where Γ ⊂ ∂Ω, T 0 ∈ (0, T ) and n is the outward-pointing unit normal vector defined on ∂Ω.
Regarding the applications, this model can represent the evolution of a pollutant in the atmosphere. The source in the right hand side corresponds to a spill of pollutant and we want to localise it. This model can also be viewed as a simplified model to represent the evolution of the electrical potential in the heart (we refer to [8] for a detailed presentation of this application domain and more precisely to Section 2.9.7 of [8] for cubic-like reaction models). In our model, the natural propagation of the potential is initiated by the initial condition and the source in the right hand side may correspond to a secondary undesirable source that we want to identify.
Let σ max > 0 be a fixed constant. We assume that the source term σ that we want to reconstruct belongs to L ∞ (Ω) and satisfies the following a priori bound: For this problem, according to Bukhgeim-Klibanov method, σ is uniquely determined by the measurements and a Lipschitz stability estimate holds under appropriate assumptions on the data (the precise result is stated in Prop. 2.6). Bukhgeim-Klibanov method [5] is a classical theoretical method to prove the uniqueness and stability for parameter identification problems. For a presentation of this method which relies on Carleman estimates and for a survey on its applications, we refer to [18] and in particular to Section 3.3 on parabolic equations. For the inverse problem of coefficients identification in nonlinear parabolic equations, let us in particular mention that [3,9] deal with the theoretical stability of the reaction term in a semi-linear PDE.
In this paper, our aim is to tackle the numerical reconstruction of σ and to propose for this nonlinear problem a globally convergent method. Our work is drawn from a numerical method presented in [1] for the identification of a potential in a wave equation. The method strongly relies on Carleman inequalities and it consists of an iterative algorithm minimizing at each iteration a cost functional involving Carleman weights. The main strength of this numerical method is that it globally converges to the exact solution i.e. it converges independently of the initialization. In particular, contrary to classical minimization techniques like Tikhonov methods [16], it is not necessary to add a priori knowledge on the source term through the data of a background state to convexify the cost functional.
As pointed out in the introduction of [2], this method induces several numerical challenges. In particular, the classical Carleman weights have very strong variations due to the presence of a double exponential involving large coefficients. That is why, as in [2] for the wave equations, we need to construct new Carleman weights for the heat equation which involve single exponentials (these weights are given by (2.3) and (2.4)).
The presence of a nonlinearity in our PDE leads to additional difficulties in the study of the numerical methods. In particular, the strong convexity properties of the cost functional are restricted to bounded spaces. Moreover, the operator appearing in the cost functional has to be modified by adding truncation operators in the nonlinear terms to tackle these terms in the proof of the convergence of the numerical method. At last, contrary to [2] where the PDE is linear, introducing a conjugate variable e sϕ z does not allow to overcome the fact that, even with single exponentials, the minimization of the cost functional is challenging.
Let us mention that we could have considered general cubic functions of the form a 3 u 3 + a 2 u 2 + a 1 u with a 3 > 0 instead of a simple cubic monomial in this semi-linear parabolic partial differential equations (1.1). By this way, the study includes the bistable equation or Allen-Cahn equation. Moreover, we refer to Remark 2.4 for some remarks on the case of other boundary data (Neumann boundary conditions instead of the second equation in (1.1) and boundary measurements on the solution itself).
Using Carleman estimates to solve numerically inverse problems has been first considered in the paper [17] by Klibanov. This method called convexification method has been applied in several papers. In [19], the authors are interested by the reconstruction of a coefficient in a parabolic equation and present a gradient method applied to a strictly convex cost functional involving Carleman weights. In [22], the authors consider the reconstruction of the initial condition in a nonlinear parabolic equation. We also refer to [20] for the most recent paper which applies this convexification method.
For numerical studies applying Carleman estimates to controllability problems, we refer to [7] for the numerical controllability of the wave equation and [12] for the numerical controllability of the heat, Stokes and Navier-Stokes equation.
The paper is organized as follows. Sections 2 give some preliminary results. First, in Section 2.1, we present a Carleman estimate for the heat operator with Dirichlet boundary conditions. In this estimate, we consider two kinds of Carleman weights: the classical weights for the heat equation with a double exponential and new weights involving single exponentials which are introduced for numerical purposes. Then, in Section 2.2, we state a regularity result satisfied by the solution of equation (1.1). The proofs of the Carleman estimate and the regularity result are presented in Appendices A and B respectively. At last, in Section 2.3, we state the stability inequality associated to our inverse problem. Section 3.1 is the core of the paper and presents the numerical reconstruction method of the source term. The latter is an iterative process which requires at each iteration the minimization of a functional based on the Carleman estimate. This section states the global convergence of the method (Theorem 3.4). In Section 3.3, we establish properties satisfied by the functional to minimize at each step. In particular, the existence of a global minimizer of the functional is stated in Lemma 3.7 and the strong convexity on bounded set is proved in Lemma 3.8. In the last section (Sect. 3.4), we prove the global convergence property. Finally, Section 4 is devoted to the implementation of the algorithm and the numerical results obtained for several 1D and 2D test cases.

Carleman inequality for the heat equation
Without loss of generality, from now on, we assume that T 0 = T 2 . In this section, we state a Carleman inequality for the heat equation in two cases. The first case corresponds to the classical weights with a double exponential while, in the second case, the weights only involve single exponentials as in Section 3 of [27]. Let us specify these two cases: -Case 1: For λ > 0, we define θ and ϕ by: for all (t, x) ∈ (0, T ) × Ω where η 0 satisfies the following properties: and where x 0 is an arbitrary point in R d \ Ω and ρ is a constant satisfying 0 < ρ < 1. We notice that θ > 0 and ψ < 0. In this case, we assume in addition that x 0 and Γ are such that Let us mention that the spatial part ψ of the Carleman weight in Case 2 resembles the one proposed in [1] for the wave equation. Moreover, the geometric condition (2.5) which is classical for the wave equation (see [14,23]) is unusual for the heat equation and is linked to this new choice of weights. With these weights, we have less flexibility in the computations and we need an extra condition on the measurement domain compared to the classical weights corresponding to Case 1. On the other hand, if we take the classical weights, the presence of a double exponential in the functional to minimize (see (3.11)) is prohibitive to address numerical applications (we refer to Rem. 2.3 for additional comments). In all our numerical tests presented in Section 4.2, we have considered the weights given by Case 2.
Let us now formulate the Carleman inequality in Case 1 and Case 2.
Here and in all the paper, we denote by C a positive constant which depends on T and Ω, λ in Case 1 and ρ in Case 2, unless specified otherwise where appropriated. The proof of this theorem is given in Appendix A. A consequence of Theorem 2.1 is the following lemma: Under the same assumptions as Theorem 2.1, there exist s 0 > 0 and C > 0 such that, for all s ≥ s 0 : 0 Ω e 2sϕ (s 2 θ + 2s|∂ t ϕ|)|z| 2 dxdt 0 Ω e 2sϕ s 3 θ 2 |z| 2 dxdt where we have used that |∂ t ϕ| ≤ Cθ 2 . Thus, the result follows from (2.6).

Remark 2.3.
To better design Carleman weights for numerical purposes, it would be interesting to make a comprehensive comparison between different possible choices of Carleman weights for the heat equation. In particular, in such a study which is beyond the scope of our paper, it would be necessary to spell the lower bound on s in the associated Carleman inequality.
Remark 2.4. We could have considered other kinds of boundary data by completing the first equation of (1.1) with Neumann conditions instead of Dirichlet conditions and by replacing the first measurement in (1.2) by a measurement on a part of the boundary of u itself. In this case, following [11,13], we still have a Carleman inequality with the classical weights (2.1) and we can still prove the global convergence of the numerical method. For the numerical tests, it would be interesting to see if we can get a Carleman inequality with weights similar to the ones of Case 2.

Regularity result
Let us give a regularity result for problem (1.1). The proof of this result is presented in Appendix B.
Let us note that, in the above proposition, the regularity assumed for g is not optimal, it would indeed be sufficient to assume that g ∈ H 1 (0, T ; H 3/2 (∂Ω)) ∩ H 2 (0, T ; H κ (∂Ω)) with κ > 0 (see [24], Chap. 1, Sect. 9.2). In this result, if we do not make the assumption that h(0, ·) = 0 in Ω, it is necessary to assume that σ belongs to H 1 (Ω) (since we need an initial condition in H 1 (Ω) for the problem satisfied by ∂ t u). But this additional regularity assumption on σ leads to difficulties in the construction of the iterations in Algorithm 3.2.

Stability inequality
In this paragraph, we state a Lipschitz stability inequality for our inverse problem. This result asserts in particular that the unknown σ is identifiable from the measurements given by (1.2). It is obtained thanks to a direct application of Bukhgeim-Klibanov method [5] and relies on the Carleman inequality given by Theorem 2.1 and the regularity result given by Proposition 2.5. We do not give the proof here and refer to [15] for a closely related result. Proposition 2.6. We assume that u • ∈ H 3 (Ω), g ∈ H 1 (0, T ; H 3/2 (∂Ω)) ∩ H 2 (0, T ; H 1/2 (∂Ω)) and h ∈ H 1 (0, T ; L ∞ (Ω)) is such that h(0, ·) = 0inΩ and |h(T 0 , ·)| ≥ β > 0inΩ. We consider σ 1 and σ 2 in L ∞ (Ω) which satisfy (1.3). Then, for i = 1, 2, if we denote by u i the solution of (1.1) associated to σ i , we have the following inequality: there exists C > 0 such that

Presentation of the algorithm and convergence
In this subsection, we construct a sequence (σ k ) k∈N which approximates the unknown σ and we state the convergence of this sequence. We make the following assumptions: and -The weights θ and ϕ are given by Case 1 or Case 2 described at the beginning of Section 2.1. In Case 1, the parameter λ is fixed and large enough.
First, we initialize the sequence with σ 0 = 0 (or any guess such that σ 0 L ∞ (Ω) ≤ σ max ). Now, let us assume that we are at step k and that we have constructed σ k which satisfies We denote by u k the solution of (1.1) associated to σ k and by u σ the solution of (1.1) associated to the unknown σ. Moreover, we set v k = u σ − u k . We then use Proposition 2.5 and we denote by M > 0 a fixed constant depending on T , Ω, σ max , u • H 3 (Ω) , h H 1 (0,T ;L 2 (Ω)) and g H 1 (0,T ;H 3/2 (∂Ω))∩H 2 (0,T ;H 1/2 (∂Ω)) such that v k C([0,T ]×Ω) + v k C 1 (0,T ;H 1 (Ω)) + v k H 2 (0,T ;L 2 (Ω)) + v k H 1 (0,T ;H 2 (Ω)) ≤ M . (3.4) The function v k is solution of where we have set q[v, u] = 3u 2 + 3uv + v 2 . Let us differentiate the equation with respect to time. We introduce w k = ∂ t v k which satisfies: x ∈ Ω, (3.6) where, for all (t, x) ∈ (0, T ) × Ω, (3.7) Let us now explain the core idea of the numerical method that we will introduce below. We notice in (3.5) that x ∈ Ω. (3.8) Hence, if w k (T 0 , ·) was known, then σ could be directly computed, because we assume that h(T 0 , ·) satisfies (3.2) and the other terms in (3.8) are given observations thanks to (1.2). However, since f k defined in (3.7) depends on σ, w k is unknown. Thus, the idea is to use Z k obtained via the minimization step (Step 2) of Algorithm 3.2 as a proxy for w k . In Section 3.4, we will estimate the discrepancy between Z k and w k (both seen as minimizers of functionals) with respect to f k . For the constant M > 0 introduced in estimate (3.4), we consider the following function: The properties satisfied by T M are given in Section 3.2. For any µ in L 2 ((0, T ) × Γ) and for s large enough, we introduce the functional J 0,k [µ] by where By this way, since v k satisfies (3.4), T M (v k ) = v k and P k (w k ) corresponds to the left hand side of the first equation of (3.6). We consider the functional J 0,k [µ] on the function space , endowed with its natural norm. The next iteration σ k+1 is defined by following four steps: where m is the measurement defined in (1.2) and u k is the solution of (1.1) associated to σ k . (3.14) • Step 4 -At last, we define where Π σmax is given by Remark 3.3. Let us give some comments regarding the different steps of Algorithm 3.2.
-According to Lemma 3.7, J 0,k [µ k ] (defined in (3.11)) admits a global minimizer in E if s is large enough. Hence, this result in particular ensures the existence of a critical point which is needed in Step 2 of the algorithm. For the numerical implementation of this step, we refer to Remark 4.2. Let us also notice that, as the functional J 0,k [µ k ], Z k depends on s but we drop this dependence to avoid heavy notations. -In Step 3, σ k+1 is well-defined because h satisfies the positivity condition (3.2) and h(T 0 , ·) belongs to L 2 (Ω). Moreover, in this expression, v k (T 0 , ·) is known and given by v k (T 0 , ·) = r − u k (T 0 , ·) where r is the measurement defined in (1.2). Since u k (T 0 ) and v k (T 0 ) belong to H 2 (Ω) and Z k (T 0 ) belongs to L 2 (Ω), σ k+1 belongs to L 2 (Ω). -Step 4 is needed to ensure that σ k+1 satisfies (3.3) at step k + 1. Now we state the main theoretical result which gives the global linear convergence in the weighted L 2 -norm of the sequence (σ k ) k∈N : Thus, for s large enough, (σ k ) k∈N tends to σ when k goes to +∞.
This theorem will be proved in Section 3.4.
Remark 3.5. Let us notice that our method may also be applied to the identification of a source in the simpler case of the linear heat equation. In this case, the inverse problem is linear and thus the properties of our method (in particular the global convergence) are much more classical. If we simply consider the least square method, the functional is quadratic and, thanks to the stability estimate, we can prove that this functional is strongly convex. Thus, a classical gradient descent method globally converges and it is not necessary to introduce our algorithm which is more complex.

Properties satisfied by the function T M
Proposition 3.6. The function T M defined by (3.9) belongs to C 2 0 (R) satisfies the following properties: a) For all X ∈ R, where χ A is the characteristic function of a set A. c) For all X 1 , X 2 ∈ R, which implies in particular that T M is a Lipschitz operator. d) There exists C > 0 such that, for all X 1 , X 2 ∈ R, b) By definition (3.9), T M (X) = 0 for all |X| ≥ 2M . Moreover, for all |X| ≤ 2M , we have c) This is a direct consequence of (3.17) and the mean value inequality. d) By the same arguments as in b), we show that Hence, (3.19) is a direct consequence of (3.20) and the mean value inequality.

Properties satisfied by J 0,k [µ]
The following lemma ensures the existence of a global minimizer of J 0,k [µ] in E. This result in particular ensures the existence of a critical point Z k introduced at Step 2 in Algorithm 3.2.
Proof. According to Proposition 2.5, there exists a constant M > 0 such that (3.21) Using this estimate and (3.18), we have the continuity of J 0,k [µ] in E. Moreover, since J 0,k [µ] is positive, it admits an infimum in E and we can introduce a sequence (z n ) n∈N such that Let us study the convergence properties of the sequence (z n ) n∈N . First, we notice that we can write P k z n under the form where r k,n and s k,n only depend on u k , ∂ t u k and T M (y n ). Using inequality (3.21) and the property (3.16), we deduce that r k,n is bounded in L ∞ ((0, T ) × Ω) and s k,n is bounded in L 2 ((0, T ) × Ω) by some constant M . Hence, writing that According to the Carleman inequality given by (2.6) and using the fact that se 2sϕ θ ≤ C in (0, T ) × Ω for the third term in the right hand side, we deduce that, for s large enough, (3.22) By construction of (z n ) n∈N , the sequence (J 0,k [µ](z n )) n∈N is bounded and thus the left hand side of this last inequality is bounded. According to the definitions of θ and ϕ which are given by (2.1) or by (2.4), we have in (0, T ) × Ω |∂ t θ| + |∂ t ϕ| ≤ Cθ 2 and |∇θ| + |∇ϕ| + |D 2 θ| + |D 2 ϕ| ≤ Cθ and thus (e sϕ θ −1/2 z n ) n∈N is bounded in H 1 (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)). We deduce that, (e sϕ θ −1/2 z n ) n∈N weakly converges to some element in H 1 (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)) that we denote e sϕ θ −1/2z (all the convergence results given in this proof are valid up to a subsequence but we do not specify it in order to lighten the writing).
Moreover, since H 1 (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)) is compactly embedded in L 2 ((0, T ) × Ω), and, by identification of the limit, since θ −1 belongs to L ∞ ((0, T ) × Ω), we also have e sϕ θ 3/2 z n e sϕ θ 3/2z weakly in L 2 ((0, T ) × Ω) and e sϕ θ 1/2 ∇z n e sϕ θ 1/2 ∇z weakly in L 2 (0, T ; H 1 (Ω)). (3.24) Let us now prove that lim and (e sϕ P k z n ) n∈N weakly converges in L 2 ((0, T ) × Ω) and it is sufficient to identify their weak limits. The fact that directly comes from (3.24). To identify the limit of (e sϕ P k z n ) n∈N , we will prove that We first consider in the definition (3.12) of P k the three first terms which correspond to the linear part. The weak convergence of ( Now we define, for all t ∈ (0, T ) For the other terms in the operator P k , let us first prove that (e sϕ θ −1/2 T M (y n )) n∈N strongly converges to By definition (2.1) or (2.4) of ϕ and θ and since T 0 = T 2 , we have, for all t between T 0 and t, for all x ∈ Ω This implies that Thus, according to (3.23), (e sϕ θ −1/2 y n ) n∈N strongly converges to e sϕ θ −1/2ỹ in L ∞ (0, T ; L 2 (Ω)) and since T M satisfies (3.18), this implies that We can now study the limit of the remaining terms of e sϕ P k z n when n tends to +∞ : using (3.16), (3.21) and (3.28), we have and Let us now prove that The strong convergence (3.23) of (e sϕ θ −1/2 z n ) n∈N implies the almost everywhere convergence of (z n ) n∈N tõ z and the existence of a function z b in L 2 ((0, T ) × Ω) such that, for all n ∈ N Moreover, the strong convergence of (e sϕ θ −1/2 y n ) n∈N in L 2 ((0, T ) × Ω) implies the almost everywhere convergence of (y n ) n∈N toỹ. Thus, we deduce that And these two properties imply (3.31) according to Lebesgue's dominated convergence theorem. At last, we use the same arguments to prove that (3.32) Finally, gathering (3.26) and (3.29) to (3.32), we obtain (3.25) and we conclude thatz is a minimizer of Due to the nonlinearities in our equation, we can not state the strong convexity of J 0,k [µ] in E. Nevertheless, it is interesting to notice that the strong convexity property holds if we consider a smaller space than E including some boundedness hypotheses which allow to deal with the nonlinearities (similar results are obtained in [19])). We state this property in the following lemma: This Lemma is proved in Appendix C. Contrary to the wave equation where the weights stay far from 0 (see [1], Sect. 4), our weights, as usual for the heat equation, vanish at 0 and T and it is not clear that a minimizer of J 0,k [µ] in E will belong to E C for some C > 0. Therefore, it is not possible to deduce the uniqueness in E from the strong convexity in E C .
This convexity property is not used in the proof of Theorem 3.4 but it is an important property for the convergence of numerical minimization methods (we refer to Rem. 4.2 for a discussion on the numerical methods and the fact that the property shown here does not correspond exactly to the numerical framework).

Proof of the convergence result stated in Theorem 3.4
For µ k = ∂ t (m − ∇u k · n) on (0, T ) × Γ, we define the functional where P k is given by (3.12) and f k is defined by (3.7). We notice that w k , solution of the equation (3.6) minimizes Let us now compute the Gâteaux derivative of P k at point w, for any w ∈ E. Let z ∈ E, Then, w k satisfies the first order optimality condition given by Similarly, Z k satisfies the first order optimality condition Let us define z k = w k − Z k . We compute the difference between (3.35) and (3.36) and take z = z k . We get For what follows, we define, for all t ∈ (0, T ) We will estimate separately the different terms of this inequality. We divide the computations in several steps.
-Step 1. Let us first find a lower bound for the first term in the left-hand side of (3.38).
Using (3.16), (3.18) and (3.21), we can estimate R 1,k Moreover, from (3.34), we can write DP k (Z k )(z k ) = ∂ t z k − ∆z k + R 2,k where, according to (3.16) and (3.21) We deduce from these inequalities that According to (3.17) and using that Z k = w k − z k and that we get for the last two terms that ≤ M e sϕ y k 2 L ∞ (0,T ;L 2 (Ω)) .
We have, for all t ∈ (0, T ) using inequality (C.5) for ϕ. By this way, we deduce that Thus, using again (3.43), we get (3.45) -Step 3. To bound the last term of (3.38), we notice that

This implies that
(3.47) Using Theorem 2.1, we can eliminate the last term in the right hand-side of (3.47) for s larger than some constant s 0 . Thus, using inequality (2.7), we get the following bound on z k (T 0 ): In the left hand-side of this inequality, we have z k (T 0 , x) = w k (T 0 , x) − Z k (T 0 , x) for x ∈ Ω and, using (3.8) and (3.14), we get that In the right hand-side of (3.48), since f k = (σ k − σ)∂ t h and h is assumed to be bounded in H 1 (0, T ; L ∞ (Ω)), we have Using ( Thus, for s large enough we deduce from this inequality the convergence of the sequence (σ k ) k∈N to σ. This concludes the proof of Theorem 3.4.

Numerical methods
In this subsection, we present the discretization procedure and the numerical methods used in our numerical simulations. To simplify the presentation, we explain the discretization scheme in the one-dimensional case and assume that Ω = (0, L) for L > 0 and Γ = {x = L}.

Generation of the data
In this article, we work with synthetic data. To discretize the reaction-diffusion equation (1.1) for the exact source σ, we use a finite differences scheme based on the three-point backward Euler scheme and a linearization of the cubic term. We denote by N x ∈ N the number of discretization points in the interior of [0, L] and by N t ∈ N the number of discretization points in the interior of [0, T ]. The space and time steps are denoted by ∆x = L N x + 1 and ∆t = T N t + 1 respectively and we define, for 0 ≤ j ≤ N x + 1 and 0 ≤ n ≤ N t + 1, u n j a numerical approximation of the solution u(t n , x j ) with t n = n∆t and x j = j∆x. The approximated solution is computed in the following way: For 0 ≤ n ≤ N t , knowing u n , compute u n+1 as the solution of the linear system: where the time implicit cubic term (u n+1 j ) 3 has been approximated by its first order Taylor expansion (u n j ) 3 + 3(u n j ) 2 (u n+1 j − u n j ). Then, we compute the counterpart of the continuous measurements r and m given in (1.2) as follows: with n 0 is the integer part of N t /2 + 1.
On the computed data, we may add a Gaussian noise: where N (0, 1) satisfies a centered normal law with deviation 1 and α is the level of noise (i.e. α = 0.01 corresponds to a noise of 1%).

Discrete algorithm
We present in this subsection the discrete version of Algorithm 3.2. In order to lighten the notations, we will denote byσ the source term at the current iteration (previously denoted by σ k in the continuous framework).
Algorithm 4.1. Initialisation: Start withσ = 0. Iteration : Until the convergence criteria is reached, do and set v j = r j −ū n0 j .
The iterative loop is stopped when two consecutiveσ are closer than a fixed relative tolerance ε or when the maximal number of iterations is reached. In the absence of knowledge of the exact solution σ, the quality of the converged solution is measured thanks to the following criteria that should be of the order of the noise level on the observations. If the exact solution σ is known, we can also compute the relative error Remark 4.2. For the numerical implementation of Step 2 in Algorithm 4.1, we determine a critical point of J 0,k [µ] by applying the Newton method or Newton-Krylov method [4,21]. This last method belongs to the family of inexact Newton methods and consists of solving at each step the linear system in a Krylov subspace. The Newton minimization method is globally convergent if the functional is strongly convex. To be in this framework, it would be necessary to prove that the discrete functional J 0,k (4.5) is strongly convex under boundedness assumptions (like in Lemma 3.8 for the functional in the continuous setting) and to prove that the discrete minimization sequence satisfies these bounds. A complete study of the discretized algorithm could be tackled in the future and would involve in particular a Carleman inequality for the discretized heat equation. In our numerical simulations, we have taken the initial guess of the iterative Newton-Krylov method equal to 0 and checked that the convergence does not depend on this initialization. Thus, we observe numerically global convergence properties for the minimization of J 0,k [µ].

Remark 4.3.
In order to avoid the inverse crime, we introduce a bias by taking different schemes for the direct and the inverse problems. Hence, we solve (4.1) associated to σ thanks to a linearized implicit scheme and we use an explicit scheme for the nonlinear term in equation (4.3) withσ = σ k .

Numerical challenges
One of the main drawbacks of the numerical method presented in Algorithm 3.2 is that we have to differentiate in time the observation m in (4.4) and to take the Laplacian of the observation r in (4.7). Thus, even a small perturbation (noise) on the observations may induce a large perturbation on its derivatives. In order to partially remedy this problem in the presence of noise, we first regularize the data (m, r) thanks to a 3-order low-pass Butterworth filter [6] associated to a cutoff frequency ω. We also replace the classical finite difference formulae in (4.4) and (4.7) that generate instabilities by a Savitzki-Golay formula [25] associated with a cubic polynomial and a window size of 5 points.
As already mentioned previously, another difficulty is the presence of the exponential weights in the functional that leads to severe numerical difficulties when performing the minimization for s large. Indeed, to ensure the strong convexity of the functional J 0,k (see Lem. 3.8) and the convergence of Algorithm 3.2 (see Thm. 3.4), s has to be large. In [2], this difficulty was solved by choosing a functional that only depended on the conjugate variable e sϕ z and the corresponding conjugate operator. But this was possible because the considered operator was linear. Here, we managed to deal with this difficulty by introducing the new weight functions (2.4). In Figure 1, we plot e sϕ in (0, T ) × (0, L) for s = 1 and s = 100. Notice that even for s large, the function does not vanish at the observation time T 0 = 0.5 what allows a good reconstruction of the source term in the whole domain Ω. Numerically, we observe that for s = 1, the minimisation step is slow (5202 seconds for s = 1 versus 17 seconds for s = 100 for the test case of Fig. 3a) and in some cases the convergence of the algorithm is not achieved (for example in the test case of Fig. 3b).

Numerical results
This subsection is devoted to the presentation of some numerical examples to illustrate the properties of the numerical reconstruction method and its efficiency. All simulations are executed with Python. The source codes are available on request. Table 1 gathers the numerical values used for all the following examples, unless specified otherwise where appropriate. Moreover, we construct the function Φ introduced in (3.10) in the form: Figure 2 presents some examples of data generated by the direct problem. In all the figures presenting the numerical results, the exact source that we want to recover is plotted by a red line, whereas the numerical source recovered by our method is represented by a dotted black line. The convergence informations (number of iterations, running time, convergence errors) are reported in Table 2.

Simulations from data without noise
In Figure 3, we present the successive results obtained at each iteration of Algorithm 3.2 in the case of the reconstruction of the source σ(x) = sin(πx) for two different choices of h. One can observe that in both cases the convergence criteria (4.8) for ε is achieved in less than 20 iterations. In Figure 4, several results of reconstruction of sources obtained using Algorithm 3.2 in the absence of noise are given.   16 554 0.7% 0.1% 0.8% Figure 5a 3 87 1% 0.3% 2% Figure 5b 3 91 1% 0.3% 4% Figure 5c 3 97 3% 0.5% 9% Figure 6b 4 497 0.05% 0.1% 0.05% Figure 6d 6 802 0.1% 0.1% 0.01% Figure 5 shows the results for σ(x) = sin(πx) with different levels of noise in the measurements (α = 1%, 2% and 5%). In Table 2, we report the corresponding errors on the reconstructed source. In fact, we observe that a noise of level α in the measurements gives rise to an error of order 2α in the recovered source.

Simulations in two dimensions
We also performed some reconstructions in two dimensions where Ω = (0, 1) 2 , x 0 = (−0.3, −0.3) and Γ = {0} × [0, 1] ∪ [0, 1] × {0} . By this way, assumption (2.5) is satisfied. Figure 6 presents the results obtained    If we are in Case 1 with the first choice of weight (2.1), this result is proved in an identical way as Lemma 1.2 in [13] which considers the case of internal measurements. Assume now that we are in Case 2 where θ and ψ are given by (2.4).
Let us give some properties on ϕ which will be useful in what follows: In the proof, we assume that z belongs to C 2 ([0, T ] × Ω) and satisfies z = 0 on (0, T ) × ∂Ω. A density argument allows to come back to the regularity hypotheses of the theorem.
For all s > 0, we set w = e sϕ z and we introduce the conjugate operator Q defined by If we set f = ∂ t z − ∆z, we have Qw = e sϕ f.

Some computations give
where the operators Q + and Q − are defined by In a classical way, we write that The main part of the proof consists of bounding from below the terms in the right hand side by positive and dominant terms and a negative observation term located in (0, T ) × Γ. For the sake of clarity, we divide the proof in several steps.
-Step 1 -Explicit calculation of the cross-term. We set where I i,k is the integral of the product of the ith-term in Q + w and the kth-term in Q − w.
Integrations by parts in time give easily since w(0) = w(T ) = 0 in Ω and w = 0 on (0, T ) × ∂Ω. An integration by parts in time gives for I 21 We compute in the same way, by integrating by parts in space ∇w · n(∇ϕ · ∇w) dγdt and Since ∆ϕ is independent of x, by integration by parts in space, we have At last, we have Gathering all these computations and using the second property in (A.1), we get ∇w · n(∇ϕ · ∇w) dγdt.
(A.7) For the second term in the right hand side, we notice that the main part in s 3 is given by For the boundary terms in (A.7), we notice that, since z = 0 on (0, T ) × ∂Ω, ∇w = e sϕ ∇z. In particular, ∇w · τ = 0 on (0, T ) × ∂Ω. Thus, we get We divide this last integral as follows According to (2.5), the second integral is positive and the first integral corresponds to an observation integral. Gathering these estimates, (A.7) becomes, for s large enough Step 2 -Bounds on ∆w and ∂ t w.
From the definition of Q − (A.5), we have In the same way, Thus, coming back to (A.6) and, gathering (A.8) and these last two estimates, we get, for s large enough Step 3 -Back to the variable z.
Appendix B. Proof of the regularity result given by Proposition 2.5 We split the proof in several steps.
where the coefficients (α im ) 1≤i≤m being to be determined by the conditions: From Picard-Lindelöf theorem (see, for example [26]), the system (B.5)-(B.6) of nonlinear ordinary differential equations, admits a unique local in time solution (α im ) 1≤i≤m in C 1 defined on a maximal interval (0, T m ). -Step 3 -A priori estimates.
Next, let us consider w = ∂ t u which is, according to (1.1) and (3.1), formally the solution of in Ω.
Appendix C. Proof of Lemma 3.8 On E, we consider the norm · s defined by For any fixed C > 0, the set E C is convex and closed in ( E, · s ). Let z 1 and z 2 be given in E C . We set z = z 1 − z 2 . Then, we have 0 Ω e 2sϕ P k z 2 DP k (z 2 )(z) dxdt − s T 0 Γ e 2sϕ θ(∇z 2 · n − µ)(∇z · n) dγdt = T 0 Ω e 2sϕ (P k z 1 − P k z 2 )DP k (z 2 )(z) dxdt + T 0 Ω e 2sϕ P k z 1 (DP k (z 1 )(z) − DP k (z 2 )(z)) dxdt + s T 0 Γ e 2sϕ θ|∇z · n| 2 dγdt. (C.1) To estimate the two first terms, we will follow similar computations as in Section 3.4 since we can notice that these terms also appear in (3.37) if we replace z 1 , z 2 and z respectively by w k , Z k and z k . Since the assumptions on the functions are different (in Section 3.4, w k = ∂ t v k where v k satisfies (3.4) whereas here z 1 and z 2 belong to E C ), we detail the arguments below when they are different from the ones in Section 3.4. For the first term, we follow the computations made in Step 1 of Section 3.4 and the counterpart of (3.39) is: 0 Ω e 2sϕ |y| 2 (|∂ t u k | 2 + |z 1 | 2 + |z 2 | 2 )dxdt.