pythonscipynumerical-methodsnumerical-integration

SciPy's solve_bvp going solving for the wrong parameter values

I am solving a ordinary differential equation with scipy's solve_bvp to get a family of eigenstates phi_n and eigenvalues omega_n. The equation looks something like y'' = -(omega_n - f(x)) ** 2 * y, and I want it to obey the boundary conditions y(0)=y(1)=0. f(x) is a given function that depends on some physical parameters.

The problem is that for certain 'critical' f(x) the system starts behaving very weirdly, and a bit unstable. Near these instabilities, solve_bvp either doesn't converge, or the omega_n go to wrong eigenvalues. For example if the eigenvalues are pi, 2pi, 3pi, ... and I input as a omega_n pi-0.3, it might converge to 2pi. This is very relevant, as having twice repeated solutions screws over some physical quantities I want to calculate from there.

Also, the solution and eigenvalue guesses are very good. Their either calculated by hand or calculated from previous solve_bvp solutions.

My questions is:

  1. What is the integration method solve_bvp uses? Is there a method that might work better here?

I can't find this in solve_bvp documentation, and I fear I can't change it.

  1. What are things that usually affect the convergence of the parameters? Increasing max_nodes changes the error from 'repeated eigenvalue'/'solution did not converge' to MemoryError. Decreasing tol to have better accuracy doesn't seem to help much either.

I can gladly post the MWE, however I do not find it necessary/useful since:

  1. It is about 400-500 lines long,

  2. I am asking about the general behaviour of solve_bvp

Solution

  • The problem was in the numerical method used. The function whose roots I want to find is basically y(omega)(1), i.e. find omega such that the boundary conditions are met.

    When studying this function, it presented an extremum fairly close to the root. When solving this problem using Newton's method i.e. approximating the function to a straight line using the derivative of the function, if a guess omega came too close to the extremum, it would just shoot out very far away from the actual root. I fixed it using the bisection method, since I know approximately the shape of the function y(omega)(1).

    And no, I didn't need specifics of the functions I was using because explaining these choices and the code I am using is what I am writing my thesis for, and I don't plan on writing a second thesis for a SO question.