Solving two coupled second order boundary value problems

Question

I have solved a single second order differential equation with two boundary conditions using the module solve_bvp. However, now I am trying to solve the system of two second order differential equations;

U'' + a*B' = 0

B'' + b*U' = 0

with the boundary conditions U(+/-0.5) = +/-0.01 and B(+/-0.5) = 0. I have split this into a system of first ordinary differential equations and I am trying to use solve_bvp to solve them numerically. However, I am just getting arrays full of zeros for my solution. I believe I am implementing the boundary conditions wrong. It is not clear to me how to handle more than two equations from the documentation. My attempt is below

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import solve_bvp
alpha = 1E-8
zeta = 8E-3
C_k = 0.05
sigma = 0.01
def fun(x, y):

    return np.vstack((y[1],-((alpha)/(C_k*sigma))*y[2],y[2], -(1/(C_k*zeta))*y[1]))

def bc(ya, yb):

    return np.array([ya[0]+0.001, yb[0]-0.001,ya[0]-0, yb[0]-0])

x = np.linspace(-0.5, 0.5, 5000)
y = np.zeros((4, x.size))
print(y)


sol = solve_bvp(fun, bc, x, y)
print(sol)

In my question I have just relabeled a and b, but they're just parameters that I input. I have the analytic solution for this set of equations so I know one exists that is non-trivial. Any help would be greatly appreciated.

Answer 1

It is most times really helpful if you state at least once in a comment or by assignment to specifically named variables how you want to compose the state vector.

By the form of the derivative return vector, I would think you intend

U, U',  B, B'

which means that U=y[0], U'=y[1] and B=y[2],B'=y[3], so that your derivatives vector should correctly be

return y[1], -((alpha)/(C_k*sigma))*y[3],  y[3], -(1/(C_k*zeta))*y[1]

and the boundary conditions

return ya[0]+0.001, yb[0]-0.001, ya[2]-0, yb[2]-0

Especially your boundary condition should throw the algorithm in the first step because of a singular Jacobian, always check the .success field and the .message field of the solution structure.

---

Note that by default the absolute and relative tolerance of the experimental solve_bvp is 1e-3, and the number of nodes is limited to 500.

Setting the initial node number to 50 (5000 is much too much, the solver refines where necessary), and the tolerance to 1-6, I get the following solution plots that visibly satisfy the boundary conditions.