How do you feed Scipy's bvp only the BCs you have?

Question

The only example/docs I can find are on the Scipy docs page.

To test, I'm looking at a time-independent Schrod eq in a 1d infinite potential well. This has a neat analytic solution found by solving the DE, and inserting boundary conditions of ψ(0) = 0, ψ(L) = 0, and that the func soln to 1, but this question applies to solving any DE where the BCs we know aren't for the initial value.

You can solve it numerically with Scipy's solve_ivp by starting with ψ(0) = 0, and cheating to place ψ'(0) appropriately using the analytic soln. Can use shooting method to find an appropriate E value, eg the normalization condition above.

These are two sets of BCs: ψ(0) = 0 for both, normalization for both, and a second value of ψ for the analytic approach, and an initial value of ψ' for the ivp approach. Scipy's solve_bvp seems to offer a solution using the first set of BCs numerically (since we're cheating by inserting ψ'), but i can't get it working. This pseudocode describes the problem, and is how I expect the API to behave:

bcs = {0: (0, None), L: (0, None)} # Two BCs on ψ; no BCs on derivative
x_span = (0, L)

sol = solve_bvp(rhs, bcs, x_span)

In reality, the code looks something like this, and I can't get it to work:

def bc(ψ_a, ψ_b):
    return np.array([ψ_a[0], ψ_b[0]])

x_span = (0, L)
x_eval = np.linspace(x_span[0], x_span[1], int(1e5))

x_guess = np.array([0, L])
ψ_guess = np.array([[0, 1], [0, -1]])

res = solve_bvp(rhs_1d, bc, x_guess, ψ_guess)

I've no idea how to build the bc function, and don't know why the guesses are set up the way they are. And unsure how I can guess for the value of ψ without also inserting a guess for ψ'. (The docs imply you can) Also of note, the docs shows an example implying you can use solve_bvp for a normalization BC as well, but not sure how to approach. (Example is too sparse)

The equivalent and working ivp code, for ref: (Compare to my solve_bvp pseudocode)

Python code:

ψ_0 = (0, sqrt(2/L) * n*π/L)
x_span = (0, L)

sol = solve_ivp(rhs_1d, x_span, ψ_0)

Answer 1

For the eigenvalue problem

-u''+V(x)u = c*u

with boundary conditions

u(0)=0=u(L)

and normalization

int(u(x)^2, x=0 to L)=1 

set up the integral as third component. With the eigenvalue as parameter these are 4 dimensions allowing for 4 boundary conditions, the additional 2 are that the integral at 0 is zero and that the integral at L has value 1.

# some length
L = 10;

# some potential function
def V(x): return 1+(2*x-L)**2;

# the ODE function
def odesys(x,y,p):
    u,v,S = y; c=p[0]
    return [v, (V(x)-c)*u , u**2 ]

# the boundary conditions
def boundary(y0, yL, c):
    return [ y0[0], yL[0], y0[2], yL[2]-1 ]

With the initial guess you select approximately what eigenfunction/eigenvalue you will get, more or less.

n=11;
w = (np.pi*n)/L
x_init = np.linspace(0,L,4*n+1);
u_init = np.sin(w*x_init);
v_init = np.cos(w*x_init)*w;
y_init = [ u_init, v_init, x_init/L ]

There is no need to put too many points into the guess, just enough that the structure of the first component is faithfully represented.

Then call the solver with the prepared data, take notice that the default tolerance is 1e-3, if you want better you have to allow for a finer subdivision. If everything runs fine, plot the solution.

res = solve_bvp(odesys, boundary, x_init, y_init, p=[w**2], max_nodes=10000, tol=1e-6)
print res.message

if res.success:
    x_disp = np.linspace(0,L,3001)
    y_disp = res.sol(x_disp)
    plt.plot(x_disp, y_disp[0])
    plt.title("eigenfunction to eigenvalue $\lambda=%.6f$"%res.p[0]);
    plt.grid(); plt.show()