New travelling wave solutions of the Porous-Fisher model with a moving boundary

. We examine travelling wave solutions of the Porous-Fisher model, ∂ t u ( x, t ) = u ( x, t ) [1 − u ( x, t )] + ∂ x [ u ( x, t ) ∂ x u ( x, t )], with a Stefan-like condition at the moving front, x = L ( t ). Travelling wave solutions of this model have several novel characteristics. These travelling wave solutions: (i) move with a speed that is slower than the more standard Porous-Fisher model, c < 1 / √ 2; (ii) never lead to population extinction; (iii) have compact support and a well-deﬁned moving front, and (iv) the travelling wave proﬁles have an inﬁnite slope at the moving front. Using asymptotic analysis in two distinct parameter regimes, c → 0 + and c → 1 / √ 2 − , we obtain closed-form mathematical expressions for the travelling wave shape and speed. These approximations compare well with numerical solutions of the full problem.


Introduction
Travelling waves arise in many fields, including ecology [1][2][3][4], cell biology [5][6][7][8][9][10][11], and industrial applications involving heat and mass transfer [12][13][14][15]. Such processes are often modelled using reaction-diffusion equations and, depending on the choice of reaction and diffusion terms, three broad classes of monotone travelling waves are commonly reported ( Figure 1). The most commonly-observed travelling wave is a smooth front (Figure 1a), whereby the concentration, u(x, t), is a monotone decreasing function decaying to zero as x → ∞. For example, travelling wave solutions of the Fisher-KPP model [1,2,7], are smooth. Unfortunately, smooth fronts do not have compact support, which makes defining the "edge" of the moving front ambiguous [8,9,16]. This feature of the Fisher-KPP model means that it can be hard to apply to practical problems, such as cell invasion [7][8][9], where well-defined fronts are often observed.

