schrodinger_infinite_well_bisect.py

def deviation(psi_computed, psi_actual):
    return psi_computed - psi_actual

def find_energy_eigenstate(E_low, E_high, integration_params):
    psi_right = 0.0

    psi_low = integrate(E_low, **integration_params)
    psi_high = integrate(E_high, **integration_params)
    dev_low = deviation(psi_low[-1], psi_right)
    dev_high = deviation(psi_high[-1], psi_right)

    def helper(low, high, d_low, d_high):
        E_mid = (low + high) / 2.0
        psi_mid = integrate(E_mid, **integration_params)
        if (high - low) < 1e-8:
          return E_mid, psi_mid

        d_mid = deviation(psi_mid[-1], psi_right)
        if (d_low * d_mid) < 0.0:
            return helper(low, E_mid, d_low, d_mid)
        else:
            return helper(E_mid, high, d_mid, d_high)

    return helper(E_low, E_high, dev_low, dev_high)

  
E_approx, psi = find_energy_eigenstate(E_low=1.0, E_high=1.5, integration_params=infinite_well)