IDENTIFICATION OF NONLINEARITIES IN TRANSPORT-DIFFUSION MODELS OF CROWDED MOTION

. The aim of this paper is to formulate a class of inverse problems of particular relevance in crowded motion, namely the simultaneous identiﬁca- tion of entropies and mobilities. We study a model case of this class, which is the identiﬁcation from ﬂux-based measurements in a stationary setup. This leads to an inverse problem for a nonlinear transport-diﬀusion model, where boundary values and possibly an external potential can be varied. In spe- ciﬁc settings we provide a detailed theory for the forward map and an adjoint problem useful in the analysis and numerical solution. We further verify the simultaneous identiﬁability of the nonlinearities and present several numerical tests yielding further insight into the way variations in boundary values and external potential aﬀect the quality of reconstructions.

The mathematical modeling of many-particle systems is of central importance for many areas of science, e.g. physics, biology and increasingly also in socio-economics. The classical approach of passing from microscopic models (at the single particle level) such as Newton's equations of motion to macroscopic partial differential equations relies on closure relations like Boltzmann's "Stosszahlansatz" (cf. [5]), which can be made rigorous under certain scaling assumptions on particle number N and radius σ (cf. [11]). A classical example is the Boltzmann-grad limit N σ 2 tending to a constant, which applies well to dilute gases (cf. [12]). The crowded case N σ 3 tending to a constant on the other hand is much more subtle and the macroscopic limit is not unique. Using different assumptions and closure relations (in some cases even different microscopic models), a variety of macroscopic equations have been derived. In a diffusive scaling as thought of in particular for biological applications, a canonical form is where u = u(x, t) denotes the particle density, D = D(u) a nonlinear mobility, E = E(u) an entropy functional, V = V (x) an external potential and W = W (x) a symmetric interaction kernel to model long-range forces such as electrostatics or chemotaxis (W * u denoting the convolution of W with the density). Such equations are (metric) gradient flows of the entropy-energy functional In the case of no-flux boundary conditions stationary solutions are characterized as stationary points of E subject to non-negativity and a total mass constraint. Several special cases of (1) have been proposed in different models for crowded motion (cf [8,57,15,52,10,56,9]), with different assumptions and approximations difficult to be verified. We thus propose a different approach in this paper, namely to identify the model components directly from available data related to the particular applications. A variety of data is collected already in the situations mentioned earlier, e.g. video surveillance, traffic counts, time-resolved microscopy, or electrophysiological measurements (cf. [29,54,14,40]). So far these data are mainly used for statistical evaluations not coupled to appropriate macroscopic models, but only phenomenological descriptions like anomalous (linear) diffusion based on analyzing root-mean-square distance over time. Here we propose to use such data for inverse problems related to (1). Although inverse problems are of interest for the potential V (mainly control) and interaction kernels (identification) as well, we focus on the two nonlinearities D and E appearing from the limiting procedures (note that V and W can be translated directly from microscopic forces and could thus be considered also in inverse problems for the particle system). In the classical limit without crowding, the mobility is linear and the entropy is logarithmic, i.e.
In the crowded case several modifications have been proposed, e.g. additional quadratic terms in E (short-range pair repulsion, cf. [48]) or even more nonlinear ones (cf. [58,3]). In particular saturating mobilities D tending to zero as u reaches a maximal density value (cf. [28,57]) have been of interest for many years. As one may imagine from the multitude of applications, there are indeed several different mathematical expressions of the identification problems to be solved. Distinctions concern the type of measurements (density-or flux-related), the design of experiments (single or multiple, stationary or transient situations), prior information (non-negativity and concavity of D, convexity of E), and the use of information about the stochastic nature of the limit (e.g. in the likelihood functional). There is plenty of literature on the identification of nonlinear diffusion coefficients, either from distributed measurements (cf. [25,17]) or from boundary measurements (cf. [51,41] for the parabolic and [27,43] for the elliptic case). There are also some works about the reconstruction of the nonlinear flux term in conservation laws, cf. [34,33]. However, the simultaneous reconstruction of nonlinearities has, to the best of our knowledge, not been considered in the literature so far. The paper closest to our approach appears to be [26], where they identify two coefficients simultaneously, but not both of them are nonlinearities.
This paper is organized as follows: In Section 2 we will discuss the modeling, measurements, and possible formulations of inverse problems. Section 3 is devoted to the properties of the direct problem. In the remaining part of the paper we shall focus on a particular case, namely the identification from boundary flux measurements in a stationary situation. We consider the uniqueness of the inverse problem as well as differentiability of the parameter-to-output map (Section 4), provide an explicit reconstruction formula using linearizations (Section 4.3). In Section 5 we present a least-squares approach to the fully nonlinear problem and conclude by presenting the results of numerical experiments (Section 6).

