PaperlayerView exact revision on arXiv

Source-linked AI summary

Hybridizable Discontinuous Galerkin methods for the Keller-Segel model

Silvia Bertoluzza, Gabriella Bretti, Micol Pennacchio, Christophe Prud'homme

arXiv:2610.33260v2math.NAcs.NA

TL;DR

Keller–Segel solutions can be singular and sensitive to parameters and numerical step sizes, making reliable approximation challenging. This paper analyzes an HDG discretization with two cellular-flux degree choices and establishes convergence rates for them. The analysis and numerical experiments support those convergence findings, including tests in a blow-up regime.

  • Problem

    Keller–Segel solutions may be singular and sensitive to model parameters and numerical step sizes, while numerical schemes must capture sharp density growth without nonphysical oscillations or instability.

    Source: §1 · p. 2

  • Method

    The paper analyzes an HDG discretization of Keller–Segel with two choices for the cellular-flux polynomial degree and a mixed formulation with flux unknowns.

    Source: §1 · p. 2, §3.3 · p. 6

  • Results

    HDG(k,k) gives order k+1 convergence for all variables; HDG(k−1,k) gives order k+1 for all except the cellular flux, which converges at order k, with ε reductions in the mixed Neumann–Robin case.

    Source: §1 · p. 3

  • Takeaways & Limitations

    Numerical experiments confirm the convergence findings and also test the method in a blow-up regime.

    Source: §1 · p. 2, §5.1 · p. 25

    For HDG(k−1,k) with mixed Neumann–Robin boundaries, the hidden constant may depend on ε; constants may also depend on the data and specified solution norms.

    Source: §3.3 · p. 8

Abstract

from arXiv · show

We propose a hybridizable discontinuous Galerkin (HDG) finite element method to approximate the solution of the time-dependent Keller-Segel model, a drift-diffusion system describing chemotactic cell movement. We develop and analyze the corresponding HDG scheme, and establish theoretical convergence results for the proposed discretization. Finally, we present numerical experiments that confirm the theoretical findings.

1. Introduction

The Keller–Segel model captures chemotactic movement and aggregation, but its nonlinear dynamics and numerical sensitivity make reliable approximation challenging. This paper analyzes an HDG method for the model and establishes convergence results for two discretization choices.

1. Introduction

  • The analysis studies the bounded solution regime, while numerical tests also examine the method in the challenging blow-up regime.The biological application motivating the analysis has a slightly different formulation in which blow-up does not occur.Source: §1 · p. 2
  • Keller–Segel solutions can be singular and sensitive to model parameters and numerical step sizes, so approximations must capture sharp density growth without nonphysical oscillations or instability.For no-flux boundaries, the desired scheme should also preserve nonnegativity and total mass.Source: §1 · p. 2
  • The paper analyzes HDG for Keller–Segel, addressing the proliferation of degrees of freedom that can hinder discontinuous Galerkin methods.The formulation differs from earlier drift–diffusion HDG work in its time-dependent chemoattractant equation, mixed Neumann–Robin boundaries, and cellular-flux definition.Source: §1 · p. 2
  • In HDG(k,k), all variables converge at order k+1; in HDG(k−1,k), all except the cellular flux converge at order k+1 (or k+1−ε for mixed Neumann–Robin boundaries), while that flux converges at order k (or k−ε).The paper also reports numerical tests in the blow-up regime.Source: §1 · pp. 2–3

2. Mixed form of the Keller–Segel equations

The Keller–Segel system is recast in mixed form to prepare for its HDG discretization. The formulation introduces cellular and chemoattractant flux unknowns and allows mixed Neumann–Robin boundary conditions.

2. Mixed form of the Keller–Segel equations

  • The system is considered on a convex polygonal domain over a finite time interval with initial conditions.Source: §2 · p. 3
  • Mixed Neumann–Robin boundary conditions divide the boundary into Neumann and Robin portions, each formed from whole polygon edges.Either portion may be empty; setting the Robin data εφ to zero recovers the usual no-flux condition on the Neumann portion.Source: §2 · p. 3
  • The mixed formulation introduces two additional unknowns, j and ψ, as fluxes for the cellular density and chemoattractant equations.Source: §2 · p. 3

