Singular BVP with integrable boundaries

[NoFishLikeIan]

Issue

In many applications, one would want to solve boundary value problems (BVP), with singular boundary conditions. For example, in the study of reputation games in continuous time (Faingold & Sannikov, 2011), one finds equations of the form

$$ u' = \frac{v}{t (1 - t)} $$

$$ v' = \frac{1}{t (1 - t)} \left(v + \frac{u - w(t, v)}{\gamma(v)^2}\right) $$

on $t \in [0, 1]$ with $v(0) = v(1) = 1$, and $u(0)$ and $u(1)$ finite. Currently, defining this class of problems as a TwoPointsBVProblem and employing a shotting or a RK solver yields unstable results.

Possible Solution

I think MATLAB can handle these using a finite difference scheme, in the HOFiD_bvp (Hohenegger et al., 2024).

Context

This issue was prompted by this discord thread.

[ChrisRackauckas-Claude]

Thanks for raising this issue. Handling singular BVPs is indeed an important capability.

The MATLAB bvp4c and bvp5c solvers have built-in support for a class of singular BVPs of the form:

$$y' = \frac{S}{x} y + f(x, y)$$

where S is a constant matrix and the problem is posed on [0, b]. They handle this via a SingularTerm option that separates the singular part from the regular ODE evaluation.

Key papers describing these methods:

1. bvp4c (ACM TOMS 2001): Kierzenka, J. and Shampine, L.F. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27(3), 2001, pp. 299–316.

   • https://dl.acm.org/doi/10.1145/502800.502801

2. bvp4c Tutorial (2004): Shampine, L.F., Reichelt, M.W., and Kierzenka, J. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c."

   • PDF: https://classes.engineering.wustl.edu/che512/bvp_paper.pdf

3. bvp5c (JNAIAM 2008): Shampine, L.F. and Kierzenka, J. "A BVP Solver that Controls Residual and Error." J. Numer. Anal. Ind. Appl. Math. Vol. 3(1-2), 2008, pp. 27–41.

   • PDF via ResearchGate: https://www.researchgate.net/publication/228748602_A_BVP_solver_that_controls_residual_and_error

The singular BVP approach requires that the singularity be at the left endpoint and uses the smoothness condition S·y(0) = 0. For problems like the one you describe with singularities at both endpoints (t=0 and t=1), one approach is to split the domain and match solutions in the middle, or to reformulate the problem.

We could implement the SingularTerm style interface in BoundaryValueDiffEq.jl following the approach in these papers.