Solving basic barystochrone problem in python

Question

I am trying to solve $\frac{u''}{1+u'^2} - \frac{1}{2(1-u)} = 0$ subject to $u(0)=1, u(1)=0$.

If I understand how to do this properly, I first do the variable substitutions:

$u = y$, $y_1 = y; y_2 = y'$

yielding:

$$\frac{y_2'}{1+y_2^2} - \frac{1}{2(1-y_1)} = 0 \iff y'_2 = \frac{1 + y_2^2}{2(1-y_1)}$$

And the boundary conditions give $y(0) + y(1) - 1 = 0$

To solve the transformed problem I tried using scipy this way:

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

def fun(x, y):
    denom = 2.0 * (1.0-y[0])
    num = 1.0 + y[1]**2
    result = num / denom

    return np.vstack((y[1], result))

def bc(ya, yb):
    return ya + yb - 1.0

x = np.linspace(0, 1, 10, endpoint=False)
y = np.full((2, x.shape[0]), 0.01)

res = solve_bvp(fun, bc, x, y, verbose=2)

x_plot = np.linspace(0, 1, 100)
y_plot = res.sol(x_plot)
e_plot = 1.0 - np.sqrt(-(x_plot - 2.0)*x_plot)

plt.plot(x_plot, y_plot[0], label='$y$  approx')
plt.plot(x_plot, e_plot, label='$y$ exact')
plt.legend(loc='lower right')
plt.show()

But that is giving me:

Which is not even close to correct.

Answer 1

Your analytical solution $$ u = 1 - \sqrt{x (2-x)} $$ does not solve your original ODE. If you try plugging it into the original ODE, you get $$ \frac{u''}{1+u'^2} - \frac{1}{2 (1-u)} = \frac{1}{2 \sqrt{x(2-x)}} $$ Your solution does solve $$ \frac{u''}{1+u'^2} - \frac{1}{1-u} = 0 $$

The problem with your use of solve_bvp is that your boundary conditions are not the boundary conditions you want. The two boundary conditions you provided were $$ y_1(0) + y_1(1) - 1 = 0\\ y_2(0) + y_2(1) - 1 = 0 $$ Technically these two are true for the modified ODE, however they do leave open other solutions, and solve_bvp found one of those. I have no idea if the second BC would fit the original ODE.

What you really want are $$ y_1(0) - 1 = 0\\ y_1(1) = 0 $$

So the bc function should be

def bc(ya, yb):
    return np.array([ya[0]-1,yb[0]])

Here's the numerical solution I get with the correct BCs, both for your original ODE and the one your analytical solution actually solves.

edit note:

as pointed out by Lutz Lehmann you also don't specify x correctly as it's missing the correct endpoint, giving the wrong BC on the right side. That alone doesn't seem enough to find a reasonable solution, however it does slightly change the numerical solution. There's also potentially some numerical challenges on the left BC where $u=1$ causes a divide by zero. This doesn't seem to affect the solution too badly in this case, however it is something to be aware of.