2.1.
Motivation of the Inverse Problem. In this Section we further detail the identification problems for interacting particle systems and their continuum limit. We illustrate the continuum limit process for pair interacting Brownian motions, a commonly used model in biological applications, and comment on the particular challenges of this process for crowded setups. The motion of pairwise interacting particles, which undergo Brownian motion, is given by where X i (t) ∈ R d denotes the position of the i-th particle at time t, V ext = V ext (x) is a potential due to external forces, and V int = V int (x) an interaction potential (e.g. Coulomb or gravitational forces in physics or attractive-repulsive forces in animal group motion). The noise terms W i denote independent Wiener processes. The dynamics of the N -particle system can be characterized by considering the evolution of the joint probability density which can be computed from the Liouville-Fokker-Planck equation In order to obtain feasible macroscopic descriptions, the single-particle density is used. Assuming that the initial states are equal and independent, i.e.
and that the standard propagation of chaos property holds, i.e.
the mean field (Vlasov) limit of the process is of the form This derivation can be made rigorous in terms of the BBGKY-hierarchy (consisting of all k-particle distribution functions) under suitable assumptions like regularity of V int and hence the continuum limit is rather clear. Identification issues in this aspect mainly concern the external and interaction potentials, which might correspond to unknown interaction rules (cf. [45,20]) and the diffusion coefficient σ.
In the crowded setup two additional challenging issues arise in the derivation of the above macroscopic limit, namely: • In order to treat the limit of infinitely many particles with finite sizes rescaling has to be carried out in the interaction distance, typically of the form ). • The typical interaction potentials are not regular, they even might have singularities. Often even hard-core potentials (zero outside and infinite inside the particle volume) are used to enforce non-overlapping of particles.
The first issue results in nonlinear diffusion terms in the equation in the limit N → ∞, e.g. into where E is a local entropy functional, cf. [49,50,48]. However, even in the most natural case of hydrodynamic scaling (α = 1) the exact form of the limiting entropy E cannot be characterized, it can only be shown to be bounded above by a quadratic function (cf. [49,50]). The second issue, i.e. the addition of singular interaction potentials, can destroy the propagation of chaos property (if particles become strongly correlated locally). The limit rather becomes with a correlation term C N that cannot be obtained explicitly. Different assumptions and approximations have to be used to obtain the limiting equations of the form (1). Typical measurements collected for such particle systems (related inherently to Lagrangian respectively Eulerian viewpoints) are: • Tracking: Here M single-particle trajectories X i (t) are recorded over a time interval [0, T ], but often M N . This procedure is routinely carried out in cell biology using fluorescence and phase-contrast microscopy (cf. [60,46,14]), it seems to become of increasing importance in behavioral biology and human crowd motion based on video surveillance and other methods (cf. [29,54,14,40]). • Flux Measurements: Here the number of particles passing a codimension one manifold per unit of time is recorded. This is the standard procedure in car and pedestrian traffic counting, where the manifold is an imaginary line on the street or pathway, as well as in electrical measurements, where it corresponds to the electrodes recording charged particles.
A standard approach in model calibration for tracking data is to use log-likelihood estimation by minimizing or a penalized version subject to the unknown parameters. This is based on the interpretation of u = u(x, t) as the probability density of a particle being at x at time t. This is well justified asymptotically in the standard mean-field limit due to (5), which implies that the log-likelihood of the M particles indeed tends to the sum of the independent log-likelihoods. The latter is highly questionable in the crowded case if C N in (7) does not tend to 1. We leave this inherent modeling question (or the development of alternative approaches) as a challenge to be considered and hopefully clarified in the future. We finally mention that due to limited resolution in the imaging processes (single particle trajectories not to be analyzed easily), it might also be an interesting route to directly consider macroscopic measurements of the density u or the flow density where U ∈ R corresponds to a possibly variable strength. This quantity can be obtained from image analysis techniques such as optical flow estimation (cf. [44,30,6,53]).

