Numerical or exact solution for a system of differential algebraic equations

Question

I am looking for ways to solve the following system of boundary value implicit ODEs over the interval $t \in [0, 1]$: \begin{equation} \lambda (Fg - Gf)^3 + 4 FGfg(g-f) = 0 \\ fg(Fg-Gf) + 2FG(gf' - fg') = 0, \end{equation} where $f = F'$ and $g=G'$ and $f, -g \geq 0$. The boundary values are $F(0) = G(1) = 0$ and $F(1) = G(0) = 1$.

The equations arise as the Euler Lagrange equations of the following constrained optimization problem: \begin{equation} \text{minimize}_{F, G} \int_0^1 \frac{dx}{F(x) / f(x) - G(x) / g(x)}, \end{equation} subject to: \begin{equation}\frac{1}{2} \int_0^1 [f(x) \log \frac{2f(x)}{f(x) - g(x)} - g(x) \log \frac{2g(x)}{g(x) - f(x)}] dx = J_0 \end{equation}

The constraint is simply the Jensen Shannon divergence of the two distributions on $[0, 1]$, whose densities are given by $f$ and $-g$ respectively. The objective function is another similarity function between the two distributions related to the Jaccard distance. Note that only the minimization problem above is nontrivial. The maximization problem can be solved with a pair of distributions on a 2-point space. Finite dimensional numerical solutions suggest the minimizing $F$ and $G$ are strictly monotonic over $[0, 1]$, thus describe distributions on an infinite space.

Any analytic insight (not necessarily leading to a solution) would be greatly appreciated. Alternatively, suggestion on what commercial or open source solvers to use would be just as awesome: I have looked into implicit Runge-Kutta implementation in R and matlab but they seem pretty limited in scope, in particular, only capable of solving index-1 initial value problems, unless there is a way to convert BVP to IVP that I missed?

Answer 1

If you assume that $f>0$ and $g < 0$ on $[0,1]$, then one can integrate the equations explicitly.

Assume $0<t<1$, so that $F$ and $G$ are positive in $(0,1)$. Let $p = f/F = (\log F)' >0$ and $q = g/G = (\log G)' <0$. Moreover, we have $p' = f'/F- p^2$ and $q' = g/G -q^2$. Then the second equation, after dividing by $(FG)^2$, reduces to $$ 2(qp'-pq')-pq(q-p) = 0. $$ Dividing this equation by $pq<0$, this becomes $$ 2p'/p-2q'/q-q+p=0, $$ which can be logarithmically integrated to yield $$ p^2F/(q^2G) = c^2 $$ for some constant $c >0$. Consequently $p\sqrt F + c\,q\sqrt G = 0$ (since $p>0$ and $q <0$), and this can be integrated again to yield $$ \sqrt F + c\sqrt G = a $$ for some constant $a$. Now, the boundary conditions imply that $a=1$ and $c=1$, so we have $$ \sqrt F + \sqrt G = 1. $$ If we now write $F = x^2$ and $G = (1-x)^2$ for some function $x$, we have $f = 2xx'$ and $g=-2x'(1-x)$, we can now substitute these into the first equation, which becomes $$ -8(\lambda-4)\,x^3(1-x)^3 (x')^3 = 0. $$ Thus, there is no solution satisfying the boundary conditions if $\lambda\not=4$, while, if $\lambda=4$, there are infinitely many solutions, namely $$ (F,G) = \bigl(x^2,(1-x)^2\bigr) $$ where $x$ is any strictly increasing function on $[0,1]$ that satisfies $x(0)=0$ and $x(1) = 1$.