pith. sign in

arxiv: 2607.00713 · v1 · pith:OB7RCAYSnew · submitted 2026-07-01 · 🧮 math.NA · cs.NA

A linear, decoupled and positivity-preserving time-staggered block-centered finite difference method for the multi-species Keller-Segel chemotaxis system

Ao Zhang , Bingyin Zhang , Hongfei Fu This is my paper

Pith reviewed 2026-07-02 07:51 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords multi-species Keller-Segel systemblock-centered finite differencepositivity-preservingmass conservationsecond-order convergencetime-staggered discretizationnon-uniform gridschemotaxis
0
0 comments X

The pith

A prediction-then-projection block-centered finite difference scheme for the multi-species Keller-Segel system is positivity-preserving, mass-conserving, and second-order accurate on non-uniform grids.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops a linearly implicit, second-order block-centered finite difference prediction-then-projection scheme for the multi-species Keller-Segel chemotaxis system on non-uniform spatio-temporal grids. It combines Crank-Nicolson time marching with an L2 projection step and time-staggered discretization to decouple cell densities from the chemoattractant while enforcing positivity and mass conservation. Mathematical induction and energy analysis establish unique solvability together with second-order convergence of cell densities in the discrete L2 norm and of the chemoattractant in the discrete H1 norm. The variable time-step and grid features support adaptive resolution near blow-up without loss of the preserved properties.

Core claim

The proposed scheme is a linearly implicit, second-order block-centered finite difference prediction-then-projection method that integrates a standard Crank-Nicolson time-marching algorithm with an L2 projection step. On non-uniform spatio-temporal grids with variable time steps and time-staggered discretization, it fully decouples the multi-species cell density variables from the chemoattractant concentration. The scheme is uniquely solvable, positivity-preserving and mass-conserving, with cell densities converging at second order in the discrete L2 norm and the chemoattractant at second order in the discrete H1 norm, as proved by mathematical induction and energy analysis.

What carries the argument

The time-staggered block-centered finite difference discretization combined with an L2 projection step that enforces positivity and mass conservation while maintaining linear decoupling.

If this is right

  • The scheme is uniquely solvable at each time step.
  • Cell densities remain positive and both cell densities and chemoattractant conserve mass.
  • Cell densities converge at second order in the discrete L2 norm and the chemoattractant converges at second order in the discrete H1 norm.
  • Variable time steps and non-uniform grids permit adaptive resolution near blow-up while retaining positivity and mass conservation.
  • Time-staggered discretization decouples the multi-species densities from the chemoattractant, yielding a linear system at each step.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same projection-plus-staggering construction could be applied to other positivity-constrained reaction-diffusion systems that require decoupling for efficiency.
  • Non-uniform grids with this method may improve resolution of localized singularities in related pattern-formation models without custom mesh refinement techniques.
  • Variable time-step control could support stable continuation past initial aggregation events in long-time chemotaxis simulations.
  • The energy-analysis technique used for convergence might extend to other block-centered schemes on staggered grids.

Load-bearing premise

The L2 projection step combined with time-staggered discretization on non-uniform grids preserves the second-order accuracy of the Crank-Nicolson scheme without order reduction or loss of decoupling.

What would settle it

A blow-up test case in which the computed cell densities become negative, total mass is not conserved, or observed convergence rates in the discrete L2 and H1 norms fall below second order.

Figures

Figures reproduced from arXiv: 2607.00713 by Ao Zhang, Bingyin Zhang, Hongfei Fu.

