<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 ENERGY EIGENVALUES OF SPIKED HARMONIC <br> OSCILLATOR POTENTIALS USING AIMpy CODE
M. Karakoc
Presented at https
Step1:  Definitions
Variables
Step2:  $\lambda_0$ and $s_0$
Step3:  Case
Step4:  Initialize AIM solver
Step5:  Calculation necessary coefficients
Step6:  The Solution
Step7:  Case
Step8:  Case
Step9:  Case
Step10:  Eigenvalues and Eigenfucntions for $\alpha = 4, \gamma = 5, \lambda = 100$
<ASSISTANT_TASK:>
Python Code:
# Python program to use AIM tools
from asymptotic import *
Explanation: ENERGY EIGENVALUES OF SPIKED HARMONIC <br> OSCILLATOR POTENTIALS USING AIMpy CODE
M. Karakoc
Presented at https://technoscience2020.weebly.com/live-hall.html
Reference:
ASYMPTOTIC ITERATION METHOD FOR SINGULAR POTENTIALS
B. Champion, R. L. Hall, N. Saad, Int. J. Mod. Phys. A 2008.23:1405-1415.
https://www.worldscientific.com/doi/abs/10.1142/S0217751X08039852
Import AIM library
End of explanation
# symengine (symbolic) variables for lambda_0 and s_0 
En, r = se.symbols("En, r")
r0, m, λ, γ = se.symbols("r0, m, λ, γ")
Explanation: Definitions
Variables
End of explanation
# lambda_0 and s_0
l0 = 2*r - 2*se.sqrt(λ)/r**(m+1) - (1 + m)/r
s0 = se.expand(2 + m + 2*se.sqrt(λ)/r**m + ((2*γ+1-m)*(2*γ+1+m))/(4*r**2) - En)
Explanation: $\lambda_0$ and $s_0$
End of explanation
nr0 = o* 3
nλ  = o* 1/10
nm  = o* 1/2
nγ  = o* 2
nα  = 2* (nm + 1) 
print('nα =', nα)
# parameters of lambda_0 (pl0) and s_0 (ps0)
pl0 = {λ: nλ, m:nm}
ps0 = {λ: nλ, m:nm, γ:nγ}
Explanation: Case: $\alpha = 3, m = 1/2, \gamma = 2, \lambda = 0.1, r_0 = 3$
Numerical values for variables¶
End of explanation
%%time
# pass lambda_0, s_0 and variable values to aim class
singular = aim(l0, s0, pl0, ps0)
singular.display_parameters()
singular.display_l0s0(0)
singular.parameters(En, r, nr0, nmax=51, nstep=10, dprec=1000, tol=1e-101)
Explanation: Initialize AIM solver
End of explanation
%%time
# create coefficients for improved AIM
singular.c0()
singular.d0()
singular.cndn()
Explanation: Calculation necessary coefficients
End of explanation
%%time
singular.get_arb_roots(showRoots='+r', printFormat="{:25.20f}")
Explanation: The Solution
End of explanation
nr0 = o* 3
nλ  = o* 1/10
nm  = o* 1
nγ  = o* 0
nα  = 2* (nm + 1) 
print('nα =', nα)
# parameters of lambda_0 (pl0) and s_0 (ps0)
pl0 = {λ: nλ, m:nm}
ps0 = {λ: nλ, m:nm, γ:nγ}
# pass lambda_0, s_0 and variable values to aim class
singular = aim(l0, s0, pl0, ps0)
singular.display_parameters()
singular.display_l0s0(0)
singular.parameters(En, r, nr0, nmax=51, nstep=10, dprec=500, tol=1e-101)
# create coefficients for improved AIM
singular.c0()
singular.d0()
singular.cndn()
singular.get_arb_roots(showRoots='+r', printFormat="{:25.20f}")
Explanation: Case: $\alpha = 4, m = 1, \gamma = 0, \lambda = 0.1, r_0 = 3$
End of explanation
nm  = o* 1
nα  = 2* (nm + 1) 
print('nα =', nα)
nr0 = o* 3
for nλ in [o* 1000/10**gi for gi in range(7)]:
    print('\033[91m\033[1m{linea} λ ={fnλ:>8s} {linea}\033[0m'.format(linea = '='*30, fnλ = str(nλ)))
    for nγ in [o* gi for gi in range(6)]:
        print('{linea} γ ={fnγ:>2s} {linea}'.format(linea = '-'*30, fnγ = str(nγ)))
        
        # parameters of lambda_0 (pl0) and s_0 (ps0)
        pl0 = {λ: nλ, m:nm}
        ps0 = {λ: nλ, m:nm, γ:nγ}
        
        # pass lambda_0, s_0 and variable values to aim class
        singular = aim(l0, s0, pl0, ps0)
        #singular.display_parameters()
        #singular.display_l0s0(0)
        singular.parameters(En, r, nr0, nmax=51, nstep=10, dprec=500, tol=1e-101)
        
        # create coefficients for improved AIM
        singular.c0()
        singular.d0()
        singular.cndn()
        
        
        singular.get_arb_roots(showRoots='+r', printFormat="{:25.20f}")
