WEAK SOLUTIONS FOR A SIMPLE HYPERBOLIC SYSTEM

The model studied concerns a simple ﬁrst-order hyperbolic system. The solutions in which one is most interested have discontinuities which persist for all time, and therefore need to be interpreted as weak solutions. We demonstrate existence and uniqueness for such weak solutions, identifying a canonical ‘ exact ’ solution which is everywhere deﬁned. The direct method used is guided by the theory of measure-valued diﬀusions. The method is more eﬀective than the method of characteristics, and has the advantage that it leads immediately to the McKean representation without recourse to Itˆo’s formula. We then conduct computer studies of our model, both by integration schemes (which do use characteristics) and by ‘random simulation’.


Introduction
After a change of the original notation with u replacing 1 − u, the FKPP equation (Fisher [6], Kolmogorov, Petrovskii & Piskunov [9]) reads: McKean [12,13] showed that one can prove existence and uniqueness results for certain [0, 1]valued solutions by using martingales to describe u as an explicit functional of a certain branching Brownian motion. He was thereby able to obtain results on convergence to travelling waves for suitable initial data f . Neveu [15] gives an important completion of McKean's treatment.
Two recent papers -Champneys, Harris, Toland, Warren & Williams [2] and Lyne [10] have each presented both analytic and probabilistic studies of simple extensions of the FKPP equation. These papers used the probabilist's golden rule that Itô's formula leads to martingales (see, for example, Rogers & Williams [17]). But the uses of Itô's formula involved the 'formal generator' of the branching process in a way which might cause some unease to analysts. The model studied by Lyne concerns a simple first-order hyperbolic system (of similar nature to those studied in Dunbar [4], Othmer, Dunbar & Alt [16], Holmes [8] and Hadeler [7] -see the discussion in Lyne [11] for the relations between the models), a generalization of which we consider here. The solutions in which one is most interested have discontinuities which persist for all time, and therefore need to be interpreted as weak solutions. Section 2 settles existence and uniqueness for such weak solutions, identifying a canonical 'exact' solution which is everywhere defined. The direct method used is guided by the theory of measure-valued diffusions, MVDs, (see, for example, Dawson [3] and Dynkin [5]). (We stress that no knowledge of MVD theory is assumed here.) The method is more effective than the method of characteristics, and has the advantage that it leads immediately to the McKean representation without recourse to Itô's formula.
Having (we hope) satisfied the requirements of analysis in Section 2, we turn in Section 3with more freedom -to the question of computer studies of our model, both by integration schemes (which do use characteristics) and by 'random simulation'. Numerical analysts would wish for more rigour in Section 3. We settle for a certain amount of cross-checking.
Notational point.W enever use u t to denote ∂u ∂t ,r a t h e ru t denotes u at time t.
2 Theoretical results

Our hyperbolic system
Let I be a finite set (with the discrete topology); and let B and R be functions of I with R ≥ 0. Let Q be an I × I matrix with non-negative off-diagonal elements and zero row sums.
Let f be a Borel function on R × I with 0 ≤ f ≤ 1. Let u, written (t, x, j) → u t (x, j)a n d regarded as a column vector in j when multiplied by Q, be a Borel function on [0, ∞) × R × I with 0 ≤ u ≤ 1. Suppose that u is a weak solution of In full, the first equation reads: By the statement that u is a weak solution, we mean that for t>0 and a test function ϕ ∈ C 1,1,0 and in space) thought of as a row vector in j,

Analytic statement of some results
Shortly we shall shortly reformulate these results probabilistically.
Introduce the unique one- , where θ ∈ R and g is a function (or column vector) on I,t h e n In regard to the existence of {P −R t : t ≥ 0}, see the discussion around equation (5) below. By the Riesz representation theorem, {P −R t : t ≥ 0} has a canonical extension to a semigroup on B b (R × I), the space of bounded Borel functions on R × I.

Equation
(2) implies that for each t>0, u satisfies: for almost every x. This is proved as follows. Standard Fourier theory shows that for ψ 0 (·, ·)i n C ∞ K (R × I), Now take ϕ s (x, j)=ψ t−s (x, j) in equation (2).
Then it is almost immediate by the usual Picard/Gronwall argument that u * := lim u (n) exists monotonically and uniformly on each [0,t]×R ×I,a n ds og i v e sthe exact solution of equation (4): it is a solution in which there are no 'exceptional sets'; and if v is another exact solution to equation (4) with 0 ≤ v ≤ 1, then v is equal to u * everywhere.M o r e o v e r ,u * is a weak solution of equation (1).