3. Hybridized discontinuous Galerkin discretization

The HDG discretization uses elementwise polynomial spaces, skeleton traces, and numerical fluxes to formulate the Keller–Segel system. Two scheme variants differ in the polynomial degree assigned to the cellular flux and its trace, with stated convergence guarantees under regularity assumptions.

3.1. Notation.

  • The notation defines a quasi-uniform, shape-regular triangulation and broken function spaces and operators for elementwise gradients, divergences, and normal traces.Source: §3.1 · pp. 4–5

3.2. HDG formulation of the KS equation.

  • HDG introduces independent skeleton traces for u and φ and uses polynomial test spaces to write the mixed system in elementwise form.The discretization assigns degree k to most unknowns, with j and the trace of u assigned degree k1∈{k−1,k}.Source: §3.2 · pp. 5–6

3.3. The HDG scheme.

  • The scheme uses polynomial degrees k1 for the cellular flux jh and its trace, and k for uh and the chemoattractant block; the two choices are k1=k−1 and k1=k.The cellular-flux definition differs from, where the diffusive flux is separate and the drift term remains in the conservation equation.Source: §3.3 · p. 6
  • Under the theorem’s regularity and parameter assumptions, the semidiscrete problem has a unique maximal solution that either persists to T or has an unbounded uh L8 norm at its finite endpoint.The convergence estimate’s hidden constants may depend on data and solution norms; in the mixed Neumann–Robin case with k1=k−1, they may also depend on ε.Source: §3.3 · pp. 7–8
  • The cited passage identifies Table 1 as reporting the convergence orders attained by both schemes, with the comparison made under sufficiently high smoothness of the continuous solution.Source: Table 1 · p. 8

4. Proof of Theorem 3.5

This section completes the proof of Theorem 3.5 by establishing the auxiliary and pointwise estimates needed for its error analysis, with arguments and bounds that distinguish the two polynomial-degree cases.

4. Proof of Theorem 3.5

  • Mixed boundary conditions make the bounds for k1 = k − 1 suboptimal, whereas the flux definition yields an optimal result for k1 = k.The proof also differs from the analogous result in because of the time derivative and mixed Neumann–Robin conditions.Source: §4 · p. 8

Scheme of the proof.

  • The proof constructs intermediate projections, bounds the chemoattractant error in terms of mesh size and u − uh, then uses Gronwall’s lemma and continuation to reach T.The two degree cases diverge only in the projection step; the constants generally depend on equation parameters and the L∞ norm of ∇φ.Source: Scheme of the proof · p. 9, §4 · p. 9

4.1. Preliminary bounds.

  • The preliminary analysis introduces discrete norms and seminorms and establishes an integration-by-parts identity for the HDG bilinear forms.The seminorm associated with Au differs from the norm associated with Aφ, which is a norm because a > 0.Source: §4.1 · p. 10

4.2. Approximation bounds.

  • Approximation bounds cover the stabilization-parametric HDG projection and L2 projection, including improved chemoattractant bounds when ψ = −ν∇φ.For γ sufficiently large, the HDG projection errors are bounded by mesh powers capped by k + 1; the stated corollary gives h^min(r−1,k+1) regularity dependence for φ and ψ.Source: §4.2 · p. 11

4.3. The elliptic projection onto Vh.

  • The elliptic projection is well posed for sufficiently large M, and its uniform-in-time error is bounded in terms of u, ∂tu, j, and ∂tj.The bound applies for M > M0 and h < h0 under Theorem 3.5’s assumptions, with approximation orders determined by the regularity indices.Source: §4.3 · pp. 12–13

4.4. Estimating Φ − Φh as a function of u − uh.

  • The chemoattractant error is bounded using the projection estimates, while the cell and flux error relative to the elliptic projection is controlled up to pTh before Gronwall yields a bound.The argument’s constants depend on the L∞ norms of ψ and u; the section defines pTh as an interval endpoint within the discrete solution’s existence interval.Source: §4.4 · pp. 14–15

4.4. Estimating Φ ´ Φh as a function of u ´ uh.

  • The estimate in subsection 4.4 follows after Lemma 4.7 bounds the final term in the preceding inequality.Source: §4.4 · p. 16