2.2.
Inverse Problems from Total Flux Measurements. In this paper we shall consider the inverse problem based on macroscopic flux data, i.e. the expected value of the flux. Apparently such macroscopic data can be treated in a more straight forward manner than trajectories. We will restrict ourselves to applications without interaction potential, i.e. we assume W = 0 in the following. The macroscopic limit of the expectation of the flux over a manifold Γ can be carried out as the limit in the partial differential equation and yields where n denotes the outward unit vector and with j as in (9). In this paper we consider the inverse problem in the special setup of a stationary solution. In this case u is the solution of with Ω ⊂ R d , d = 1, 2, 3 and Dirichlet boundary conditions u = f on Γ D ⊂ ∂Ω and no-flux boundary conditions on ∂Ω \ Γ D . As above, the external potential is of the form U V , with variable strength U ∈ R. The no-flux boundary conditions correspond to impenetrable walls while the Dirichlet conditions prescribe the density of particles at entrance and exit. From a modeling point of view it would also be reasonable to prescribe a constant flux at the entrance and a flux proportional to the density at the exits. However, since we focus on the theoretical treatment of the problem in this paper, we shall prescribe Dirichlet conditions on the whole boundary from now on, namely Here, 0 < f < 1 is meant pointwise since Ω is at most of dimension three and thus f is continuous by Sobolev embedding. Hence, the inverse problem is to find D and E from measurements over a part of the boundary Γ ⊂ ∂Ω given by for a class of boundary conditions f and potentials U V . In most instances we shall consider the one-dimensional version of this inverse problem, which is wellmotivated by applications like ion channels and traffic. In this case Ω corresponds to an interval, Γ D to its end points with Γ being one of them.

Basic Properties of the Stationary Model
In this Section, we shall examine basic properties of our model given by where Ω is bounded and V = V (x) is a given potential. For simplicity we shall rewrite (14) as We remark that no information is lost since E can be recovered from G once D is known and E is anyway only determined up to a constant. The flow density (9) in normal direction then reads as where n is the outward normal to ∂Ω. As mentioned in Section 2, we shall for simplicity supplement (14) with Dirichlet conditions of the form (12). Furthermore, we assume that the domain is such that and introduce the nomenclature: We denote by PS(G, D; U, V, f ) the problem of solving (15) for given G, D, U, V and subject to the boundary conditions (12) for the unknown u. In one space dimension, we shall write PS(G, D; U, V, u L , u R ), where u L , u R ∈ R + denote the left and right boundary value, respectively.
The following assumptions on the nonlinearities D, E and V will remain valid throughout the whole paper: Here I = [0, 1] is the interval between zero and the maximal density, normalized to 1. Due to the continuity of D and G, assumption (A2) implies G(u) ≥ 0 and D(u) ≥ 0 for u ∈ [0, 1]. In some Sections we will need additional assumptions on D, G. These will be stated explicitly.
3.1. Transformation to Entropy Variables. We consider the entropy functional It is easy to verify that, if u is a solution of the transient version of (14) with no-flux boundary conditions, this functional decreases in time. The entropy functional is also important in the stationary case as it allows to transform (15) into entropy variables ϕ given by Due to (A2), the function E is invertible for 0 < u < 1 and we can solve this equation for u. Using (15) this yields An important feature of (18) is that it satisfies a maximum principle which, in original variables, is only true if ∆V = 0, i.e. if the potential is harmonic.