Probabilistic interpretations/proofs
Let {η t : t ≥ 0} be a Markov chain on I with Q-matrix Q. Define and set where E x,j is the measure corresponding to starting position (ξ 0 ,η 0 )=(x, j). Let us check that this agrees with the semigroup of equation (3). Fix θ ∈ R. For a vector g on I, define If we write . = to signify equality modulo differentials of local martingales, then (see Rogers & Williams [17]) d{g(η t )} . =(Qg)(η t )dt, In fact, the difference between the integrals of the two sides of equation (6) is bounded on each [0,t], and so is not just a local, but a true, martingale. Hence where S t is the linear map on vectors on I defined by Since S 0 is the identity, equation (3) now follows.
Guided by McKean, we now construct a branching Markov process related to the hyperbolic equation (1). Time 0 sees the birth of one particle, labelled 1, which has 'type' Y 1 (0) in I and 'position' X 1 (0) in R.A t t i m e t ≥ 0, there are N (t) particles which, when labelled in order of birth, have 'types' Y 1 (t),...,Y N (t) (t)i nI,a n dp o s i t i o n sX 1 (t),...,X N (t) (t)i nR.T h et y p e of each particle behaves (independently of previous history, of the behaviour of other particles currently alive, etc) as a Markov chain on I with Q-matrix Q. A particle of type j moves on R with constant speed B(j), and gives birth to a new particle of its own type with rate R(j), so that in small time h, independently of 'everything else', it gives birth with probability R(j)h +o(h). Particles live forever, once born. We write P x,j and E x,j for the probability and expectation corresponding to the situation when X 1 (0) = x and Y 1 (0) = j.
With f as our 'initial value for u 0 ', let We now utilise an obvious argument. Let T be the time of the first birth after time 0, so that T is the birth-time of particle 2. Because of the rôle of R as birth-rate function, we have From time T on, the family tree of particle 2 evolves independently of its complement in the family tree of particle 1. We therefore have Hence, for s ≤ t, Since we see that v satisfies equation (4) exactly. Hence v = u * everywhere, and we have the McKean representation

Convergence to travelling waves: the easy case
Suppose that

Condition ( †)
We consider the very special situation in which Q is irreducible and there is a state j 0 ∈ I with We will refer to these conditions as ( †).
When ( †) holds there will almost surely exist at some time a particle of type j 0 which has an infinite 'line of descent' consisting entirely of particles of type j 0 . Thus there will be a random interval [σ, ∞), which we choose to be maximal, such that for some random constant A, Then and w(x + tB(j 0 ),j) is a travelling-wave solution of equation (1).
In this case, For numerical and simulation studies of such a case, see section 3 below.

The difficult cases
It is hoped to make the difficult cases when ( †) fails to hold the subject of another paper giving direct proofs for this simple situation of results Indeed, a(t) may behave like a multiple of log t or of log log t. Such results follow from deep results in existing literature. The classic paper on the 'logarithmic correction'f o rt h eF i s h e r equation is that of Bramson [1].

