I am trying to implement a shooting method for a Sturm-Liouville problem on the interval [0,1] given by$$\frac{d}{dx}(P(x,\lambda^2)\frac{dy}{dx}) - Q(x,\lambda^2)y = 0$$with$$P(x,\lambda^2) = (1+sx)\lambda^2, \\Q(x,\lambda^2) = k^2((1+sx)\lambda^2 + sg),$$and parameters $k^2, s, g$. After providing a guess for the eigenvalue $\lambda^2$ the code needs to find the optimal $\lambda^2$ given the knowledge of Sturm's oscillation theorem and the shooting method. By chance I figured out a pretty good set of parameters. When I evaluate the ODE with my Runge-Kutta solver using the parameters and the boundary values $y_1(0) = 0$, $y_1'(0) = y_2(0)/P(0,\lambda^2) = 1$ it gives a good trajectory, but not the best one. The user is then asked to provide another guess for the eigenvalue and the system is again evaluated. When still not hit correctly in a user defined tolerance, a new eigenvalue is calculated using the secant method. This is the part where everything seems to go wrong. The optimal eigenvalue according to the secant method is of the order 1e-15, which is way to small, and occurs immediately after one iteration. As a consequence when I plot the corresponding eigenvector it is also small. It should be similar to the plot generated by the Runge-Kutta procedure, but with a better hit to $y(1) = 0$. I am aware that the secant method can converge pretty fast but why to such small values?  I have tried to spot a possible error in the past 2 days without success. I would be happy if someone could explain to me why this happens and moreover take his/her time to look into the, rather long, code.   import sysimport timeimport numpy as npimport matplotlib.pyplot as pltdef f(u, x, k, s, g, lam):    Define the ODE of the problem     u = (y1, y2)       f = f(y2/P, Q*y1)    -----------------------        du0/dx = f0 = u1/P(x,ww)        du1/dx = f1 = Q(x,ww)*u0    with u0(0) = 0.0 and u0'(0) = 1 = u1(0)/P(0,ww)          y1, y2 = u    return np.array([y2/P(p0,x,s,lam), Q(p0,x,k,s,g,lam) * y1])# ODE coefficientsdef p0(x, s):    '''Linear density profile    '''    return 1 + s*xdef P(p0, x, s, lam):    '''Defines the function P(x,w^2) in the ODE system'''    p0 = p0(x, s)    return p0*lamdef Q(p0, x, k, s, g, lam):    '''Defines the function Q(x,w^2) in the ODE system, note that it contains       the first derivative of the density profile    '''    p0 = p0(x, s)    return k*(p0*lam + s*g)My Runge-Kutta solver:def RungeKutta4(f, U0, N, tbegin, tend, **kwargs):        General function to solve generic first order OR first order coupled ODEs     of the form u' = f(u, t) with Runge-Kutta's 4th order Method    INPUT:    f           - function containing the first order ODE --> np.array when                       coupled ODEs    U0          - intial condition(s) for the ODE given as list    N           - # points to be generated on domain of interest = # iterations    tbegin      - begin point of integration domain    tend        - end point of integration domain    **kwargs    - (optional) dictionary containing constants needed in ODE                   function(s)                  NOTE: kwargs is just convention name, only aesterisc **                      needed followed by variable name    OUTPUT:    u       - 1D or 2D numpy array containing the derivative(s) corresponding to                   the values in the time array    t       - 1D array containing the time points         dt = tend/N #step size    u = np.zeros((N+1, len(U0)))    t = np.linspace(tbegin, tend, len(u))    u[0] = U0    for i in range(N):              k1 = dt*f(u[i], t[i], **kwargs)        k2 = dt*f(u[i] + 0.5*k1, t[i] + 0.5*dt, **kwargs)        k3 = dt*f(u[i] + 0.5*k2, t[i] + 0.5*dt, **kwargs)        k4 = dt*f(u[i] + k3, t[i] + dt, **kwargs)        u[i+1] = u[i] + (k1 + 2.0*(k2+k3) + k4)/6.0     return u, tMy shooting method:def ShootingMethod(f, bvs, fbvs, iv, err, **kwargs):        Function to solve boundary value problems by using the shooting method on a     first order (coupled) ODE of the form u' = f(u, t)     Making use of the Secant Method it possible that the algorithm not     converges, in order to overcome this we impose a maximum number of     iterations    INPUT:    f           - function containing the first order ODE --> np.array when                       coupled ODEs     bvs         - tuple containing the two boundaries a & b in which we                   integrate such that u(a) = X, u(b) = Y    fbvs        - tuple containing the function values X & Y at the boundaries    iv          - first guess for u'(0) at begin boundary such that u'(a) = iv    err         - accuracy we desire to get at the end boundary for accepting a                   solution    **kwargs    - (optional) dictionary containing constants needed in ODE                   function(s)    OUTPUT:    u       - 1D or 2D numpy array containing the derivative(s) corresponding to                   the values in the time array which satisfy boundary value problem    t       - 1D array containing the time points          #  maximum number of iterations before termination    Nmax = 20    # unpack the tuples of bvs and fbvs    (bv1, bv2) = bvs    (fbv1, fbv2) = fbvs    iv1 =  iv    # integrate the ODE to the other boundary using the first guess    u, t = RungeKutta4(f, U0=[fbv1,iv1], N=1000, tbegin=bv1, tend=bv2,                                 **kwargs)    # in order to look if the boundary is hit or missed we need last value    u1last = u[-1,0]    # determine the zero of function at boundary, which also gives error between     # numerical and exact solution at boundary    fiv1 = u1last-fbv2    print('\nTrial %d: iv1 = %10.8f, error = %10.8f' % (0, iv1, abs(fiv1)))    time.sleep(0.5)    # the probability that the first guess gives the solution within pre-chosen     # accuracy is small, so we do at least a second iteration with new guess    # we start with the second guess and follow the same procedure as earlier,     # if the second guess is not within the pre-chosen accuracy then we apply the     # Secant Method to refine our guess and start over again    iv2 = float(input('Please provide a new value for the guess: '))    print()    # it is forbidden that the second guess has same value as first and program     # terminates    if iv2==iv1:        raise ValueError('Second guess value cannot be the same as the first!')    Niter = 1       while Niter<Nmax:        u, t = RungeKutta4(f, U0=[fbv1,iv2], N=1000, tbegin=bv1, tend=bv2,                                   **kwargs)             u2last = u[-1,0]        fiv2 = u2last-fbv2        print('Trial %d: iv%d = %10.8f, error = %10.8f' % (Niter, Niter+1, iv2,                                                            abs(fiv2)))        time.sleep(0.5)            if abs(fiv2)<err:              print('\nSolution found!')            print('Exiting shooting algorithm to go plotting...')            time.sleep(3.0)            break        else:            # apply the Secant Method given our two failed guesses iv1, iv2 for             # searching a better root inbetween them            # our new point uv3 will be determined by the straight line through              # the points (iv1, f(uv1)) and (iv2, f(iv2)) with f(iv1) = fiv1 &             # f(iv2) = fiv2            # the calculation will be done now with iv2 and iv3 as initial values            # in order to accomplish this, iv1 will be overwritten by iv2 and iv2            # will be overwritten with our newly found root iv3            # of course f(iv1) needs also to be overwritten by f(iv2)              iv3 = iv2 - fiv2 * (iv2 - iv1)/(fiv2 - fiv1)            iv1 = iv2            iv2 = iv3            fiv1 = fiv2        Niter += 1    if Niter==Nmax:        print('\nNo convergence at the boundary after %d iterations'                  %Nmax)          print('Error is still: %.8f' %(abs(fiv2)))        print('Aborting program...')        time.sleep(3.0)        sys.exit()    return u, tThe main program:#k = k^2 and lam = lam^2k, s, g, lam = 2.0, 1.0, -9.8, 2.0# put parameters in dictionary to feed into functionsvars = ['k', 's', 'g', 'lam']values = [k, s, g, lam]params = dict(zip(vars, values))   # define the boundaries, function values at boundary, initial value and # accuracy where the solution has to be inbvs = (0.0, 1.0)fbvs = (0.0, 0.0)iv = P(p0, 0, s, lam)err = 1e-4# integrate with RK4, only for checking behaviouru, x = RungeKutta4(f, U0=[0.0,iv], N=1000, tbegin=0.0, tend=1.0,                          **params)# unpack the derivativesy1 = u[:,0]y2 = u[:,1]# plotting from Runge-Kutta                      plt.plot(x, y1)#plt.plot(x, y2)plt.show()  # apply shooting algorithm to do real problemu, x = ShootingMethod(f, bvs, fbvs, iv, err, **params)# unpack the derivativesy1shoot = u[:,0]y2shoot = u[:,1]# plotting from Shootingplt.plot(x, y1shoot) #plt.plot(x, y2shoot)plt.show()
Very small eigenvalues using shooting method on Sturm-Liouville problem
numerical analysis;python;computational physics