SciPy: solve_bvp Problem 2nd Order Diff. Eq

Question

Trying to solve a 2nd order diff. eq. with 2 boundary conditions and nothing I try seems to work and I can't find a tutorial which includes all/similar terms to what I have in my expression and, at least to me, the scipy documentation doesn't really explain how to use solve_bvp clearly.

I have the equation: y'' + 2/r * y' = A * y + b * y^3 where y is a function of r.

I've rewritten it as follows:

y1 = y(r)

y2 = y1'

so that y2' = -2/r * y2 + y1(A + b * y1^2)

With conditions that y'(0) = 0, y(r=10) = constant

A, b and the constant are known.

I have the following code but it doesn't seem to work, and as said, I'm a little confused by the documentation so any help would be appreciated!

def fun(x, y):
    return np.vstack((y[2], -2/x*y[1]+y[0]*(A+b*y[0]*y[0])))

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


x = np.linspace(1, 10, 10)
ya = np.zeros((3, x.size))
yb = np.zeros((3, x.size))

sol_1 = solve_bvp(fun, bc, x, ya)
sol_2 = solve_bvp(fun, bc, x, yb)

Thanks!

==========EDIT========================= There is an analytical solution, for which I have calculated, it was just seeing if I could also find the same solution numerically, I think the major issue is that there are two separate regions for the solutions, one where r < R (in) and one where r > R (out). This leads to two different solutions (only valid in their respective domain) with the conditions that y_in(R) == y_out(R) and y_in'(R) == y_out'(R). Full 2-part solution where Radius = 1, a=99, b = 1 and constant = 1, y(inf) = constant

From Lutz Lehmann's solution, it gets the right shape (at least for the inner region, although not on the right scales).

I'm just unsure how you'd go about coding all the equating solutions and I suppose even getting their solutions in the first place, although Lutz's answer is an amazing point in the right direction. Thanks

Answer 1

Issues

• The order of the equation is 2, thus the dimension of the state vector is 2, the value is always y[0], the derivative y[1], there is no y[2], probably a remnant from the translation from Matlab.

• Also in the boundary conditions, there is no ya[2], the derivative value is ya[1], and the function value in the second one is yb[0].

• The initial solution guess has to have the same number 2 of state components.

• Why do you compute two solutions with identical data?

• Remark: It is not necessary to convert the return values into numpy types, the solver checks and converts anyway.

BVP framework with singularity treatment

The ODE is singular at r=0, so one has to treat the first segment in a special way. The mean value theorem gives

(y'(r)-y'(0))/r->y''(0)  for  r->0,

so that in that limit r->0 you get

3*y''(0) = a*y(0) + b*y(0)^3`. 

This allows to define the first arc as

y(r) = y0 + (a*y0 + b*y0^3)*r^2/6
y'(r) = (a*y0 + b*y0^3)*r/3

up to order r³ and r². So if you want precision 1e-9 in y(r), then the first segment should not be longer than 1e-3.

Instead of trying to eliminate y0 from the equations for y(h) and y'(h) to get a condition connecting ya[0] and ya[1], let the solver do also this work and add y0 as parameter to the system. Then the boundary conditions have 3 slots corresponding to the added virtual dimension for the parameter, which can be naturally filled with the equations y(h)=ya[0] and ya[1]=y'(h) and the right boundary condition.

Altogether you can define the system as

h = 1e-3;

def fun(r, y, p):
    return  y[1], -2/r*y[1]+y[0]*(a+b*y[0]*y[0]) 

def bc(ya, yb, p):
    y0, = p
    yh = y0 + y0*(a+b*y0*y0)*h*h/6;
    dyh = y0*(a+b*y0*y0)*h/3
    return ya[0]-yh, ya[1]-dyh, yb[0]-c


x = np.linspace(h, 10, 10)
ya = np.zeros((2, x.size))

sol = solve_bvp(fun, bc, x, ya, p=[1])
print(sol.message,f"y(0)={sol.p[0]}");
plt.plot(sol.x, sol.y[0]);

so that with example parameters a, b, c = -1, 0.2, 3 you get a convergent solver call with y(0)=2.236081087849196 and the resulting plot