Numerical analysis
We use both probabilistic simulation of the branching process and finite difference methods (the most effective being upwinding along characteristics) to study the initial value problem for the case where I = {1, 2}, that is, we have a pair of coupled first order PDEs (as studied in Lyne [10]). As in Dunbar [4] and Holmes [8], these can be rewritten as a single second order PDE which is a generalization of the telegraph equation. It is convenient to change to moving coordinates (moving at a speed of 1 2 (B(1) + B(2))) and then re-scale time so that the coefficients of ∂u ∂x are 1 and −1. This is possible unless B(1) = B(2) -this case reduces to a pair of first order ODEs and is dealt with by Lyne [10]. For the remainder of this paper, we shall denote u t (x, 1) by u(t, x)a n du t (x, 2) by v(t, x)a n ds e tB(1) = 1,B(2) = −1. Particularly interesting is Heaviside initial data, that is, We present first the skeleton of the C program used to produce the numerical solution plotted in Figures 1 and 4. It implements a naive Euler method along the characteristics of the system, and a modification of the Euler method. The figures were produced using the modified method, but output of the two methods is practically indistinguishable.
/* EULER METHODS This is only part of a program.
Use of naive Euler methods for a simple hyperbolic system du/dt = du/dx + f(u,v), f(u,v) = q1 (v-u) + r1 u(u-1); dv/dt = -dv/dx + g(u,v), g(u,v) = q2 (u-v) + r2 v(v-1); with Heaviside initial data u(0,x) = v(0,x) = 1 for x > 0, 0f o rx< =0 . Several finite difference methods were implemented on a rectangular lattice. These all proved to be less effective than the Euler method used along the characteristics (via the customized lattice, that is, using the characteristics to build the grid). The most effective of these schemes was the Lax-Wendroff scheme, as implemented in the following program. Mitchell and Griffiths [14,Chapter 4] give a good discussion of the Lax-Wendroff scheme and hyperbolic equations in general; see also Strikwerda [18].
A solution plotted from this program is presented in Figure 3. It is fairly similar to the plots in Figures 1 and 2 for the same parameter values using the other methods investigated, but it suffers from typical Gibbs phenomena, and does not maintain the sharp discontinuities actually present in the true solution along each characteristic. However, for the parameter values given, these discontinuities decay exponentially to zero, and for longer times the solutions of all three methods agree very well. In cases where a discontinuity does not decay to zero (but instead to a finite size between 0 and 1, as in Figures 4 and 5), the Lax-Wendroff method is visibly worse because the discontinuity is smeared out over several grid points. Away from the discontinuity agreement is good. /* LAX-WENDROFF SCHEME This is only part of a program.
nu denotes the next u-array, that is, u one time step later.

Probabilistic simulation
The key to the probabilistic simulation is a simple recursive function called Life. This function tracks the path of an individual particle, updating the record of the left-most position yet reached by any particle each time the record is broken, and storing in arrays the time, place and type of each birth. Then, upon completion of a particle's run (since we only run each particle up until a pre-specified maximum time), we check to see if we have any more particles to do, and run the Life function on them.
Each particle is dealt with in a series of segments (the function DoSegment). An exponential random variable is generated (by Rexp)-the length of time until either the particle splits into two or the particle changes type. The position of the particle is updated by simply adding the length of time multiplied by its speed in its type for the segment to the current position. If the length of the segment takes us past our maximum time, we have completed the life story for that particle, and move on to the next one (incrementing c, the number of our current particle). Then we determine which event it was that did actually occur -birth or mutation. For birth we call the function Create, which stores the time, position and type in the arrays tt, xx and yy respectively. We also increment the counter n to inform us there is one more particle to be dealt with later. We then do another segment. If we change type, then we flip the type variable y and do the next segment. Once we reach a point where we have completed the life story of a particle and there are no more sets of birth information unused (i.e. c greater than n), we have finished the simulation run.
If one is interested in calculating solely, for example, the left-most particle position, much calculation can be saved. Since a particle can only travel at speeds 1 and −1 we have immediate bounds on its future position (and identical bounds on all its future descendants). Thus, if the maximum time we are running until is T and the current record for left-most position is X,w e can discard any particle in the simulation (denoting its position and time by (t, x)) for which x − (T − t) >X; particles satisfying this inequality have no chance of changing the record. The fact that their descendants also cannot break the record means we can discard the parent and save even generating the descendants.

/* SIMULATION OF BRANCHING PROCESS MODEL
This is only part of a program. It shows how to extract information on the left-most particle from the simulation (into the arrays leftu and leftv), and save some computation if this is all we are interested in (by pruning away particles far away from the left-most). While in type 1 particle moves at speed b [1] (wlog set to 1), mutates at rate q [1] and breeds at rate r[1] -in type 2 particle moves at speed b [2] (wlog set to -1), mutates at rate q [2] and breeds at rate r[2] -q and r user inputs.
The program notes the position of the left-most particle in each of NUMRUN (1000) simulations starting from 1 particle at the origin (doing 1000 runs for that particle being type 0, and 1000 for type 1), observing each simulation at T (6) points, each GAP (user input) units apart.
*/ #define NUMRUN 1000 #define MAXPART 1000000 /* limit on particle numbers, warning if exceeded */ #define T 6 double b [3] Figures 2 and 5. These agree very well with those produced by the Euler method for corresponding parameters in Figures 1 and  4. When more simulation is done the probabilistic plots are smoother and agreement is even better. This simulation method could of course easily be extended to the n-type case.

Discussion of figures
Consider the situation in which there is just one particle at time 0, of type 2 and with position x. We know that v(t, x) is the probability that all particles are to the right of 0 at time t.S i n c e no particle can travel left at speed greater than 1, v(t, x)=1forx>t. However, if x = t,t h e n v(t, t) is the probability that up to time t our initial particle has no line of descent consisting only of particles of type 2: in other words, that every descendant of our initial particle spends some time before t moving right (in which case it can never get to 0 at time t). It is therefore clear that there is a positive jump 1 − v(t, t)i nv(t, ·)a tt i m et, and that this jump may be calculated by regarding any particle of type 1 as 'dead' and ignoring it and its descendants. Precisely, the jump 1 − v(t, t) is the probability that a continuous-time branching process starting from 1 particle, and with birth-rate r 2 and death-rate q 2 (per individual) survives until time t.I f q 2 >r 2 , the situation in Figures 1, 2 and 3, then this survival probability tends exponentially fast to 0. If r 2 >q 2 , the situation in Figures 4 and 5, then v(t, t) → 1 − π, where the extinction probability π satisfies π = q 2 q 2 + r 2 + r 2 q 2 + r 2 π 2 , so that π = q 2 /r 2 .
Two other features of the pictures are worthy of comment. Firstly, the fast convergence to the travelling wave in Figures 4 and 5 illustrates the case discussed in section 2.4. Secondly, the almost linear nature of v(t, ·) for small t, clearly apparent in Figures 1-3   q 1 =1.00,q 2 =2.00,r 1 =2.00,r 2 =1.00.
Distribution of left-most particle for t =0.0(0.400)2.400 1000 Runs each from type 1, and from type 2, initial particle   follows. Again, consider the situation in which there is just one particle at time 0, of type 2 and with position x. For small t, the dominant contribution to v(t, x) will arise from cases where there is a random time S before t at which the particle changes type. Thus the particle moves left for a time S and right for a time t − S, ending up at x − S + t − S.T h u s ,f o rs m a l lt, v(t, x) ≈ P(x − S + t − S>0) = P S< 1 2 (x + t) =1− e − 1 2 q 2 (x+t) ≈ 1 2 q 2 (x + t).