Accurate simulation of pulled and pushed fronts in the nonautonomous Fisher-KPP equation

TL;DR

The paper tackles accurate simulation of front propagation in the nonautonomous Fisher-KPP equation on an effectively infinite domain by introducing a Green's function boundary condition (GBC) that couples a nonlinear interior to a linear leading-edge region. The method relies on a moving frame and a boundary Green's function with the constraint $c(t)=\gamma d(t)$ to yield exact half-line solutions, enabling accurate front velocity and profile measurements on small domains. It demonstrates that pulled fronts under time-dependent parameters can exhibit nonlinear deviations from linear-theory predictions, and shows that pushed fronts retain exponential convergence with high precision, while nonautonomous forcing can induce dynamic pushed–pulled transitions with measurable delays. Overall, GBC provides a robust numerical tool for analyzing front propagation in nonautonomous reaction-diffusion systems and clarifies complex interactions between parameter changes and front dynamics. The approach lays groundwork for rigorous analysis, buffer-size optimization, and extensions to more general settings beyond one dimension.

Abstract

We introduce a novel numerical method for direct simulation of front propagation in the Fisher-KPP equation with a time-dependent parameter on an infinite domain. The method computes a time-dependent boundary condition that accurately captures the leading-edge dynamics by coupling the nonlinear simulation region to a linear approximation region in which the dynamics can be solved exactly via the Green's function of the linearized equation. This approach enables precise front velocity measurements on relatively small computational domains for a variety of nonautonomous regimes and initial conditions for which existing numerical methods break down. We apply the method to pulled and pushed fronts in the Fisher-KPP equation with quadratic and quadratic-cubic nonlinearities, finding that it improves the accuracy of the simulated front velocity even for constant parameters and a fixed domain size. For pulled fronts with a diffusion coefficient that increases algebraically in time, our results reveal a deviation from the natural asymptotic velocity predicted by linear theory, whose explanation requires nonlinear theory. For pushed fronts with constant parameters, the method reproduces the exponential convergence to the theoretical asymptotic front speed and profile with improved precision. For a slowly time-varying linear growth parameter, we find that the pushed front velocity follows the changing parameter adiabatically if the asymptotic pushed velocity remains faster than the natural asymptotic pulled velocity. As the growth parameter moves toward the pushed--pulled transition point, the competition between the pushed and pulled fronts can result in both delayed and even premature onset of the pushed--pulled transition, depending on the form of parameter growth. The numerical method presented here proves to be an effective tool for analyzing front propagation in nonautonomous systems.

Figures (15)