3.2.
Existence in the case: V harmonic. For a given potential with ∆V = 0, (15) features a maximum principle which implies 0 < u < 1 due to (12). Then G > 0 and therefore H ∈ C 1 (Ω) defined by H = G is strictly monotone and H −1 exists and is also in C 1 (by means of the implicit function theorem). This allows us to define w = H(u) and F (w) = D(H −1 (w)) with the new variable w satisfying the elliptic equation Since H maps the interval [0, 1] to [a, b] for some constants a, b with 0 = H(0) < a ≤ b < H(1) < ∞, it is meaningful to supplement (19) with the Dirichlet boundary condition We shall prove existence of solutions of (19), (20) using a fixed-point argument. We consider the linearized equation (21) ∆w = U F (w)∇w · ∇V, withw given. Note that F is well-defined here since D ∈ C 1 (Ω). The fact that this equation features an obvious comparison principle (and its solution is therefore uniformly bounded in L ∞ (Ω)) motivates the definition of the set We consider the solution operator to equation (21) w → w, w solution to (21), (20), where I H 2 →L ∞ denotes the embedding operator. Using Schauder's fixed point theorem (cf. [55]) we thus conclude: Theorem 3.2 (Existence in the case of a harmonic potential). Assume d = 1, 2, 3, ∆V = 0 and furthermore D, G such that F (w) = (D(H −1 (w))) is continuous and bounded for all 0 < a ≤ w ≤ b < 1. Then, there exists a solution w to (19), (20) in Proof. As noticed above the existence of a solution follows from Schauder's theorem and moreover the argument implies that w = H(u) ∈ H 2 (Ω). Since H has a continuously differentiable inverse and w is bounded, we easily infer that ∇u = (H −1 ) (w)∇w ∈ H 1 (Ω). Since the transform between u and w is one-to-one from , it suffices to verify uniqueness of the solution w of (19) with given Dirichlet boundary condition. This follows for U sufficiently small since the linearized differential operator in (21) is uniformly elliptic on H(M) ∩ H 1 (Ω).

3.3.
Existence in the case: ∆V = 0. Next we present an existence result in the case of general nonlinear potentials in one space dimension. The above strategy fails as equation (15) no longer satisfies a comparison principle. However, a maximum principle still holds for the equation in entropy variables (see Section 3.1), i.e.
with boundary conditions To prove that there exists a solution to this equation, we again construct a fixedpoint argument, this time in the set which is well defined and continuous due to assumptions (A1)-(A3). Next we consider the linearized equation with a given coefficient b ∈ L ∞ ((0, 1)) and define the corresponding solution oper-atorR This operator is well-defined (cf. [21]) and its continuity can be shown using standard estimates, cf. [10,Lemma 4.2]. Using the fact that the embedding of H 1 ((0, 1)) into L ∞ ((0, 1)) is compact we conclude that the operator maps N into N and that the set R(N ) is precompact. Using again Schauder's fixed-point theorem we conclude that: Since this results has been obtained by Schauder's theorem we cannot deduce uniqueness in general. However, for U small, we can construct a contractive fixedpoint operator: Proof. By adding and subtracting the term Due to our assumptions, (D •(E ) −1 ) positive and thus it is clear that the linearized elliptic equation with givenṽ ∈ L 2 ((0, 1)), has a unique solution. If we chose U sufficiently small and exploit the fact that under our assumptions D•(E ) −1 and K −1 are continuous, the operator solving this equation is a contraction on H 1 ((0, 1)) and thus, by Banach's fixed-point theorem, we conclude the existence of a unique solution.

