Solve system of coupled differential equations using scipy's solve_bvp

Question

I want to solve a boundary value problem consisting of 7 coupled 2nd order differential equations. There are 7 functions, y1(x),...y7(x), and each of them is described by a differential equation of the form

d^2yi/dx^2 = -(1/x)*dyi/dx - Li(y1,...,y7) for 0 < a <= x <= b,

where Li is a function that gives a linear combination of y1,...,y7. We have boundary conditions for the first-order derivatives dyi/dx at x=a and for the functions yi at x=b:

dyi/dx(a) = Ai,
yi(b) = Bi.

So we can rewrite this as a system of 14 coupled 1st order ODEs:

dyi/dx = zi,
dzi/dx = -(1/x)*zi - Li(y1,...,y7),

zi(a) = Ai,
yi(b) = Bi.

I want to solve this system of equations using the Python function scipy.integrate.solve_bvp. However, I have trouble understanding what exactly should be the input arguments for the function as described in the documentation (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html).

The first argument that this function requires is a callable fun(x,y). As I understand it, the input argument y must be an array consisting of values of yi and zi, and gives as output the values of zi and dzi/dx. So my function would look like this (pseudocode):

def fun(x,y):
    y1, z1, y2, z2, ..., y7, z7 = y
    return [z1, -(1/x)*z1 - L1(y1,...,y7),
            ...,
            z7, -(1/x)*z7 - L7(y1,...,y7)]

Is that correct?

Then, the second argument for solve_bvp is a callable bc(ya,yb), which should evaluate the residuals of the boundary conditions. Here I really have trouble understanding how to define such a function. It is also not clear to me what exactly the arrays ya and yb are and what shape they should have?

The third argument is x, which is the 'initial mesh' with a shape (m,). Should x consist only of the points a and b, which is where we know the boundary conditions? Or should it be something else?

Finally the fourth argument is y, which is the 'initial guess for the function values at the mesh nodes', and has shape (n,m). Its ith column corresponds with x[i]. I guess that the 1st row corresponds with y1, the 2nd with z1, the 3rd with y2, etc. Is that right? Furthermore, which values should be put here? We could put in the known boundary conditions at x=a and x=b, but we don't know how the function looks like at any other points. Furthermore, how does this y relate to the function bc(ya,yb)? Are the input arguments ya,yb somehow derived from this y?

Any help with understanding the syntax of solve_bvp and its application in this case would be greatly appreciated.

user_3088138 · Answer

For the boundary condition, bc returns the residuals of the equations. That is, transform the equations so that the right side is zero, and then return the vector of the left sides. ya,yb are the state vectors at the points x=a,b.

For the initial state, it could be anything. However, "anything" in most cases will lead to failure to converge due to singular Jacobian or too large mesh. What the solver does is solve a large non-linear system of equations (see collocation method). As with any such system, convergence is most assured if you start close enough to the actual solution. So you can try for yi

constant functions with value Bi
linear functions with slope Ai and value Bi at x=b

with the values for zi using the derivatives of these functions. There is no guarantee that this works, especially as close to zero you have basis solutions close to the constant solution and the logarithm, the closer a is to zero, the more singular this becomes.

Solve system of coupled differential equations using scipy's solve_bvp

1 Answers