Numerical challenges for the understanding of localised solutions with different symmetries in non-local hyperbolic systems

We consider a one-dimensional nonlocal hyperbolic model introduced to describe the formation and movement of self-organizing collectives of animals in homogeneous 1D environments. Previous research has shown that this model exhibits a large number of complex spatial and spatiotemporal aggregation patterns, as evidenced by numerical simulations and weakly nonlinear analysis. In this study, we focus on a particular type of localised patterns with odd/even/no symmetries (which are usually part of snaking solution branches with different symmetries that form complex bifurcation structures called snake-and-ladder bifurcations). To numerically investigate the bifurcating solution branches (to eventually construct the full bifurcating structures), we first need to understand the numerical issues that could appear when using different numerical schemes. To this end, in this study, we consider ten different numerical schemes (the upwind scheme, the MacCormack scheme, the Fractional-Step method, and the Quasi-Steady Wave-Propagation algorithm, combining them with high-resolution methods), while paying attention to the preservation of the solution symmetries with all these schemes. We show several numerical issues: first, we observe the presence of two distinct types of numerical solutions (with different symmetries) that exhibit very small errors, which might initially suggest that we have reached a steady-state solution, but this is not the case (this also implies an extremely slow convergence); second, in some cases, none of the investigated numerical schemes converge, posing a challenge for the development of numerical continuation algorithms for nonlocal hyperbolic systems; lastly, the choice of the numerical schemes, as well as their corresponding parameters such as time-space steps, exert a significant influence on the type and symmetry of bifurcating solutions. To conclude we emphasize that if we want to construct numerically bifurcation diagrams for these localised solutions with different symmetries, the resulting bifurcations may vary when different numerical schemes and/or corresponding parameters are employed.