Existence and Uniqueness for a Linearized Equation.
In this Subsection we assume again ∆V = 0 and consider the linearized equation where 0 < u < 1 is solution of (15), g ∈ H −1 (Ω), H = G and supplemented with Introducing v =ṽ − v b , we obtain the equation x ∈ Ω and with homogeneous boundary conditions. This can be rewritten as and hence we define with H −1 defined as the dual of H 1 0 . The formal adjoint to equation (30) in the scalar product is given by with an arbitraryF ∈ L 2 (Ω). Dividing by G(u) > 0, using the transformation ϕ =φ − ϕ b and applying ∆ −1 we obtain the equation with homogeneous Dirichlet boundary conditions. The operator is continuous and compact (due to the boundedness of D , G > 0 and ∇V ∈ L ∞ guaranteed by assumptions (A1) and (A3)). As the solutions to the homogeneous version of (32) is unique due to the maximum principle, we obtain a unique solution ϕ for every F ∈ L 2 (Ω) by Riesz' Theorem as in Section 3.2. Now, the Fredholm alternative asserts that the solution of (29) is unique. Thus the existence of a solution v can again be obtained by Riesz' Theorem, since the operator is obviously continuous and compact as well. Furthermore, note that the operator norm U is proportional to U . Thus, for U sufficiently small we can apply a perturbation argument for linear operators, cf. [37, Thm. 1.16, Remark 1.17], to the equation where ∆V = 0, I denotes the identity, and conclude the existence of a unique solution. To summarize, we have Lemma 3.5. Let u be a solutions of (15), (12) in the sense of Theorem 3.2. Then, there exists for every g ∈ H −1 (Ω) a solutionṽ ∈ H 1 (Ω) of (27), (28). The solution depends continuously on g and is unique if U is sufficiently small.
Proof. The existence ofṽ follows from the above considerations. The fact that the solution is in H 1 (Ω) can be easily obtained exploiting the fact that G(u) > 0. The continuity follows from standard arguments, see [22,Chapter 8].

Inverse problem
This Section is devoted to inverse problems related to the nonlinear convection diffusion equation (15). We are interested in identifying the functions D and G from flux measurements in different situations.

4.2.
Differentiability. First we would like to study the differentiability of the parameter-to-solution map. Note that throughout this Section we use the letter C as well as C 1 , C 2 , . . . for generic, not necessarily equal constants. We impose the following additional assumption (B1) G ∈ C 1 (I) (i.e. E ∈ C 3 (I)), then we have: Condition (2): The Fréchet derivative with respect to u is given by We start again by examining the continuity with respect to G (for D similar arguments hold) As before, this estimate yields the continuity with respect to G and G as h → 0.
Using the same arguments we obtain where C depends on G C 1 (Ω) , D C 1 (Ω) and U ∇V L 2 (Ω) . Conditions (1) In these cases we are able to deduce explicit solutions for the identification of D and G, which serve as a starting point and validation for the nonlinear problems. These linearizations can be motivated from different applications, like ion channel or crowd motion models. First we start with the following assumptions being reasonable for flux measurements in ion channels: let u R = u L and V (x) = x. Then the equilibrium solution for u = const satisfies Therefore we can identify D(u) from given flux measurements, i.e.
We conclude that by using the definition of c. Furthermore we assume that u R is a small perturbation of u L , i.e. u R = u L + and obtain We conclude that it is possible to identify the functions D and G in situations close to equilibrium. Equations (42) and (43) are not valid for non-equilibrium situations, but they serve as a starting point for fully nonlinear identification problems.

Remark 1.
We would like to mention a different application, namely crowd motion models. Here equation (15) includes an additional geometry information via an area function a(x). The function a(x) corresponds to the "width" of the domain at a point x. Then the one dimensional reduction of (15) reads as A common assumption in crowd motion modeling is that the velocity is constant, i.e. the potential is linear. Here the constant U corresponds to the typical walking speed of a human which is around 1.3 m/s, cf. [59]. Then one can consider the same linearization approach as above, but it is not possible to obtain an explicit reconstruction formula for G anymore.

