BIOLOGICAL AND INDUSTRIAL MODELS MOTIVATING NONLOCAL CONSERVATION LAWS: A REVIEW OF ANALYTIC AND NUMERICAL RESULTS

This paper is devoted to the overview of recent results concerning nonlocal systems of conservation laws. First, we present a predator – prey model and, second, a model for the laser cutting of metals. In both cases, these equations lead to interesting pattern formation.


1.
Introduction.Several problems lead to the study of conservation (or balance) laws where the source and/or the flux are nonlocal functions of the solution.A first example is related to memory effects in elastodynamics, see [9], leading to balance laws with source terms depending on the past history.Then, in [13,28], the description of radiating gases led to source terms containing nonlocal operators in the space variable.In vehicular traffic, the realistic assumption that drivers react to what happens in a (possibly, forward) neighborhood of their position naturally leads to conservation laws with fluxes that are nonlocal in the space variable, see [15,7].Nonlocal interactions are even more relevant in crowd dynamics models, where a variety of nonlocal models were recently introduced, see for instance [2,12,17,21,34] and the review [6].Nonlocal models for the dynamics of granular materials are studied, e.g., in [3,14,25].Finally, we recall the wide research area related to structured populations, where standard models consist in initial -boundary value problems for balance laws with nonlocal boundary conditions, see for instance the book [33] or also the approach in [8], based on measure valued solutions to nonlocal balance laws.
In this paper, we present two models based on nonlocal balance laws.The first is a mixed hyperbolic -parabolic predator -prey model, introduced in [19].The latter is devoted to the dynamics of melted metal during the cutting of a metal plate by means of a laser beam, introduced in [16].In both cases, we recall the recently obtained well posedness results.Again in both cases, numerical integrations display qualitative features that are still lacking an analytic treatment.Indeed, both models lead to the formation of patterns, common also to some of the other aforementioned 50 RINALDO M. COLOMBO, FRANCESCA MARCELLINI AND ELENA ROSSI nonlocal models, see [12,17].In the case of the predator -prey system, these patterns are due to the reaching of a dynamic equilibrium between the hyperbolic drift of predators and the parabolic diffusion of prey, see Section 2. In the latter case of the laser model, this pattern provides a first description of the formation of ripples, see Section 3.
At present, a rigorous analytic definition of these patterns is missing, as well as detailed results on the qualitative behavior of solutions to balance laws with nonlocal fluxes.Below we exploit numerical integrations to investigate the role of specific parameters in the development of qualitative properties of solutions.Indeed, these parameters may well play the role of controls used to pursue specific goals, such as maximizing the predators' growth or minimizing the formation of ripples.In both cases of the predator-prey system and of the laser model, we single out two parameters.The first one, namely the prey natality γ in Section 2 and the laser speed v L in Section 3, enters the equations in a rather standard way.We expect that the L 1 -Lipschitz continuous dependence of solutions on these parameters can be obtained suitably refining the techniques in [16,19].
A less standard role, again in both cases, is played by the radius of the support of the convolution kernel, which is the second parameter we investigate.Indeed, nonlocal systems of conservation laws are known to lead to well posed Cauchy problems, see [17].As the convolution kernels converge to a Dirac delta, we obtain systems of conservation laws in several space dimensions, whose well posedness is a well known hard, and currently completely open, problem.