• Figure 1: Naïve choices for the boundary condition do not accurately capture the dynamics at the leading edge of a front in the F-KPP equation. (a) The exponential leading edge of a front with steepness $\lambda$. A desirable "infinity-simulating" boundary condition would give correct values for $u(x_1)$ and $u_x(x_1)$ of this tail, where $x_1$ is the boundary of the direct numerical simulation. (b) A zero Dirichlet boundary condition $u(x_1)=0$ results in a front that is too steep and therefore too slow. (c) A zero Neumann boundary condition $u_x(x_1) = 0$ results in a front that is too shallow and therefore too fast.
• Figure 2: In the Green's (function) boundary condition (GBC) method, the domain is split into two regions: the nonlinear (NL) simulation region and the linear approximation (LA) region. These regions overlap in a buffer region of size $\delta$ required for regularization. At each time step, (1) the value $u(x_1-\delta, t)$ acts as a boundary forcing for the LA region to predict $u(x_1, t+dt)$. Then, (2) the updated value $u(x_1, t+dt)$ is the new boundary condition for an implicit step to compute $u(x,t+dt)$ within the NL region. The left boundary condition is a zero-flux condition $u_x(x_0, t)=0$.
• Figure 3: The GBC method benchmarks well for the simplest case of the autonomous Fisher equation \ref{['eq:fisher_time_independent']} with a step function initial condition. Simulations use $dt=0.025$, $t\in[0,1000]$, $dx=0.01$, $x\in[-30, 30]$, with a moving frame of speed $c=v^*=2$. A buffer region of size $\delta = 10$ is used for the GBC nonlinear-linear transition (see Fig. \ref{['fig:numerical_approach']}). (a) Lab-frame visualization of the propagating front. The solid, black curves at $t=0$, $t=20$, $t=40$, and $t=60$ are the numerically simulated profiles obtained with the GBC method. The dotted, gray line with solid points indicates the tracked front position at $u=0.1$. (b) Front position $x(t)$ in the lab frame with the GBC method (solid, black) and a naïve zero Dirichlet condition (solid, gray). The difference in this plot is not perceptible. (c) Front position $\Delta x(t) = x(t) - v^*t$ in the moving frame for the two cases. Here, we clearly see that the GBC front correctly follows the theoretical $\ln t$ scaling (fitted with a red dashed line), while the DZ (Dirichlet zero) front does not. (d) Front velocity $v(t)$ in the lab frame, computed by averaging over a time window $\Delta t = 5 = 200\,dt$. The asymptotic value $v^*=2$ (dashed, red) is approached from below in the GBC method, but the velocity undershoots in the DZ method. This slowdown in the DZ case is expected (cf. Fig. \ref{['fig:boundary_condition']}). (e) Front velocity deviation from the asymptotic value $\Delta v(t) = v(t) - v^*$ over $1/t$ approaches the $-3/2t$ line (dashed, red) as $t\to\infty$, in agreement with theory Bramson1983convergencevan_saarloos_front_2003avery_universal_2022.
• Figure 4: The GBC method also benchmarks well for Eq. \ref{['eq:fisher_time_independent']} with an exponential initial condition $e^{-\lambda_0|x|}$. Here, $\lambda_0=2/\sqrt{6}$ which is less than the critical steepness $\lambda^*=1$. We set the moving frame speed $c$ to be the asymptotic speed $v=\lambda_0 + 1/\lambda_0 = 5/\sqrt{6}$. (a) Numerically simulated front (solid, black) with an exponential initial condition propagates through the lab frame. This particular case has the known exact asymptotic solution given in Eq. \ref{['eq:shallow_exact']}. At $t=60$, the numerical solution shows good agreement with this exact solution (magenta, dashed). (b) Front velocity over time shows that the GBC front (solid, black) correctly attains the theoretical asymptotic velocity $v=5/\sqrt{6}$ (magenta, dashed), while the DZ front (solid, gray) eventually falls below $v^*=2$ (red, dashed). (c) Local steepness $\lambda = -u_x/u$ in the moving frame coordinate $z=x-\int_0^t c(t')\,dt'$ at $t=200$ shows that the GBC front steepness correctly approaches $\lambda_0 = 2/\sqrt{6}$ at the leading edge while the DZ front does not.
• Figure 5: The GBC method applied to the Fisher equation with linearly increasing time-dependent diffusion coefficient $d(t) = 1+\epsilon t$, $\epsilon = 0.001$. The numerical simulation parameters are $dt=0.025$, $t\in[0,1000]$, $dx=0.01$, $x\in[-90, 90]$ with a moving frame $c(t)=\gamma d(t)$ with $\gamma=1.65$ and a buffer region of size $\delta = 10$. (a) The numerical front velocity $v(t)$ (solid, black) deviates from the natural asymptotic velocity $v^{**}(t)$ (red, dashed) as $t\to\infty$ in disagreement with linear theory. (b) The velocity difference $\Delta v(t) = v(t) - v^{**}(t)$ reveals this disagreement more clearly. Simulations with $dt\in \{0.2, 0.1, 0.05, 0.025\}$ are shown in progressively darker shades of gray, with the $dt=0.025$ simulation depicted in (a) drawn in black. (c) A log-log plot of $\Delta v(t)$ suggests the existence of short- and long-time asymptotic algebraic regimes.