4.4.
Identifiability. From now on, we restrict ourselves to one spatial dimension and the domain Ω = [0, 1] and assume that the potential V (x) is harmonic. Next we shall prove identifiability of D and G in the fully nonlinear case. We will use the adjoint method that has been introduced in [16]. The key ingredient is an integral identity relating the nonlinearities to solutions of an adjoint problem defined below.
Here we need the additional assumption that D is concave, i.e. (B2) D (u) ≤ 0. This assumption is natural since all standard examples as well as many analytical results require the concavity of D.
Remark 2 (Flux in One Dimension). We mention that in spatial dimension one, the flux in the stationary equation is constant in space, hence it does not matter whether we speak about total flux or flux density, and we will simply use the nomenclature of flux j in the following. Let us also comment on the rule of thumb about dimensions of data and reconstructed quantities: We can change both boundary values separately, yielding the flux as a function of two variables, i.e. two-dimensional measurements. Since we aim at identifying two functions of a single variable, the problem might be overdetermined. On the other hand there may be correlations between different settings of boundary values. However, for U = 0 we obtain at least independent measurements if we fix one boundary value, e.g. at 1 2 and vary the other. Hence, the effective dimensionality of the data is certainly not smaller than the dimensionality of the unknowns. Overdetermination should further hold if we vary U in addition.
Remark 3 (Identifiability of D). The strategy to identify D by applying equal Dirichlet boundary values used in (42) in Section 4.3 works for the fully nonlinear problem considered here. Therefore we will always assume that D is known in the following.  For given (D, G 1 ) and (D, G 2 ) consider the corresponding solutions u 1 , u 2 of (15), (12) with associated fluxes j 1 , j 2 . Then the following relation holds: where λ is a solution of and with Proof. We subtract equation (15) with (u 1 ; D, G 1 ) and (u 2 ; D, G 2 ) from each other and obtain For H i with H i = G i , i = 1, 2 we rewrite the equation as Multiplying with λ and integrating by parts yields to Note that the last line is zero, due to (46a) and that the boundary term (H 2 (u 1 ) − H 2 (u 2 )) λ x | x=1 x=0 , vanishes because u 1 and u 2 satisfy the same Dirichlet boundary conditions. Imposing boundary conditions (46b) on λ we obtain (45).

Remark 4.
Since p, q ∈ L ∞ (Ω), p > 0 in Ω and V x ∈ L ∞ (Ω), the existence of a unique solution λ ∈ H 1 (Ω) follows from standard results for linear elliptic equations, see again [22,Theorem 8.3 and Theorem 8.12]. Furthermore, since for identifiability we always need to assume j 1 = j 2 , the term (j 1 − j 2 ) will always be zero in the proceeding proofs. Thus, we are free to chose arbitrary Dirichlet boundary conditions for λ. Here we apply (46b) to ensure that (46a) has a non-trivial solution with specific properties, see below.
We continue by proving the followings results about the sign of u x and λ x : is equivalent to solve two Cauchy problems on the intervals I 1 = [0, x 0 ] and I 2 = [x 0 , 1], both with boundary conditions u(x 0 ) = c, u x (x 0 ) = 0. Then, in one of the intervals I 1 , I 2 there exists a pointx such that u(x) < u(x 0 ). However, Hopf's Lemma applied to this interval asserts that u x (x 0 ) is either strictly positive or strictly negative which is a contradiction. Since by smallness of U we know that u is unique, this completes the proof.
This strategy also applies to the linear adjoint equation: Proof. The proof is the same as in Lemma 4.3. Since Equation (46a) is linear it features a maximum principle and we conclude that there exists only one solution.
Before we state the main identification result we would like to introduce the notion of distinguishability. In the following we denote by u i a solution obtained by solving (15) with G = G i and Dirichlet boundary conditions (u L , u R ) for i = 1, 2. For each u i we use j i for the corresponding flux. Proof. If G 1 and G 2 would be distinguishable, there would exist a partition of I into intervals I k G , k = 1, . . . , N such that either G 1 > G 2 or G 2 < G 1 on each interval. In each interval I k G , we can choose boundary values u L and u R such that the values of u 1 lie in this interval (due to the maximum principle). Then we have Since due to Lemma 4.3 and 4.4, both (u 1 ) x and λ x have a sign, the integrand is either strictly positive (G 1 (u 1 ) > G 2 (u 1 )) or negative (G 1 (u 1 ) < G 2 (u 1 )). Thus, we obtain a contradiction and G 1 = G 2 on I k G . In the case that G 1 = G 2 on an interval, there is nothing to show.
We remark that the above Theorem contains the implicit assumption that there exist only two different functions G 1 , G 2 which would create the same output data. Of course, this is not true in general and there could be an infinite number of different functions G i creating the same data. This also explains why the above Theorem contradicts the rule of thumb that the dimension of the measured data should be the same as of object one wants to identify. However, the finite measurements depend on the specific G 1 , G 2 . Given an infinite number of measurements, with all boundary values in an interval, Theorem 4.6 can be applied to an arbitrary pair and thus the final identifiability result is given by the following corollary: Corollary 1. Given a set of infinitely many flux measurements I k (generated by infinitely many combinations of U , u L and u R ), it is possible to identify the function G, up to distinguishability, on the interval I.