Fisher-Stefan Fisher-KPP
Porous-Fisher-Stefan To obtain travelling wave solutiuons with a well-defined front, two main modifications of the Fisher-KPP model have been proposed. The first modification involves incorporating nonlinear degenerate diffusion, whereby the concentration flux is generalized to −D(u)∂ x u, with D(0) = 0. One common choice of degenerate diffusivity is D(u) ≡ u, which leads to the Porous-Fisher model [10,11,[17][18][19][20][21][22][23]: The Porous-Fisher model supports travelling wave solutions that move with speed c ≥ c min , where c min = 1/ √ 2 [1,19,20,24]. Phase plane analysis shows that travelling wave solutions of the Porous-Fisher model are sharp-fronted, with compact support, when c = 1/ √ 2 [1,20,24]. In contrast, travelling wave solutions of the Porous-Fisher model with c > 1/ √ 2 are smooth and do not have compact support [1,24]. Interestingly, travelling wave solutions of the Porous-Fisher model with c < 1/ √ 2 have never been reported.
The second modification of the Fisher-KPP model that leads to a well-defined front is to maintain the use of linear diffusion, but to incorporate a moving boundary condition so that we consider the Fisher-KPP model on −∞ < x < L(t) [4,7,25]. This second modification of the Fisher-KPP model has been called the Fisher-Stefan model [7,25], which incorporates a Stefan-like condition at the moving boundary, L(t), to relate the concentration flux, −∂ x u(x, t), with the speed of the moving boundary. One of the limitations of the Fisher-Stefan model is that this model allows populations to become extinct [7,25], since there is an outward flux of at the leading edge, x = L(t). While the Fisher-Stefan model has the advantage that it can lead to sharp-fronted travelling waves, the physical or biological explanation of the outward flux at x = L(t) is not obvious.
Both the Porous-Fisher and the Fisher-Stefan models have the advantage that they lead to travelling wave solutions with a well-defined front (Figure 1b,c). While these two modifications of the Fisher-KPP model have been considered previously, the extension of combining nonlinear degenerate diffusion with a moving boundary condition has yet to be considered. We refer to this combination of the Porous-Fisher model with a moving boundary condition as the Porous-Fisher-Stefan model, which we define as the Porous-Fisher equation (2) on x ∈ (−∞, L(t)] with a Stefan-like condition stating that the speed of the moving front, dL(t)/dt, is proportional to the nonlinear concentration flux, −u(x, t)∂ x u(x, t), with constant of proportionality κ > 0. Unlike the Fisher-Stefan model [7,25], the incorporation of nonlinear degenerate diffusion in the Porous-Fisher-Stefan model prevents u(x, t) from being driven to extinction. Interestingly, preliminary numerical solutions of the Porous-Fisher-Stefan model suggest that travelling wave solutions exist with 0 ≤ c < 1/ √ 2 and that these sharp-fronted travelling waves have infinite slope at x = L(t) (Figure 1d). Neither of these properties have been reported or analyzed previously.
In this work, we examine sharp-fronted travelling waves that arise in the Porous-Fisher-Stefan model. By transforming this model into travelling wave coordinates, we examine the phase plane for various regimes of the wave speed c. Using asymptotic analysis, we obtain approximations of the resulting phase plane trajectories when c 1 and when c is close to the critical wave speed 1/ √ 2. Additionally, we determine the relationship between the travelling wave speed c and the Stefan parameter κ. In doing so, we determine an approximate form of the travelling wave front that matches numericallycomputed travelling wave solutions with high accuracy.

Travelling waves in the Porous-Fisher-Stefan model
We consider the non-dimensional Porous-Fisher model, describing the concentration u(x, t) ∈ [0, 1] with a Stefan-like condition at the moving boundary x = L(t): The Stefan-like condition relates the speed of the moving front, dL(t)/dt, to the concentration flux, −u(x, t)∂ x u(x, t), via the constant κ > 0. Preliminary numerical solutions of (3)-(5) suggest that travelling waves exist and that these waves move with speed 0 ≤ c < 1/ √ 2 with infinite gradient at x = L(t). Consequently, we are motivated to examine travelling wave solutions of the Porous-Fisher-Stefan model and to determine the relationship between c and κ. To do this, we define φ(x, t) = [u(x, t)] 2 to obtain a corresponding PDE with a linear Stefan-like condition at x = L(t): To study travelling wave solutions of (6)-(8), we transform the system into travelling wave coordinates via z = x − L 0 − ct, where z ∈ (−∞, 0]. Noting that when x = L(t), we have L(t) = L 0 + ct, and hence, dL(t)/dt = c. This change of coordinates gives where = d/dz. It is convenient to write (9) as a system of two first order differential equations: A general explicit solution of (11)- (12) is not obvious, so we seek to study solutions of this system using the (φ(z), ψ(z)) phase plane. There are only two fixed points of (11)-(12): (0, 0) and (1, 0). Typically, with phase plane analysis of travelling wave solutions, we explore the possibility of a heteroclinic orbit between these fixed points by considering the local behaviour of the linearized system near the fixed points [1,23]. For the Porous-Fisher model, whose phase plane is identical to (11)-(12), such heteroclinic orbits can only occur for c ≥ 1/ √ 2 [1,24]. However, for the Porous-Fisher-Stefan model, we have the addition of the Stefan-like condition in (10). Therefore, travelling wave solutions of the Porous-Fisher-Stefan model do not correspond to a heteroclinc orbit between (0, 0) and (1, 0), but rather a trajectory between (1, 0) and another point that is determined by the moving boundary condition in (10). This particular trajectory, ψ(φ; c), corresponds to the travelling wave solution for a given c with lim To determine ψ(φ; c), we divide (12) by (11) and obtain As previously mentioned, physical travelling wave solutions in the Porous-Fisher model can only occur when c ≥ 1/ √ 2 [1,24], and numerical simulations of the Porous-Fisher-Stefan model indicate that all travelling waves have c < 1/ √ 2. As a result, we examine (13) in two limiting regimes: c → 0 + and c → 1/ √ 2 − .

Travelling wave solutions for c 1
We first consider the solution of (13) in the limit where 0 ≤ c 1 by expanding ψ(φ) as a regular perturbation expansion in c, i.e. ψ(φ) = V 0 (φ) + cV 1 (φ) + O(c 2 ). Substituting this expansion into (13) provides O(c) : The solution of (14)- (15) are Evaluating these expressions at φ = 0 gives a two-term approximation for the wave speed c as a function of κ, provided that c 1: where α = 36 .02. We note that since ψ = 2uu remains O(1) as φ → 0 + , this implies that u = O(u −1 ) as u → 0 + , confirming that we are examining a new class of travelling wave solutions with infinite slope at the moving boundary.

Comparison of travelling wave solutions
To validate our asymptotic approximations, we firstly examine the trajectory ψ(φ; c), which solves (13) for a given c, in the (φ, ψ)-phase plane. To determine ψ(φ; c), we solve (13) numerically using ode45 in MATLAB. In Figure 2a, the two-term approximation V 0 (φ) + cV 1 (φ), valid when c 1, agrees very well with ψ(φ; c) up to c = 0.2. Furthermore, in Figure 2b, we have good agreement between ψ(φ; c) and the two-term approximation Ψ 0 (φ) + εΨ 1 (φ), in the limit when c = 1/ √ 2 − ε up to ε = 0.2. Since the asymptotic approximations accurately describe ψ(φ; c) in the phase plane, we now examine the asymptotic relationship between c and κ. We can estimate c(κ) from the numerically-computed ψ(φ; c) by noting that c = −κ ψ(0; c)/2. Alternatively, we can use the numerical solutions of (6)- (8) to determine the speed of the moving front L(t), once the solution has settled to the travelling wave solution (see Appendix A). Figure 3 shows that these two numerical approaches to determine c(κ) are indeed equivalent. Furthermore, we see that both asymptotic approximations for c(κ), (17) and (22), agree with the numerically-determined c(κ) relationship. To minimize the error between c(κ) and its asymptotic approximations, we recommend that (17) be used for κ < 2 and (22) be used for κ > 20. For intermediate values of κ, it appears that a numerical solution is warranted.
Finally, we compare the asymptotic approximation of the moving front u(z) = u(x − ct), shown in (25), with the numerical solution of (6)-(8), after sufficient time has passed such that the travelling wave has formed (Appendix A). Despite the fact that the asymptotic approximation is only expected to match when κ is large, we see in Figure 4a that the asymptotic approximation of the shape of the moving fronts, shown implicitly in (25), matches the numerically-determined travelling wave solution very well, even when κ = 1. Furthermore, the asymptotic approximation in (25) agrees well with the numerically-computed travelling wave fronts for a variety of wave speeds (Figure 4b).

Conclusions
In this work, we consider new travelling wave solutions for the Porous-Fisher model with a moving boundary. This model has certain features that could be considered to be advantageous over other modifications of the Fisher-KPP model. In summary, the new travelling wave solutions: (i) move with a speed that is slower than the more standard Porous-Fisher model; (ii) never lead to population extinction; (iii) have compact support and a well-defined moving front, and (iv) the travelling wave profiles have an infinite slope at the moving front. Using travelling wave coordinates, we transform the model into a single nonlinear differential equation. The solution of this differential equation, corresponding to a particular trajectory in the associated phase plane, determines the travelling wave front and can be approximated using asymptotic analysis in two limiting parameter regimes. In both cases, we obtain a good approximation of this trajectory, which can also be used to relate the speed of the moving front to κ, a constant appearing in the Stefan-like condition. Finally, we determine a highly-accurate approximate form of the sharp-fronted travelling wave with infinite slope at the moving boundary, corresponding to solutions with wave speeds 0 ≤ c < 1/ √ 2. Further extensions of this Porous-Fisher-Stefan model can also be made. For instance, the asymptotic analysis performed in this work can be extended to include higher-order terms. Another possible extension could be to consider generalizing the nonlinear degenerate diffusivity function to D(u) = u n , for some constant n > 0 [6,10,17,19,26]. We leave both extensions for future consideration.

Acknowledgments
This work is supported by the Australian Research Council (DP170100474). The authors thank the two anonymous referees for their helpful comments.

Appendix A. Numerical Solution for the Porous-Fisher-Stefan model
A MATLAB implementation of the numerical solution of (6)-(8), described below, can be found at https://github.com/nfadai/Fadai_TW2019.
To numerically compute solutions of (6)-(8), we first approximate the semi-infinite domain (−∞, L(t)] as the finite domain [0, L(t)], provided that L 0 > 0. Additionally, we transform (6)-(8) to a fixed-space domain by setting ξ = x/L(t). Consequently,(6)-(8) becomes To compute the numerical solution of (A.1)-(A.3), we must specify values for L 0 , κ and φ(ξ, 0). We obtain numerical solutions of (A.1) on a uniformly-spaced mesh of ξ ∈ [0, 1], i.e. ξ i = i∆ξ, i = 0, . . . , N , where ∆ξ = 1/N . We denote φ(ξ i , t j ) = φ n i and L(t j ) = L j for convenience, where n ≥ 1 is the nth Picard iteration estimate at time t j . Therefore, to determine φ i , we use Here, φ p i is the solution of φ i at the previous timestep, t j−1 , and ∆t is the timestep. We identify the system (A.4)-(A.5) as a tridiagonal matrix in φ n i at time t j , which we can solve efficiently using the Thomas algorithm. This solution is stored as Φ n ; if max |Φ n − Φ n−1 | < δ, where δ is some user-specified tolerance, then the Picard loop terminates and we proceed to updating the moving boundary for the next timestep. Otherwise, n = n + 1, the solution Φ n is stored as φ n−1 i , and the Picard iteration loop is performed again.
New travelling wave solutions of the Porous-Fisher model with a moving boundary 10 From the fixed-boundary PDE, the Stefan-like condition at ξ = 1 is We approximate φ(ξ, t) as φ(ξ, t j ) during the small interval t ∈ [t j , t j+1 ], allowing us to explicitly solve (A.6) to give a closed form approximation for L(t): Therefore, evaluating this expression at t = t j+1 and using a second-order finite difference approximation of ∂ ξ φ(ξ, t j )| ξ=1 , we obtain the approximation . (A.8) With these updated values of L j+1 and Φ n = φ p i , we update t = t + ∆t and j = j + 1; we then repeat the computation to integrate through the next time increment. The algorithm terminates when t + ∆t > t f , where t f is the user-specified final time. Once sufficient time has passed that the solution settles towards a travelling wave, we expect that dL(t)/dt = c, so we fit a straight line to our numerical estimate of L(t) and use the slope of that line to provide an estimate of c.