pythonfsolve

Solving TISE using solve_ivp and then a root finder

I am solving the TISE for x between 0 and L. I am given the boundary conditions that at x=0, phi(0)=0 and dphi(0)/dx=1. I also know that at x=L, phi(L)=0. I have used solve_ivp and a suggested initial energy of 134.0 which creates a correct graph.

from scipy.constants import c, m_e, e
me_ev  = m_e*c**2/e          # mass m_e in [eV]

l_bohr = 5.2918e-11          # Bohr radius [m]
convert_length = 5067730.179 # convert [m] to [1/eV] in natural units
l_ev   = l_bohr * convert_length

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import root_scalar

# TISE
def eqn (x, y, energy):
    me_ev  = m_e*c**2/e
    f_0 = y[1]    # dphi/dx
    f_1 = -2 * me_ev * energy * y[0]   # d2phi/dx2, y[0] = phi
    return np.array([f_0, f_1])

# using solve_ivp and the initial coonditions phi(0)=0 and dphi(0)/dx=1 to find phi(L)
def solve(energy, func):
    #constants
    l_bohr = 5.2918e-11                  # Bohr radius [m]
    convert_length = 5067730.179         # convert [m] to [1/eV] in natural units
    l_ev   = l_bohr * convert_length

    # x goes form [0,L]
    x = np.linspace(0., l_ev, 1000)
    
    # initial conditions
    phi0 = 0.0                         
    dphi0dx = 1.0
    init = np.array([phi0, dphi0dx])
    
    answer = solve_ivp(func, (0., l_ev), init, t_eval=x, args=(energy,))
    phiL = answer.y[0][-1]
    return phiL

This is where the code starts to have problems

# calculating ground state energy:
# use fsolve to put different energies into solve(energy, func) such that phiL is minimized:
ground_energy = fsolve(solve, 134.0, args=(eqn,))
print(ground_energy)

When I do this I get the error: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part.

Does anybody know how to fix this?

The function solve returns phi(L) so that I can optimize the energy such that phi(L)=0. This is where I run into problems. I try use fsolve.

Solution

  • If you print out f_0 and f_1 in eqn() you will find that f_1 is an array. Tracing that back, this is because the argument energy to eqn() is also an array (though neither you nor I expected it to be).

    You can get your code to run by changing the setting of f_1 to

    f_1 = -2 * me_ev * energy[0] * y[0]

    Whether it gives your intended answer requires a bit more explanation from you.

    BTW: please include all your imports so that we can run your code with minimal hunting around. I had to add the line

    from scipy.optimize import fsolve