Least Squares Approach and Regularization for the Fully Nonlinear Problem
Next we consider the fully nonlinear inverse problem. Here we would like to solve the following optimization problem : (47) min Here I δ denotes the noisy measurements and PO the observation operator defined in (39). We remark that the above integral is well defined as PO is continuously differentiable with respect to U , u L and u R . This can easily be shown using the same arguments as in the proof of Theorem 4.2. Since the potential V is assumed to be given, we will from now on neglect it when writing the arguments of PO to increase readability. Furthermore, µ = µ(U, u L , u R ) is a measure defined on the set R × I 2 and is chosen according to the type of application. A typical example would be a given potential strength U 0 and a finite number of boundary values u k L , u k R , k = 1, . . . , M . Then, µ is given by Note that (47) might be ill posed, hence we add two Tikhonov regularization terms (see [19]) and obtain J(D, G) = 1 2

Numerical Solution
Finally we shall discuss the numerical discretization of the forward problem (15) as well as the numerical simulation of the inverse solver. Note that G(u) can become small compared to D(u) and therefore the numerical solver should be capable of dealing with convection dominated problems. The proposed method has been introduced by Egger and Schöberl for linear convection diffusion problems, see [18]. Furthermore we discuss a Kaczmarz type method for the efficient solution of the inverse problem. Finally we present numerical results based on the linearization approach and the fully nonlinear problem.
6.1. Mixed HDG. We solve equation (15) using Newton's method and a hybrid mixed discontinuous Galerkin method. Introducing the new variable σ = −∇u we obtain We discretize the mixed system using a hybrid discontinuous Galerkin method, i.e.
for all test functions (q, v, µ) ∈ Q h × V h × M h . Here T denotes the elements, ∂T the element interfaces and E h the set of all facets. The finite elements spaces are given by where RT k denotes the Raviart Thomas finite element space (up to order k) and P k are polynomials up to order k. The unknown λ corresponds to the value of the function u at the element interfaces. The function u + denotes the upwind value of u, given by Note that the variational formulation is stated for space dimensions d = 1, 2, 3, although we only consider space dimension one in this manuscript. Then the elements correspond to intervals T = [x i , x i+1 ], i = 1, .....N and the interfaces to the grid points ∂T = {x i , x i+1 }.
The above system defines a nonlinear operator equation F (σ, u, λ) = 0, which can be solved using Newton's method, i.e. JF (w n )r n = −F (w n ) w n+1 = w n − αr n , where w n = (σ n , u n , λ n ) and JF denotes the Jacobi Matrix of the operator F . We solve the adjoint equation (51) using a the similar discretization as for the forward problem. The adjoint equation (51) is linear, hence we can apply the mixed hybrid DG scheme as presented in [18].
6.2. Kaczmarz methods. In most applications several flux measurements for various boundary conditions are available. Hence we choose a Landweber Kaczmarz iteration (cf. [36,42]) to solve the inverse problem efficiently. We assume that we have p measurements available, therefore the cost functional is given by We update the functions D and G using the following iterative procedure where u k,j and w k,j are solutions of We assume that the nonlinearities D and G are in H 1 (Ω), hence the can be written as where ϕ l denotes the classical hat functions in H 1 . Therefore the updates (56) and (57) correspond to the updates of the coefficients D l and G l .

Results
Finally we present first numerical results for the reconstruction of nonlinearities from flux measurements, which are based on the linearized as well as the fully nonlinear problem. These simulations confirm the feasibility of the proposed objectives and serve as a starting point for the more sophisticated reconstruction techniques in higher space dimension.   Table 1. Boundary conditions for the linearized case We generate the synthetic data by solving (15) on an equidistant mesh of size h f = 1 2 h with the Newton solver presented in Section 6.1. We used the same discretization for I in the data generation and reconstruction routines. No noise was added to the calculated currents.