Figure 1
Figure 1. Figure 1: Evolution of the extrema of u, v, c, and the total mass and iteration numbers for u and v (µ = 0) for Example 4.2. In the following simulation, the computational domain is discretized using M = 80 grid points in both the x- and y-directions, and the time stepsize in (4.1) is set to ∆t = 2.0×10−3 . The simulation results on both uniform grids (i.e., µ = 0) and non-uniform grids (i.e., µ = 0.1) are summarize… view at source ↗
Figure 2
Figure 2. Figure 2: Evolution of the extrema of u, v, c, and the total mass and iteration numbers for u and v (µ = 0.1) for Example 4.2. 0 2 4 6 8 10 t -35 -30 -25 -20 -15 -10 -5 0 Energy 0 2 4 6 8 10 t -35 -30 -25 -20 -15 -10 -5 0 Energy [PITH_FULL_IMAGE:figures/full_fig_p019_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Evolution of the energy µ = 0 (left) and µ = 0.1 (right) for Example 4.2. Under this scenario, the solutions of the 2D two-species Keller–Segel chemotaxis model (1.2) are expected to blow up in a finite time, as the total initial mass satisfies M[u 0 ] + M[v 0 ] ≈ 47.07 > 8π. Nevertheless, as long as the solution exists prior to the blow-up time, both positivity and mass conservation properties shall be pr… view at source ↗
Figure 4
Figure 4. Figure 4: Evolution of the extrema of u, v, c, and the total mass and iteration numbers for u and v for Example 4.3. 0 0.005 0.01 0.015 0.02 t 0 1 2 3 4 5 6 7 8 Maximum of u and v 106 Max u Max v 5 5.5 6 10-3 400 600 800 t = 0.0186 Max v = 3919000 t = 0.0186 Max u = 7837000 (a) maximum of u and v 0 0.005 0.01 0.015 0.02 t 1 2 3 4 5 6 7 8 9 10 Time stepsizes 10-4 0.0176 0.018 0.0184 2 4 6 10-5 t = 0.001274 Time steps… view at source ↗
Figure 5
Figure 5. Figure 5: Evolution of the maximum values of u, v, the time stepsizes, and energy for Example 4.3. In this test, we set M = 80 grid points in each spatial direction and adopt the grid partitions (4.4) along both the x- and y-directions. The adaptive time-stepping strategy (4.3) is employed 20 [PITH_FULL_IMAGE:figures/full_fig_p020_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Contour plots of u, v, c (from top to bottom) at time instants t = 0, 1.27 × 10−3 , 1.57 × 10−2 , 1.86 × 10−2 (from left to right) for Example 4.3. Example 4.4 (3D simulation). In the last example, we consider the 3D two-species Keller–Segel chemotaxis model (1.2) in a cubic domain Ω = (0, 1)3 . The initial conditions are prescribed as follows: u 0 (x) = 4 exp −100((x − 0.5)2 + (y − 0.5)2 + (z − 0.5)2 )  … view at source ↗
Figure 7
Figure 7. Figure 7: Evolution of the extrema of u, v, c, and the total mass and iteration numbers for u and v (µ = 0) for Example 4.4. 5. Conclusion This paper has introduced a fully decoupled, linearly implicit, positivity-preserving, and time￾staggered BCFD prediction-then-projection scheme for the multi-species Keller–Segel chemotaxis system. We demonstrate several key features of the proposed scheme: (i) The proposed sche… view at source ↗
Figure 8
Figure 8. Figure 8: Evolution of the extrema of u, v, c, and the total mass and iteration numbers for u and v (µ = 0.1) for Example 4.4. 0 0.05 0.1 0.15 0.2 t -2 0 2 4 6 8 10 12 Energy 0 0.05 0.1 0.15 0.2 t -2 0 2 4 6 8 10 12 Energy [PITH_FULL_IMAGE:figures/full_fig_p023_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Evolution of the energy µ = 0 (left) and µ = 0.1 (right) for Example 4.4. (iv) An adaptive time-stepping strategy (4.3), driven by the numerical solution evolution behavior, together with the time-staggered BCFD method on non-uniform spatial grids, effectively and accurately captures the blow-up phenomenon. Furthermore, extensive numerical experiments have validated the accuracy, positive-preserving and ma… view at source ↗
Figure 10
Figure 10. Figure 10: Slices of u, v, c (from top to bottom) at time instants t = 0, 8.2 × 10−3 , 4.0 × 10−2 , 0.2 (from left to right) for Example 4.4. dissipation property is unavailable and requires further investigation. CRediT authorship contribution statement Ao Zhang: Methodology, Formal analysis, Software, Writing-Original draft. Bingyin Zhang: Methodology, Formal analysis, Writing-Original draft. Hongfei Fu: Conceptua… view at source ↗
read the original abstract

In this paper, we present a linearly implicit, second-order block-centered finite difference (BCFD) prediction-then-projection scheme for the multi-species Keller-Segel chemotaxis system on non-uniform spatio-temporal grids. The proposed scheme integrates a standard Crank-Nicolson time-marching algorithm with an $L^2$ projection step to enforce positivity and mass conservation. The use of variable time stepsize and time-staggered discretization fully decouples the solutions of the multi-species cell density variables and the chemoattractant concentration variable while facilitating linearization, thereby greatly enhancing computational efficiency. Notably, the variable time-stepping algorithm and non-uniform grid BCFD discretization jointly enable adaptive resolution and local refinement near blow-up, thereby improving efficiency and accuracy without compromising the desired physical property-preserving in the simulation. Furthermore, using the mathematical induction method and the energy analysis approach, the unique solvability of the proposed scheme is rigorously proved, and we show that cell densities achieve second-order convergence in both time and space in the discrete $L^2$ norm, while the chemoattractant concentration achieves second-order convergence in the discrete $H^1$ norm. Representative numerical experiments are presented to validate the theoretical findings and demonstrate the reliability of the proposed scheme in simulating the blow-up phenomenon.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

1 major / 0 minor

Summary. The paper proposes a linearly implicit, second-order block-centered finite difference (BCFD) prediction-then-projection scheme for the multi-species Keller-Segel chemotaxis system on non-uniform spatio-temporal grids. It combines Crank-Nicolson time marching with time-staggered discretization to fully decouple the cell-density and chemoattractant equations while linearizing the system, followed by an L² projection step that enforces positivity and mass conservation. The authors prove unique solvability, positivity preservation, and mass conservation, and establish second-order convergence in the discrete L² norm for cell densities and in the discrete H¹ norm for the chemoattractant concentration, using mathematical induction together with energy estimates. Numerical experiments are presented to illustrate the scheme’s behavior near blow-up.

Significance. If the claimed convergence rates are rigorously established without order reduction, the method supplies an efficient, structure-preserving discretization that supports adaptive non-uniform grids and variable time steps for multi-species chemotaxis models. The decoupling and linearization reduce computational cost while preserving key physical properties, which is valuable for long-time simulations of blow-up phenomena.

major comments (1)
  1. [Convergence analysis (induction and energy estimates)] In the energy analysis and induction argument establishing second-order convergence (abstract and the convergence theorem), the L² projection step is asserted to preserve the O(τ² + h²) accuracy of the Crank-Nicolson predictor. On non-uniform grids the projection operator’s approximation properties in the discrete H¹ seminorm must be controlled explicitly; without a bound showing that the projection error remains O(τ² + h²) uniformly in the staggered-time, variable-grid setting, the inductive step for the chemoattractant H¹ estimate does not close at second order.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for the careful reading and constructive comments on our manuscript. We address the single major comment below.

read point-by-point responses

  1. Referee: In the energy analysis and induction argument establishing second-order convergence (abstract and the convergence theorem), the L² projection step is asserted to preserve the O(τ² + h²) accuracy of the Crank-Nicolson predictor. On non-uniform grids the projection operator’s approximation properties in the discrete H¹ seminorm must be controlled explicitly; without a bound showing that the projection error remains O(τ² + h²) uniformly in the staggered-time, variable-grid setting, the inductive step for the chemoattractant H¹ estimate does not close at second order.

    Authors: We agree that an explicit bound on the L² projection error in the discrete H¹ seminorm is needed to close the induction on non-uniform grids. The original proof invoked standard approximation properties of the projection within the BCFD setting but did not derive the required uniform estimate for the staggered-time, variable-grid case. In the revised manuscript we will add a dedicated lemma establishing that the projection error remains O(τ² + h²) in the discrete H¹ seminorm, using the block-centered structure and the non-uniform mesh regularity assumptions already present in the paper. This addition will allow the energy estimates and inductive step for the chemoattractant to close at second order without altering the overall proof architecture. revision: yes

Circularity Check

0 steps flagged

No significant circularity in the derivation chain

full rationale

The paper constructs a prediction-then-projection scheme from standard Crank-Nicolson discretization and L2 projection operators on non-uniform grids, then proves unique solvability, positivity, mass conservation, and second-order convergence via mathematical induction plus energy estimates. No equations reduce to fitted inputs renamed as predictions, no self-definitional loops appear in the scheme definition or convergence claims, and no load-bearing steps rely on self-citations or imported uniqueness theorems. The central claims rest on explicit inductive arguments and energy analysis that are independent of the target results, rendering the derivation self-contained.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The scheme rests on standard finite-difference properties and the assumption that the projection operator preserves positivity and mass without order reduction; no free parameters or new entities are introduced in the abstract.

axioms (2)
  • standard math Block-centered finite difference operators on non-uniform grids satisfy the discrete summation-by-parts identities needed for energy estimates.
    Invoked implicitly for the convergence analysis.
  • domain assumption The L2 projection step onto the positive cone preserves the second-order truncation error of the Crank-Nicolson scheme.
    Required for the stated convergence rates to hold after projection.

pith-pipeline@v0.9.1-grok · 5772 in / 1447 out tokens · 29123 ms · 2026-07-02T07:51:29.156041+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

30 extracted references · 30 canonical work pages

  1. [1]

    Keller, L

    E. Keller, L. Segel, Initiation of slide mold aggregatio n viewed as an instability, J. Theor. Biol. 26 (1970) 399–415

    work page 1970

  2. [2]

    Keller, L

    E. Keller, L. Segel, Model for chemotaxis, J. Theor. Biol . 30 (1971) 225–234

    work page 1971

  3. [3]

    Chertock, Y

    A. Chertock, Y. Epshteyn, H. Hu, A. Kurganov, High-order positivity preserving hybrid finite- volume-finite-difference methods for chemotaxis systems, A dv. Comput. Math. 44 (2018) 327– 350. 24

    work page 2018

  4. [4]

    Slotboom, Computer-aided two-dimensional analysis of bipolar transistors, Electron Devices 20 (1973) 669–679

    J. Slotboom, Computer-aided two-dimensional analysis of bipolar transistors, Electron Devices 20 (1973) 669–679

    work page 1973

  5. [5]

    S. jin, B. Yan, A class of asymptotic-preserving schemes for the Fokker-Planck-Landau equation, J. Comput. Phys. 230 (2011) 6420–6437

    work page 2011

  6. [6]

    J. Hu, X. Zhang, Positivity-preserving and energy-diss ipative finite difference schemes for the Fokker–Planck and Keller–Segel equations, IMA J. Numer. An al 43 (2023) 1450–1484

    work page 2023

  7. [7]

    J. Liu, L. Wang, Z. Zhou, Positivity-preserving and asym ptotic preserving method for 2D Keller– Segel equations, Math. Comp. 87 (2018) 1165–1189

    work page 2018

  8. [8]

    Cheng, J

    Q. Cheng, J. Shen, A new Lagrange multiplier approach for constructing structure preserving schemes, I. Positivity preserving, Comput. Methods Appl. M ech. Engrg. 391 (2022) 114585

    work page 2022

  9. [9]

    Huang, J

    F. Huang, J. Shen, Bound/positivity preserving and ener gy stable scalar auxiliary variable schemes for dissipative systems: applications to Keller–S egel and Poisson–Nernst–Planck equa- tions, SIAM J. Sci. Comput. 43 (2021) A1832–A1857

    work page 2021

  10. [10]

    F. Tong, Y. Cai, Positivity preserving and mass conserv ative projection mthod for the Poisson– Nernst–Planck equation, SIAM J. Numer. Anal. 62 (4) (2024) 2 002–2024

    work page 2024

  11. [11]

    Arbogast, M

    T. Arbogast, M. Wheeler, I. Yotov, Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences, SIAM J. Nume r. Anal. 34 (1997) 828–852

    work page 1997

  12. [12]

    H. Rui, H. Pan, A block-centered finite difference method for the Darcy–Forchheimer model, SIAM J. Numer. Anal. 50 (2012) 2612–2631

    work page 2012

  13. [13]

    H. Rui, W. Liu, A two-grid block-centered finite differen ce method for Darcy–Forchheimer flow in porous media, SIAM J. Numeri. Anal. 53 (2015) 1941–1962

    work page 2015

  14. [14]

    J. Xu, S. Xie, H. Fu, A two-grid block-centered finite diff erence method for the nonlinear regularized long wave equation, Appl. Numer. Math. 171 (202 2) 128–148

    work page

  15. [15]

    X. Wang, J. Xu, H. Fu, A linearlized mass-conservative f ourth-order block-centered finite differ- ence method for the semilinear Sobolev equation with variab le coefficients, Commun. Nonlinear Sci. Numer. Simul. 130 (2024) 107778

    work page 2024

  16. [16]

    Y. Shi, S. Xie, D. Liang, K. Fu, High order compact block- centered finite difference schemes for elliptic and parabolic problems, J. Sci. Comput. 87 (2021) 8 6

    work page 2021

  17. [17]

    X. Li, J. Shen, H. Rui, Energy stability and convergence of SA V block-centered finite difference method for gradient flows, Math. Comp. 88 (2019) 2047–2068

    work page 2019

  18. [18]

    J. Xu, H. Fu, A decoupled linear, mass-conservative blo ck-centered finite difference method for the Keller–Segel chemotaxis system, J. Comput. Phys. 526 (2 025) 113775

    work page

  19. [19]

    Zhang, J

    H. Zhang, J. Wang, X. Pan, A non-iterative fully decoupl ed second-order projection method with staggered time discretization for Keller–Segel–Navi er–Stokes system, J. Comput. Phys. 545 (2026) 0021–9991

    work page 2026

  20. [20]

    Weiser, M

    A. Weiser, M. Wheeler, On convergence of block-centere d finite differences for elliptic problems, SIAM J. Numer. Anal. 25 (1988) 351–375

    work page 1988

  21. [21]

    H. Rui, H. Pan, Block-centered finite difference methods for parabolic equation with time- dependent coefficient, Jpn. J. Ind. Appl. Math. 30 (2013) 681– 699

    work page 2013

  22. [22]

    Berikelashvili, M

    G. Berikelashvili, M. Gupta, M. Mirianashvili, Conver gence of fourth order compact difference schemes for three-dimensional convection-diffusion equat ions, SIAM J. Numer. Anal. 45 (2007) 443–455. 25

    work page 2007

  23. [23]

    Dawson, M

    C. Dawson, M. Wheeler, C. Woodward, A two-grid finite diff erence scheme for nonlinear parabolic equations, SIAM J. Numer. Anal. 35 (1998) 435–452

    work page 1998

  24. [24]

    Plemmons, M-matrix characterizations

    R. Plemmons, M-matrix characterizations. I—nonsingul ar M-matrices, Linear Algebra Appl. 18 (1977) 175–188

    work page 1977

  25. [25]

    Huang, J

    X. Huang, J. Shen, Efficient numerical schemes for a two-s pecies Keller–Segel model and inves- tigation of its blowup phenomena in 3D, Acta Appl. Math. 190 ( 2024) 10

    work page 2024

  26. [26]

    Espejo, K

    E. Espejo, K. Vilches, C. Conca, Sharp condition for blo w-up and global existence in a two species chemotactic Keller–Segel system in R2, Eur. J. Appl. Math. 24 (2012) 297–313

    work page 2012

  27. [27]

    Lin, The fully parabolic multi-species chemotaxis s ystem in R2, Eur.J

    K. Lin, The fully parabolic multi-species chemotaxis s ystem in R2, Eur.J. Appl. Math. 35 (2024) 675–706

    work page 2024

  28. [28]

    K. Wang, E. Liu, X. Feng, Optimal error estimate of uncon ditionally positivity-preserving, mass- conserving and energy stable method for the Keller–Segel ch emotaxis model, Math. Comp. 94 (2025) 2761–2793

    work page 2025

  29. [29]

    J. Shen, J. Xu, Unconditionally bound preserving and en ergy dissipative schemes for a class of Keller–Segel equations, SIAM J. Numer. Anal. 58 (2020) 1674 –1695

    work page 2020

  30. [30]

    Acosta-Soba, F

    D. Acosta-Soba, F. Guillén-González, J. Rodríguez-Ga lván, An unconditionally energy stable and positive upwind DG scheme for the Keller–Segel model, J. Sci. Comput. 97 (2023) 18. 26

    work page 2023