Search code examples
pythonpython-3.xnumerical-methods

scipy.integrate.solve_bvp for Neumann boundary conditions?


The question is what it says on the tin. I'm trying to numerically solve a boundary value problem and my friend is asking me whether the solver would work for these conditions. The page for the solver doesn't give the conditions we have i.e. bc(y(a),y(b), p) = 0 but the form of our question is y(0) = some constant value and y'(b) = 0, giving our Neumann conditions, Would you need to rewrite the function to have a first order reduction like in the shooting method?


Solution

  • I suppose that you are solving some second order ODE

    y''(t) = F(t,y(t), y'(t))
    

    For the first order system use the state vector u = [y, v] = [y, y']. Then the minimum to call the BVP solver is

    def ode(t,u): y,v = u; return [v, F(t,y,v)];
    def bc(u0, ub): return [u0[0]-y0, ub[1]];
    
    y0 = someconstantvalue;
    
    t_init = [0,b];
    u_init = [[y0, y0], [0, 0]];
    res = solve_bvp(ode, bc, t_init, u_init)
    

    Always check res.success or at least print out res.message. You might need a better initial guess, for non-linear problems different initial guesses can give different solutions. res.sol contains the solution as an interpolating function, for generating plots it is better to use this function than the sparse collection of internal steps in res.x and res.y