7.2.
Numerical results based on the linearization approach. We start by presenting numerical results obtained using the linearization approach presented in Section 4.3. We expect that this approach gives good results close to equilibrium, but becomes inaccurate for more general settings. To reconstruct D we used the boundary conditions u L = u R = 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 and for G the ones depicted in Table 1. In this example we set h = 0.01, D(u) = u(1 − u), G(u) = u and U = 0.001. This leads to the results shown in Figure 1. However, the quality of these results depends strongly on the value of U . This is reasonable, as this factor multiplies the error due to the linearization of the convection term. In fact, for larger values of U the error in the reconstruction increases strongly as shown in Figure  2. This leads to the conclusion that the linearization approach is of theoretical interest, as it demonstrates the feasibility of identifying both nonlinearities but its practical value is limited.  Table 2. Boundary conditions for the fully nonlinear problem convection. Consequently we expect a better reconstruction of either D or G in one or the other case. The exact nonlinearities are set to D(u) = u(1 − u) and G(u) = u 2 (1 − u). For the initial datum of D and G we choose the functions D * (u) = u(1 − u 2 ) and G * (u) = u 2 (1−u 2 ), which have the same values at u = {0, 1}. We choose a set of 13 different boundary conditions, given in Table 2. We reconstruct the nonlinearities for two different potentials, i.e. U = 0.3 and U = 1, see Figure 3. Here we observe that we get a good reconstruction for G(u) for U = 0.3, the reconstruction of the function D(u) is worse especially for smaller and larger values of u. Note that we only reconstruct the nonlinearities D and G on the interval [min(u R , u L ), max(u R , u L )], in the gray shaded areas the functions stay the same. For U = 1 the "convection" term dominates the equation, therefore we can reconstruct the nonlinear diffusion much better than for U = 0.3.

Example 2: Reconstruction for different sets of boundary values
The quality of the reconstruction depends heavily on the chosen boundary values. In the first example we chose a set of very different boundary conditions, e.g. large value on the right, small on the left or vice versa. If we choose very simple boundary conditions like u i (0) = u 0 + iδu and u i (1) = u 1 + iδu with i = 1, . . . m, the reconstructed solution are worse than in the first setting. All functions and parameters are the same as in the first example, except U = 0.5. Figure 4 shows the reconstructed solutions for the boundary conditions in Table 2 and (58) u i (0) = 0.15 + 0.05i and u i (1) = 0.3 + 0.05i with i = 0, . . . 13. We observe that the quality of the reconstruction depends heavily on the boundary conditions. For the first set of boundary conditions, i.e. Table 2, the nonlinearities can be reconstructed much better than for the second set, i.e.   Table 2 or (58).
on the chosen boundary conditions. In some applications it is possible to consider variations of the applied potential V , if this is not the case (e.g. in pedestrian motion where V = x and U determines the average walking speed) one should choose various boundary settings to obtain good results. Example 3: Reconstruction quality for different number of boundary conditions and basis functions Next we would like to understand the relationship between the number of basis functions ϕ l , the number of boundary conditions and the regularization parameters β = β 1 = β 2 . The nonlinearities are chosen the same as in the previous example, the strength of the potential is U = 0.5 and h = 0.005. We picked i = 5 sets of boundary conditions given by (59). Table 3 shows the L 2 norm of D(u) − u(1 − u) 2 + G(u) − u 2 (1 − u) 2 for different regularization parameters β and different mesh sizes h u in the discretization of the nonlinearities D and G. Table 4 shows the decrease in the error as the number of boundary measurements  Table 3. We observe that the discretization has strong regularization properties and that the error decreases for a larger set of boundary conditions. Example 2 and Table 4 show that the number of boundary conditions as well as their particular choice play a significant role in quality of the reconstruction. Example 4: Reconstruction with convexity constraints on the diffusivity D(u) In our final example we would like to include convexity constraints on the diffusivity D(u), since we posed this assumption throughout the paper. To include this constraint into our numerical simulations we add κ I max(D , 0) 2 dw to the minimization function (48). This additional term penalizes non-convex solutions, which can be observed in Figure 5. In the case of κ = 0.4, the diffusion D(u) looks like the "convex envelope" of the unregularized one, i.e. κ = 0. This first simple approach serves as a starting point for further research in the direction of parameter identification problems with convexity constraints.