"""
harmonic_oscillator.py
"""
#! /usr/bin/python
# -*- coding: utf-8 -*-
doc = """
# Simple Harmonic Oscillator (SHO)
# English :
# interest of this program : numerical experiment of schrodinger equation
# static solution of Harmonic Oscillator
# e > v : sinus
# e < v : exponential decreasing first but finally diverging.
# we can detect the energy of transition
# between positive and negative divergence.
# Français :
# résolution numérique de l'équation de Schroedinger
# recherche des états propres (e, phi)
# Solution indépendante du temps de l'oscillateur harmonique
# La partie e > v : oscille
# La partie e < v : décroît exponentiellement mais finit toujours par diverger
# l'énergie propre est comprise entre 2 énergies voisines qui divergent en sens contraires
"""
import math
import matplotlib.pyplot
def data( ) :
# usage : (k, f0, df0, dx, xmin, xmax, ymin, ymax) = data()
# point de départ de l'intégration numérique
k = 1 # potentiel : v = k x²
# solution paire : f0 > 0 et df0 = df0/dx = 0
f0 = 1
df0 = 0
dx = 0.01 # pas d'espace
# bornes du graphique
(xmin, xmax, ymin, ymax) = (0, 5, -3, 3)
return (k, f0, df0, dx, xmin, xmax, ymin, ymax)
def normalization_real(phi) :
# normalisation de la fonction phi (sous forme de tableau)
n = len(phi)
phi_integral = 0
for i in range(n) :
phi_integral += phi[i] * phi[i]
phi_integral = math.sqrt(phi_integral)
for i in range(n) :
phi[i] /= phi_integral
phi[i] /= phi_integral
return
def schroedinger_step(x0, fx0, dfx0, k, e, dx) :
"""
# intégration d'un pas d'espace de l'équation de Schroedinger
# i hbar part f / part t = - i hbar part^2 f / part x^2 + (V - E) f
# on pose toutes les constantes = 1
# on résoud : f''(x) = (V - E) f(x)
# avec V(x) = k x^2
# schéma numérique à "vitesse" décalée
# comme : f''(x) = [ f'(x+dx) - f'(x) ] / dx
# f'(x+dx/2) = f'(x) + (dx/2) f''(x)
# comme : f'(x+dx/2) = [ f(x+dx) - f(x) ] / dx
# f(x+dx) = f(x) + dx f'(x+dx/2)
"""
# print("x0, fx0, dfx0, k, e, dx =", (x0, fx0, dfx0, k, e, dx) )
x1 = x0 + dx
vx = k * x0**2
d2fx0 = (vx - e) * fx0
dfx1 = dfx0 + dx * d2fx0
fx1 = fx0 + dx * 0.5 * (dfx0 + dfx1)
return (x1, fx1, dfx1)
def plot_phi(x, phi, divergence, e) :
# displays phi(x)
matplotlib.pyplot.plot(x, phi)
if divergence > 0 :
matplotlib.pyplot.text(x[-2], phi[-2]-0.2, f"${e}$")
else :
matplotlib.pyplot.text(x[-2], phi[-2]-0.2, f"${e}$")
return
def plot_potential( x, k ) :
n = len(x)
v = [ None for i in range(n) ]
for i in range(n) :
v[i] = k * x[i]**2
matplotlib.pyplot.plot(x, v, color="black")
matplotlib.pyplot.text(1, 2.5, r"$V=kx^2$")
return
def schroedinger_real(e, n) :
"""
# la fonction phi(x) finit toujours par diverger
# quand elle est presque nulle, très en dessous du potentiel
# (exponentielle divergente)
"""
# state at x=0
(k, f0, df0, dx, xmin, xmax, ymin, ymax) = data()
x0 = 0
divergence = 0
# print("énergie testée : e=", e)
(x, phi, dphi) = ([0], [f0], [df0])
for i in range(n) :
(x1, f1, df1) = schroedinger_step(x0, f0, df0, k, e, dx)
# storage of i+1 for next schroedinger_step
x.append( x1 )
phi.append( f1 )
dphi.append( df1 )
# bascule :
(x0, f0, df0) = (x1, f1, df1)
if abs(f1) > ymax :
print(f" e = {e} : divergence de phi pour x={x1}")
if f1 > 0 : divergence = 1
else : divergence = -1 # fx1 < 0
break # for i
return (x, phi, dphi, divergence)
def schroedinger_dichotomie( ):
"""
# Dans l'équation de Schroedinger pour l'oscillateur harmonique,
# on teste différentes énergies de l'intervalle [e_min, e_max] :
# la fonction phi diverge plus ou moins vite
# Pour encadrer la valeur propre, il faut des divergences de sens opposées
"""
# state at x=0
(k, f0, df0, dx, xmin, xmax, ymin, ymax) = data()
n = int(xmax / dx)
print("nombre de pas d'intégration : n=", n)
# static schroedinger :
matplotlib.pyplot.title("Harmonic Oscillator")
matplotlib.pyplot.text(2, 1.1*f0, r"$E$")
matplotlib.pyplot.text(2, 0.3*f0, r"$\phi_E(x)$")
matplotlib.pyplot.axhline(0, 0, 4.5, color="black")
divergence = 0
e_min = 0.95
e_max = 1.05
print("# Test des énergies dans l'intervalle [", e_min, ",", e_max, "]")
(x, phi, dphi, divergence_min) = schroedinger_real(e_min, n)
plot_phi(x, phi, divergence_min, e_min)
(x, phi, dphi, divergence_max) = schroedinger_real(e_max, n)
plot_phi(x, phi, divergence_max, e_max)
if divergence_min * divergence_max >= 0 :
print(" Mauvais encadrement : mêmes divergences")
print(f" e_min = {e_min} ; e_max = {e_max}")
return
e0 = e_min
for j in range(10) :
e = (e_min + e_max) / 2
divergence_precedente = divergence
(x, phi, dphi, divergence) = schroedinger_real(e, n)
plot_phi(x, phi, divergence, e)
# tracé de l'énergie (droite horizontale ) comparer au potentiel
matplotlib.pyplot.plot( [x[0], x[len(x)-1]], [e, e] )
if divergence == divergence_min :
e_min = e
elif divergence == divergence_max :
e_max = e
else : # divergence = 0
print(f" e = {e} : pas de divergence avant xmax={xmax}")
break
print(f" l'énergie e est dans [ {e_min}, {e_max} ]")
plot_potential( x, k )
# normalization_real(phi)
matplotlib.pyplot.axis( [xmin, xmax, ymin, ymax] )
matplotlib.pyplot.show( )
return
schroedinger_dichotomie( )
"""
"""