2.
A mixed parabolic -hyperbolic problem.In this section we consider the following mixed system to describe predator -prey dynamic introduced in [19]    Here, u = u(t, x) and w = w(t, x) are the predator and prey density respectively, evaluated at time t ∈ R + and position x ∈ R N .The coefficient α, β, γ, δ appearing in system (1) are all positive, µ is strictly positive.Thus, prey diffuse according to a parabolic equation, while predators evolve according to a nonlocal balance law.
In particular, the non locality is hidden in the function v , which is supposed to be a nonlocal and nonlinear function of the prey density w.Lotka -Volterra equations motivate the source terms: α represents the increase in the predator density due to feeding on prey, β is the mortality rate of predators, γ is the prey birth rate and δ their mortality rate due to predators.The model by Lotka [30] and Volterra [43], see also [32], are based on ordinary differential equations and implicitly assumes a homogeneous distribution in space of the two populations, while model (1) has the advantage of allowing for spacial variations, see also [32, § 1.2] and the references therein.More precisely, the flow u v(w) accounts for the direction preferred by predators, which then can move, for instance, towards the region where the concentration of prey is higher.A typical choice of the function v can be v(w) = κ grad(w * η) where κ > 0 is the maximal predator speed, while η is a positive smooth mollifier.The space convolution product (w(t) * η) (x) yields an average of the prey density at time t around position x.Observe that the radius of the support of η represents how far predators can "feel" the presence of prey, and, hence, the direction in which to hunt.The class of model ( 1) is studied in [19], under suitable assumptions on v.More precisely, existence, uniqueness, continuous dependence from the initial datum and various stability estimates for the solutions to (1) are proved.Solutions, understood in distributional sense, are constructed in the space L 1 ∩ L ∞ ∩ BV for predators and L 1 ∩ L ∞ for prey.Moreover, all analytic results hold in any space dimension.
An explicit numerical scheme for the discretisation of system (1) in two space dimensions is introduced in [36], where the convergence of the scheme is proved.Using the algorithm introduced in [36], qualitative properties of the solutions can be shown.In particular, we underline below a remarkable pattern formation, see Figure 1.
We recall the definition of solution to (1): We now recall what weak solution means in both cases.We fix t o , T ∈ R + , with T > t o , and denote I = [t o , T ].Concerning the parabolic problem, we give the following definition of solution, inspired by [35,Section 48.3].
and w(t o , x) = w o (x).
In a way similar to [22,Section 4.3] and [40,Section 3.5] we give the following definition of solution to the balance law.
Before stating the main analytic result, observe that system (1) is defined by a few real and positive parameters and by the function v. Hence, we require the following assumptions on this map: The bound on the L ∞ norm of v(w) by means of the L 1 norm of w, see the first and the third inequality above, is typical of a nonlocal operator, e.g.convolution.
Under reasonable regularity conditions on the kernel η, Lemma 4.1 in [19] ensures that (v) is satisfied by an operator v as in (2).
Introduce also the spaces Relying solely on the assumption (v), we state the main analytic result.
Theorem 2.4.[19, Theorem 2.2] Fix α, β, γ, δ ≥ 0 and µ > 0. Assume that v satisfies (v).Then, there exists a map R : R + × X + → X + with the following properties: in the sense of Definition 2.1.In particular, for all Local Lipschitz continuity in the initial datum: for all r > 0 and for all t ∈ R + , there exist a positive L(t, r) such that for all (u 1 , w 1 ), (u 2 , w 2 ) ∈ X + with for i = 1, 2, the following estimate holds: 4. Growth estimates: for all (u o , w o ) ∈ X + and for all t ∈ R + , denote An explicit estimate of the Lipschitz constant L(t, r) is provided in [19,Equation (4.37)].The growth estimates in point 4 provide an a priori control on the total mass and on the peaks of the solution.Point 5 means that the solution to the balance law propagates with finite speed.
For the proof of Theorem 2.4, see [19,Paragraph 4.3].There, the result is proved exploiting a fixed point argument and careful estimates on the parabolic problem and, separately, on the balance law (7) with suitable assumptions on the functions a, b and c.Remark 1. Concerning the parabolic problem (6), observe that we are interested in L 1 solutions defined on the whole space, and not L 2 on a subset of R N as it is conventional in literature, see [4,31,35].Detailed proofs and estimates for the parabolic problem are provided in [19], inspired by the results in [4,23].
Remark 2. It is interesting to note that, even if both Cauchy Problems ( 6) and ( 7) are non autonomous problems and each of them generate a process, see [19, Proposition 2.5 and Proposition 2.8], the mixed system (1) generates a semigroup.Remark 3. Consider the following generalization of (1): with the assumptions f ∈ C 2 (R; R) and f (0) = 0. Existence, uniqueness and stability estimates follow by a straightforward extension of [19, Proposition 2.8] and Theorem 2.4.
2.1.Numerical integrations.We exploit the numerical algorithm introduced in [36] for the discretisation of system (8), generalization of (1), in two space dimensions in order to provide some numerical integrations of the model proposed and to show some qualitative properties of the solutions.
The approximation scheme consists of a finite difference scheme for the parabolic part and a Lax -Friedrichs type finite volume method for the balance law.The source terms of both equations are treated using a second order Runge -Kutta method, see [36,Section 2] for additional details.The classical operator splitting technique is to combine the differential parts to the source terms.
In this numerical integration, the boundary data are assigned as in [36, Paragraph 5.1].
We focus on integrating system (1) in the two-dimensional case and we use the vector field v in (2), where the kernel function η is chosen as follows In particular, we present two different situations in which an interesting asymptotic state arises.In both cases, we set v as in ( 2) and η as in (9).
The solution is computed up to time T max = 3 on a mesh of width ∆x = ∆y = 0.0025 and the result is displayed in Figure 1.
At the beginning, predators almost disappear, while moving towards the central region of the numerical domain where prey grow most, due to the chosen boundary conditions.Then, they start to increase and slowly a regular and quite robust pattern appears: predators focus in small regions, regularly distributed, surrounded by prey.
The analytic explanation of this pattern could be the following: there is a sort of dynamic equilibrium between the first order nonlocal differential operator in the predator equation and the Laplacian in the prey equation.In the zones where predators are focused, their feeding on prey causes "holes" in the prey distribution.As a consequence, due to symmetry considerations, the average gradient of the prey density almost vanishes, and hence predators almost do not move.At the same time, these "holes" are continuously filled by the reproduction and the diffusion of the prey, which keep providing new nutrient to predators.
Coherently with this explanation, numerical integrations confirm that this asymptotic pattern essentially depends on the size of the support of η.Indeed, we Figure 1.Numerical integration of (1)-( 9)-( 10) with initial datum (11).Top row: u; middle row: w.Darker colors indicate higher densities.First, predators decrease and move towards the central region.Then, a discrete regular pattern arises: predators focus in small regions regularly distributed.This solution was obtained with space mesh size ∆x = ∆y = 0.0025.On the bottom row, to enhance the visibility of the prey distribution, we plot w− w, where w is the average of w over the computational domain, and we let the color levels vary, i.e. the same colors may mean different density in different pictures.
We can see that as decreases, also the distance among pairwise nearest peaks in the asymptotic predator density u decreases, and more peaks are possible, see Figure 2.
An other interesting feature is to see what happens when one of the parameters in (10) varies, keeping the initial datum (11) fixed.For instance, let γ assume the following values: while all the other parameters are the same.The parameter γ represents the birth rate of prey w: as γ gets smaller, less prey are born.The same pattern as before appears for every values of γ, with difference in number, distribution and heights Only the solution for predators u is displayed, since the structure is more visible.Darker colors indicate higher densities.As decreases, the distance among peaks in u decreases and more peaks are possible.The space mesh size is ∆x = ∆y = 0.0025.
of the peaks, see Figure 3, where we focus only on the solution for predators u at a particular time t = 2.52.In Figure 4 we can see the evolution in time of the L 1 -norm of the densities u and w.For larger γ, the total mass of predators u grows more, see the graphs on the left of Figure 4.The same seems to happen to the total mass of prey.However, as time passes, the great amounts of predators present for larger values of γ permits to control the prey, leading to a reduction of their mass.1)-( 9)-( 10) with initial datum (11) for different values of γ, i.e. γ = 0.1, γ = 0.2, γ = 0.4, γ = 0.6 and γ = 0.8.These integrals represent the total mass of predators and prey respectively.
The solution is computed up to time T max = 3 on a mesh of width ∆x = ∆y = 0.0025 and the result is displayed in Figure 5. Figure 6 shows the evolution of the total mass, i.e. the L 1 -norm in space, of predators u and prey w over time.As in the previous example, we can see the formation of a discrete pattern.Here, predators u are in the centre of the numerical domain from the beginning, while prey w surround them, forming a ring around predators.Prey start to diffuse, filling the centre of this ring.However, their density is still greater in the original ring than in its centre, so predators move away from the centre of the domain.Prey keep on diffusing, and predators start to focus in small regions.Peaks of predator density are initially distributed along two columns, placed rather in the centre of the domain.As prey diffuse and their density increases not only in the centre of the numerical domain, see Figure 6, some predators move away form the centre, and the peaks in their density at the end are distributed along six columns, on an area wider than the initial one.1)-( 9)-( 12) with initial datum (13).Only u is displayed, at different times, the contours of w being "complementary", as in example 1, see Section 2.1.1.Darker colors indicate higher densities.Note that, to enhance the visibility of the pattern, the same color in different figures may indicate a different density.This solution was obtained with space mesh size ∆x = ∆y = 0.0025.Figure 6.The graphs display the integral of u, on the left, and w, on the right, representing the total mass of predators and prey respectively.They refer to the solution to (1)-( 9)-( 12) with initial datum (13).
2.2.Future research directions.The mixed system (1) suggests several other research directions.First of all, inspired by the outcome of the numerical integrations, we expect it is possible to prove the Lipschitz dependence of the solution to (1) from the parameters therein, i.e. α, β, γ and δ.
Secondly, it would be worth studying system (1) with different source terms, for instance considering the carrying capacity of the species.
Lastly, up to now the problem is studied on all R N .It would be interesting to study the mixed system (1) in a bounded domain.Concerning the parabolic case, this would not be a problem, since in the literature many results on parabolic equations in bounded domain can be found.As for the hyperbolic equation, it is necessary first to settle a few issues.A key reference is the classical paper by Bardos, Leroux and Nédélec [5], but more estimates are needed in anticipation of the coupling with the parabolic equation, see also [20].3. A model for laser cutting.In this section we present a new model for the cutting of metal plates by means of a laser beam, see [16].
Currently, two types of laser are used in this kind of technology: CO 2 lasers and fiber lasers.The former ones are more powerful and accurate, but also more expensive.The latter ones are much less expensive but, unfortunately, they have the drawback of the formation of a sort of ripples along the sides of the cut, see Figure 8, right.Typical tasks in laser cutting applications involve finding those cutting parameters, e.g.laser power or cutting speed, such that these ripples are minimal.
Through a phenomenological model based on a nonlocal system of balance laws in two space dimensions, we aim to the description of the formation of these ripples, whose insurgence deeply affects the quality of the cuts.
From an analytic point of view, the formulation of this model leads to a Cauchy problem whose well posedness was not yet proved.Adapting techniques taken from the theory of conservation laws is obtained a new theorem for existence, uniqueness and continuous dependence of the solutions from the initial data for a class of nonlinear systems which includes the present model.
For the evaluation of qualitative properties of solutions of the proposed model, the use of numerical integration is indispensable.Here, by using realistic numeric parameters taken from the specialized engineering literature, we obtain the formation of a geometry similar to the ripples observed in industrial cuts.We remark that in this construction neither the initial data nor the parameters in the equations contain any oscillating term.Indeed, in perspective, this model could allow to preliminary test the parameters of cutting, to optimize and to minimize the presence of ripples in real cuts.
For its industrial relevance, laser technology is widely considered in the specialized literature.In particular, for research results related to the modeling, model analysis and numerical simulation of laser cutting, we refer to [24,26,37,38,39,42,44].The paper in [10] is devoted to the problem of ripples with a framework different from the present one.A 1D system of balance laws is used to describe this phenomenon, but is difficult a treatment of the problem from the analytic point of view.These difficulties led to the formulation of the model in [16].
We describe this process as follow.We consider a 3D geometric framework, see Figure 7 left, where the plate lies on the z = 0 plain and the laser beam is parallel to the z axis.We distinguish the height h s of the solid metal and that of the melted part, denoted h m ; see Figure 7, right.We base our model on the following 2D system of balance laws x 1 Left, the description of the 3D geometric framework: the laser beam is parallel to the z axis, while the plate lies on the z = 0 plain.Right, the distinction between the melted part h m and the solid one h s .
The source term L is directly related to the intensity of the laser: it describes the net rate at which the solid part turns into melted.It is given by where i = i(t, x) is the laser intensity.The denominator is the squared cosine of an averaged incidence angle of the laser on the surface z = H(t, x), where H = h s +h m , see Figure 8.The function η is a smooth function defined in R 2 , so that Indeed, the term grad x (η * H(t)) (x) is the average gradient at position x and time t of the surface z = H(t, x).In ( 14), the vector V = V (t, x) is given by where w = w(t, x) is the wind speed, provoked around the beam, which pushes the melted material downwards.The coefficient τ g in ( 16) is related to the shear stress, inspired by [10,44].The denominator in ( 16) is due to a (smooth) normalization of the direction − grad x (η * h s ) of the average steepest descent along the surface z = h s (t, x).
For the the laser intensity function i = i(t, x) and for the wind function w = w(t, x) we consider where x = x L (t) is the center of the laser beam and both the functions I and W can be reasonably described through a compactly supported bell shaped function centered at the location of the moving focus of the laser beam.In general, the radius of the surface where the wind blows downward is a few times that of the laser beam.Therefore, the model proposed is the following: From an analytic point of view, the formulation of this model leads to consider a wider class of Cauchy Problems for systems of nonlocal balance laws in several space dimensions.In fact, in [16], a new theorem is obtained ensuring existence, uniqueness and continuous dependence for a class of nonlinear systems which includes the model in (18).
More precisely, the following family of problems is considered: Above, t ∈ [0, +∞[ is time, x ∈ R N is the space coordinate and the vector u ≡ (u 1 , . . ., u n ), with u i = u i (t, x), is the unknown.The term φ ≡ (φ 1 , . . ., φ n ), with φ i (t, x, u i , A) ∈ R N , is the flow and Φ ≡ (Φ 1 , . . ., Φ n ), with Φ i (t, x, u i , A) ∈ R, is the source term.The function θ is a smooth function defined in R N attaining as values m × n matrices, so that The above system (19) enjoys the property that the equations are coupled only through the nonlocal convolution term θ * u.This fact plays a key role in proving well posedness.Indeed, for small times, system (19) admits a unique solution u = u(t, x).Moreover, u is proved to be a continuous function of time with respect to the L 1 topology and an L 1 -Lipschitz continuous function of the initial datum ū.We recall first the definition of solution to system (19).
3.1.Numerical integration.To display qualitative properties of the solutions of the proposed model ( 18), the use of numerical integrations is indispensable.Here, using realistic numeric parameters, we obtain the formation of a geometry similar to that of ripples observed in real industrial cuts.We remark that, in the present construction, neither the initial data nor the parameters in the equations contain any oscillating term.Nonetheless, we obtain an oscillating profile, rather similar to what happens in real cuts; see for instance Figure 8, right, refer also to [16,Paragraph 3.1].We use below the numerical method presented in [1], which is specifically devoted to systems of nonlocal conservation laws.More precisely, we use a Lax -Friedrichs type algorithm for the convective part and a first order explicit forward Euler method for the ordinary differential equations arising from the source terms; see also [29,Section 12.1].We show here a numerical integration, where the laser beam describes 3 segments.More precisely, the integration is computed for t ∈ [0, 1.19], time being measured in seconds, the mesh size is 8 • 10 −5 and the domain is the rectangle [0, 20] × [−5, 4.12].All lengths are measured in millimeters.The laser beam trajectory is given by corresponding to a laser beam moving at v L = 40 mm sec .The initial hole is drilled centered at (3.6, 0), in the interior of the metal plate.At each corner, the beam stops for 0.1sec.We used the data in [41], see also [10,Table 1].In particular, we set τ g = 4 and the convolution kernel is given by η As initial datum, we choose hm (x) = 0 and hs (x) = 4.5 for all x ∈ R 2 , which represents a metal plate 4.5 mm thick.Finally, the wind and laser functions are given by (17) with corresponding to a laser beam with radius 1.2 mm, see [41].The resulting integration is in Figure 9, where we can see these sort of ripples, generated with neither parameters nor initial data that contain anything oscillating.
It is of interest to investigate the dependence of the solution from the various parameters.While a rigorous analytic study is left to a forthcoming research, here we provide an insight on the basis of numerical integrations.18)-( 22) with the data and the realistic numeric parameters listed in § 3.1, see [10,41].Displayed are the contour plots of the solid metal h s , representing the shape of the metal remaining after the cut.
A parameter of key relevance is the horizontal speed at which the vertical laser beam moves.As Figure 10 shows, there is clear reduction in ripples as this speed Figure 10.Numerical integration of ( 18)- (22) with the data and the realistic numeric parameters listed in § 3.1, see [10,41], with different laser speeds v L : from left to right, 30 mm sec , 40 mm sec and 80 mm sec .A higher speed apparently reduces ripples.
increases.Note also that the faster beam produces a narrower cut, as is to be expected.
A more subtle parameter is the radius r η of the support of the convolution kernel η.The proof of the L 1 -Lipschitz continuous dependence of the solution from r η does not appear to be easily at reach.Nevertheless, the integrations in Figure 11 display a convincing stability.
3.2.Future research directions.The present model rises various questions.First, from the modeling point of view, a reasonable formal derivation of the convective part and of the source term in (18) would be certainly interesting.
In view of specific applications of this model, the dependence of solutions from parameters is also of interest and the techniques in [18] are likely to provide the necessary starting point.However, the inverse problem, i.e., determining the values Finally, we recall the obvious relevance of the control problem for (18).In view of its industrial interest, we expect that the ability of tuning the various parameters to optimize the quality of cuts and, more generally, the efficiency of the whole process, is of obvious interest.

Figure 2 .
Figure 2. Numerical solution of (1)-(9)-(10) with initial datum (11) computed at time t = 2.52 for different values of , i.e. from the left = 0.03125, = 0.0625, = 0.125 and = 0.25.Only the solution for predators u is displayed, since the structure is more visible.Darker colors indicate higher densities.As decreases, the distance among peaks in u decreases and more peaks are possible.The space mesh size is ∆x = ∆y = 0.0025.

Figure 5 .
Figure 5. Numerical integration of (1)-(9)-(12) with initial datum(13).Only u is displayed, at different times, the contours of w being "complementary", as in example 1, see Section 2.1.1.Darker colors indicate higher densities.Note that, to enhance the visibility of the pattern, the same color in different figures may indicate a different density.This solution was obtained with space mesh size ∆x = ∆y = 0.0025.

x 1 zzFigure 8 .
Figure 8. Left, the incidence angle α of the laser beam on the surface z = H(t, x).Right, image of a cutting surface with ripples.

Figure 9 .
Figure 9. Numerical integration of (18)-(22) with the data and the realistic numeric parameters listed in § 3.1, see[10,41].Displayed are the contour plots of the solid metal h s , representing the shape of the metal remaining after the cut.