Modélisation d'un oscillateur harmonique quantique/SHM

Question

J'ai besoin d'aide pour comprendre pourquoi mon graphique de l'état fondamental pour b) semble erroné, voici la question dans son intégralité : (J'ai pensé que l'afficher en entier donnerait un contexte à la méthode que j'essaie d'utiliser)

(a) Considérons un puits de potentiel carré avec ()=0 entre deux murs infiniment hauts séparés par une distance égale au rayon de Bohr, c'est-à-dire pour tout x dans l'intervalle [0,] .

• Écrire une fonction solve(energy, func) qui prend le paramètre énergie et une fonction python . Cette fonction doit résoudre l'EDO de Schrödinger pour le cas décrit ci-dessus et ne renvoyer que la valeur finale de () à la frontière .

• Ecrivez un script utilisant la fonction solve(energy, func) pour calculer l'énergie de l'état fondamental d'un électron dans ce puits de potentiel en unités d'eV (c'est-à-dire diviser le résultat par la valeur de la charge élémentaire). Pour la condition initiale, voir le conseil technique ci-dessous. Pour le nombre de valeurs à résoudre (), la valeur 1000 est recommandée.

Le résultat de votre calcul devrait être un nombre compris entre 134 eV et 135 eV. L'un des tests évaluera votre fonction solve(energy, func) pour un puits de potentiel déformé.

(b) Considérons le potentiel harmonique ()= ₀ 2 / 2 où ₀ et =10 11 m sont des constantes. Limiter le domaine infini de la variable à l'intervalle [10,10] avec ₀ \=50 eV.

• L'oscillateur harmonique est connu pour avoir des valeurs propres d'énergie équidistantes. Vérifiez que cela est vrai, à la précision de votre calcul, en calculant l'état fondamental et les deux premiers états excités. (Indice : l'état fondamental a une énergie comprise entre 100 et 200 eV.)

• Afin de tester votre résultat, écrivez une fonction result() qui doit renvoyer la différence des valeurs propres calculées en eV sous la forme d'un nombre unique. Notez que le test avec le nombre attendu est caché et testera votre résultat avec une précision de ±1 eV.

• Fournissez un titre de tracé et des étiquettes d'axe appropriées, tracez les trois fonctions d'onde sur un seul canevas à l'aide de lignes codées en couleur et fournissez également des limites d'axe appropriées en x et en y pour rendre toutes les courbes bien visibles.

Conseil technique : Il ne s'agit pas d'un problème de valeur initiale pour l'EDO de Schrödinger, mais d'un problème de valeur limite ! Cela demande un effort supplémentaire, une méthode différente des exercices précédents sur les EDO et est donc un problème "inédit".

• Prenons une condition initiale simple pour les deux problèmes à ₀ \=0 ou ₀ \=10 , respectivement : ( ₀ )=0 et ( ₀ )/=1 . Utilisez ce point de départ pour résoudre l'EDO et faites varier l'énergie, , jusqu'à ce qu'une solution converge. La tâche consiste à évaluer la variation de la variable énergie jusqu'à ce que la deuxième condition limite (par exemple à L pour l'exercice (a)) soit satisfaite puisque la première condition limite est déjà satisfaite en raison du choix de la condition initiale.

• Cela nécessite une estimation initiale de l'intervalle d'énergie dans lequel la solution pourrait se trouver et une méthode de calcul pour la recherche de Root. Recherchez dans Scipy des méthodes de recherche de racines et sélectionnez-en une qui ne nécessite pas la connaissance de la dérivée. La recherche de racine est appropriée ici puisque la condition limite à satisfaire est ()=0.

Dans le cadre de la physique quantique, la condition aux limites pour les deux exercices est que ()=0 à chaque limite de potentiel. Ces conditions ne sont remplies qu'à des valeurs d'énergie spécifiques et discrètes, appelées valeurs propres de l'énergie, dont la plus petite est appelée énergie de l'état fondamental.

m_el   = 9.1094e-31      # mass of electron in [kg]
hbar   = 1.0546e-34      # Planck's constant over 2 pi [Js]
e_el   = 1.6022e-19      # electron charge in [C]
L_bohr = 5.2918e-11      # Bohr radius [m]
V0 = 50*e_el
a = 10**(-11)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
from scipy.integrate import odeint
from scipy import optimize

def eqn(y, x, energy):                       #part a)
    y0 = y[1]
    y1 = -2*m_el*energy*y[0]/hbar**2

    return np.array([y0,y1])

x = np.linspace(0,L_bohr,1000)
def solve(energy, func):
    p0 = 0
    dp0 = 1
    init = np.array([p0,dp0])
    ysolve = odeint(func, init, x, args=(energy,))
    return ysolve[-1,0]

def eigen(energy):

    return solve(energy, eqn)

root_ = optimize.toms748(eigen,134*e_el,135*e_el)
root = root_/e_el

print('Ground state infinite square well',root,'eV')

intervalb = np.linspace(-10*a,10*a,1000)                    #part b)

def heqn(y, x2, energy):
    f0 = y[1]
    f1 = (2.0 * m_el / hbar**2) * (V0 * (x2**2/a**2) - energy) * y[0]

    return np.array([f0,f1])

def solveh(energy, func):
    ph0 = 0
    dph = 1
    init = np.array([ph0,dph])
    ysolve = odeint(func, init, intervalb, args=(energy,))
    return ysolve

def boundary(energy):        #finding the boundary V=E to apply the b.c
    f = a*np.sqrt(energy/V0)
    index = np.argmin(np.abs(intervalb-f))

    return index

def eigen2(energy):

    return solveh(energy,heqn)[boundary(energy),0]

groundh_ = optimize.toms748(eigen2,100*e_el,200*e_el)
groundh = groundh_/e_el

print('Ground state of Harmonic Potential:', groundh, 'eV')
plt.suptitle('Harmonic Potential Well')
plt.xlabel('x (a [pm])')
plt.ylabel('Psi(x)')

groundsol = solveh(groundh_,heqn)[:,0]
plt.plot(intervalb/a, groundsol)

La forme du graphique est la suivante pour toutes les valeurs d'énergie comprises entre 100 eV et 200 eV. Je ne comprends pas où je me trompe. J'ai essayé de tester mon code autant que possible.