Explanation: Case: $\alpha = 4, m = 1, \gamma = 0\ldots5$
$r_0 = 3$
$\lambda = 0.001, 0.01, 0.1 ,1 ,10, 100, 1000$
End of explanation
# symengine (symbolic) variables for lambda_0 and s_0 
En, r = se.symbols("En, r")
r0, α, λ, γ = se.symbols("r0, α, λ, γ")
# lambda_0 and s_0
l0 = 2*(r - (1 + γ)/r)
s0 = 2*γ + 3 + λ/r**α - En
nα  = o* 1
nr0 = o* 5
for nλ in [o* 1000/10**gi for gi in range(10)]:
    print('\033[91m\033[1m{linea} λ ={fnλ:>11s} {linea}\033[0m'.format(linea = '='*30, fnλ = str(nλ)))
    for nγ in [o* gi for gi in range(1)]:
        print('{linea} γ ={fnγ:>2s} {linea}'.format(linea = '-'*30, fnγ = str(nγ)))
        
        # parameters of lambda_0 (pl0) and s_0 (ps0)
        pl0 = {λ: nλ, γ:nγ}
        ps0 = {λ: nλ, γ:nγ, α:nα}
        
        # pass lambda_0, s_0 and variable values to aim class
        singular = aim(l0, s0, pl0, ps0)
        #singular.display_parameters()
        #singular.display_l0s0(0)
        singular.parameters(En, r, nr0, nmax=51, nstep=10, dprec=500, tol=1e-101)
        
        # create coefficients for improved AIM
        singular.c0()
        singular.d0()
        singular.cndn()
        
        
        singular.get_arb_roots(showRoots='+r', printFormat="{:25.20f}")
Explanation: Case: $\alpha = 1, \gamma = 0\ldots5$
$\lambda = 1000/10^n$ where $n=0\ldots9$
End of explanation
!pip3 install scipy
from scipy.integrate import odeint 
# https://perso.crans.org/besson/publis/notebooks/Runge-Kutta_methods_for_ODE_integration_in_Python.html
%matplotlib notebook    
import scipy.optimize as so
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
V = lambda r, α, β, γ, λ: r**β + γ*(γ+1)/r**2 + λ/r**α
# potential parameters
αβγλ = 4, 2, 5, 100
# Energy Eigenvalues
Estop = 34
Espectrum = [15.54392598171621432978,
             19.92653800987177727185,
             24.26619272363845153429,
             28.57312197406556227733,
             Estop]
# Draw a line for each eigenvalue
EigenLine = []
for Ei in Espectrum:
    Vr = lambda x: V(x, *αβγλ) - Ei
    rL, rR = so.newton(Vr, 1e-1), so.newton(Vr, 6)
    EigenLine += [[[rL, rR], [Ei, Ei]]]
    
# create colors for the lines
Eimin = .7 * Espectrum[0]/Espectrum[-1]
Eimax = .7
cm_subsection = np.linspace(Eimin, Eimax, len(Espectrum)) 
colors = [ cm.hot_r(x) for x in cm_subsection ]
Scale = 4.5
plt.figure(figsize=(Scale, Scale*1.618))
for i, Ei in enumerate(EigenLine[:-1]):
    plt.axhline(Espectrum[i], color=colors[i], linestyle=':')
    plt.plot(*Ei, color=colors[i])
    
    
for i, Ei in enumerate(Espectrum[0:-1]):
    def pend(y, t):
        return np.array([y[1], -(Ei - V(t, *αβγλ))*y[0]])
    
    y0L =  1e-1, 1e-1
    y0R =  (-1)**i * 1e-1, (-1)**i * 1e-1
    
    rL, rR = EigenLine[i][0]
    rM = (rL + rR)/2
    tL = np.linspace(1e-1, rM, 101)
    tR = np.linspace(9, rM, 101)
    
    solL = odeint(pend, y0L, tL)
    solR = odeint(pend, y0R, tR)
    
    WL = solL[:, 0]/abs(solL[:, 0][-1])
    WR = solR[:, 0]/abs(solR[:, 0][-1])
    Wmax = max(max(abs(WL)), max(abs(WR)))
    
    plt.plot(tL, Ei + 1.6*WL/Wmax, color=colors[i])
    plt.plot(tR, Ei + 1.6*WR/Wmax, color=colors[i])
    
rs = np.linspace(*EigenLine[-1][0], 500)    
plt.plot(rs, V(rs, *αβγλ), color='blue')
plt.xlabel('r [arbitrary]', fontsize=14)
plt.ylabel('$\psi(r)$, V(r), E  [all arbitrary]', fontsize=14)
#plt.grid(color='lightgray')
plt.grid(axis='x', linestyle='--', color='lightgray')
#plt.grid(linestyle='--', color='lightgray')
plt.xlim(EigenLine[-1][0])
plt.show()
Explanation: Eigenvalues and Eigenfucntions for $\alpha = 4, \gamma = 5, \lambda = 100$
End of explanation
<END_TASK>