4.6. The continuation argument.

  • For sufficiently small h, the discrete solution’s maximal existence interval is [0,T], and the continuation argument completes the proof of Theorem 3.5.If the interval ended before T, the discrete solution would have to blow up there, contradicting the continuation result.Source: §4.6 · pp. 17–18
  • For h ≤ h*, the estimate bounds the φ and ψ errors and the skeleton trace discrepancy by R_{k1,k}(h).The bound is ||φ − φh||_{L∞(0,T;L2(Ω))} + ||ψ − ψh||_{L2(0,T;L2(Ω))} + h^(1/2)||φh − pφh||_{L2(0,T;L2(Eh))} ≲ R_{k1,k}(h).Source: §4.6 · p. 17

4.7. Proof of Lemma 4.6.

  • For k1 = k, the proof uses the stabilization-parametric HDG projector; for k1 = k − 1, it uses L2 projection and an Aubin–Nitsche-type argument for the remaining bound.The adjoint problem used in the latter case is well posed for sufficiently large M; under mixed Neumann–Robin boundary conditions, its solution has H^(2−ε) regularity for 0 < ε < 1/2.Source: §4.7 · pp. 18, 20, 21, 22
  • The proof invokes pointwise estimate (4.28) for u − uI.Source: §4.7 · p. 24

5. Numerical tests

The tests compare theoretical convergence against an exact solution and examine the scheme in a strongly aggregating blow-up regime. Upwind stabilization controls oscillations in that regime, with mass and nodal positivity preserved in the reported run.

5. Numerical tests

  • The experiments test polynomial degrees k = 1, 2 across HDG(k−1,k) and HDG(k,k) configurations, using implicit Euler time discretization.Unless otherwise stated, nonlinear problems use Picard iteration and stabilization parameters γ_u = γ_φ = 1.Source: §5 · p. 24

5.1. Example 1.

  • For HDG(k,k), all variables show optimal O(h^(k+1)) convergence; for HDG(k−1,k), u_h, φ_h, and ψ_h converge at O(h^(k+1)) while j_h converges at O(h^k).The exact-solution test uses spatial refinement, with temporal error reported as negligible on the finest mesh.Source: §5.1 · p. 25, Figure 1 · p. 26, Figure 2 · p. 26

5.2. Example 2: Blow-up dynamics with chemotaxis.

  • The second test targets qualitative robustness, not numerical analysis of blow-up formation, using a highly concentrated initial density with mass above the stated critical value.The initial mass is approximately 10π, above 8π for the chosen parameters.Source: §5.2 · pp. 25–26
  • Without upwind stabilization, HDG(2,2) develops spurious oscillations and negative cell-density values by t = 6×10^-5, stopping the simulation.The run uses h = 1/70 and dt = 10^-5.Source: §5.2 · p. 27, Figure 3 · p. 27
  • The method is modified by adding upwind stabilization to the numerical normal flux to control oscillations and loss of positivity.The stabilization separates diffusive and chemotactic upwind contributions.Source: §5.2 · p. 27
  • With upwind stabilization, the HDG(2,2) run reproduces increasing cell density and preserves total mass to round-off and positive values at all Lagrange nodes.The relative mass error is about 6×10^-16; the reported results also hold with the time step halved and quartered.Source: §5.2 · pp. 27–28

Appendix A. Properties of the projections

The appendix recalls properties of the projections introduced earlier in the paper.

Appendix A. Properties of the projections

  • This appendix section recalls properties of the projections introduced in Section 4.Source: Appendix A · p. 28

A.1. Properties of the HDG projection.

The excerpt continues the HDG projection proof using the numerical flux definition and projection properties, then concludes the proof.

A.1. Properties of the HDG projection.

  • The proof sets the projected triple as Π₀(W) and uses the numerical flux definition to obtain the ensuing identity.Source: §A.1 · p. 29
  • The displayed terms include a cell-gradient contribution and a boundary normal-flux contribution involving the projection error.Source: §A.1 · p. 29
  • The argument invokes polynomial-space membership to apply the projection property, then concludes the proof.The excerpt also contains a term involving the difference between u and its projection and the difference between w and its projected value.Source: §A.1 · p. 29
Loading 2610.33260v2…