In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
from sympy.vector import CoordSys3D
%matplotlib inline


    

    




---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-b55eba4a9ff6> in <module>()
      1 import numpy as np
      2 import matplotlib.pyplot as plt
----> 3 from sympy.vector import CoordSys3D
      4 get_ipython().magic('matplotlib inline')

ImportError: cannot import name 'CoordSys3D'

In [ ]:

    

hbar = 1.
m_0 = 1.
e_0 = 1.
eps_0 = 1.


    

In [ ]:

    

from pylab import *
from scipy . optimize import minimize
from scipy . integrate import quad
from scipy . special import hermite
from math import factorial
def Psi_ho (x , l , n ):
    H = hermite ( n )
    x = x / l
    return H ( x )* exp ( - x **2/2)/ sqrt ( sqrt ( pi )* factorial ( n )*2** n * l )
def V (x , Lambda ):
    return - Lambda *( Lambda -1)/2/ cosh ( x )**2
Lambda =5.5
n_max =40
l =0.5 # parmeter for basis function set
H = zeros (( n_max , n_max ))
for i1 in range (0 , n_max ):
    print ( i1 )
    for i0 in range (0 , i1 +1):
        E_pot = lambda x : \
            Psi_ho (x , l , i0 )* V (x , Lambda )* Psi_ho (x , l , i1 )
        H [ i0 , i1 ]= H [ i1 , i0 ]= \
            quad ( E_pot , -inf , inf , epsabs =1e-03 , epsrel =1e-03)[0]
        if abs( i0 - i1 )==2:
            n = min( i0 , i1 )
            H [ i0 , i1 ]+= -1/2*( sqrt (( n +1)*( n +2))/(2* l **2))
            H [ i1 , i0 ]+= -1/2*( sqrt (( n +1)*( n +2))/(2* l **2))
        if i0 == i1 :
            n = i0
            H [ i0 , i1 ]+= -1/2*( -(2* n +1)/(2* l **2))
E , psi = eigh ( H )
N =256
x = linspace ( -8 , 8 , N )
psi_0 = zeros ( N )
psi_1 = zeros ( N )
psi_2 = zeros ( N )
for i in range (0 , n_max ):
    psi_0 += psi [i , 0]* Psi_ho (x , l , i )
    psi_1 += psi [i , 1]* Psi_ho (x , l , i )
    psi_2 += psi [i , 2]* Psi_ho (x , l , i )
    if psi_0 [ N // 2] <0:
        psi_0 *= -1
    if psi_1 [ N // 2] <0:
        psi_1 *= -1
    if psi_2 [ N // 2] <0:
        psi_2 *= -1
plot (x , psi_0 , label = r'$ \ Psi_0 ( x ) $')
plot (x , psi_1 , label = r'$ \ Psi_1 ( x ) $')
plot (x , psi_2 , label = r'$ \ Psi_2 ( x ) $')
xlabel ( r'$x$ ')
ylabel ( r'$ \ Psi_n ( x ) $')
legend ()
from tight_layout ()


    

    




0
1






    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: overflow encountered in double_scalars
  # This is added back by InteractiveShellApp.init_path()
/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:11: RuntimeWarning: overflow encountered in cosh
  # This is added back by InteractiveShellApp.init_path()






    




2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

In [13]:

    

from sympy import *
init_printing()

x, y, z, kappa, R = symbols('x y z kappa, R')


    

In [18]:

    

psi_minus = exp(- kappa * sqrt(x ** 2 + y ** 2 + (z  - R) ** 2))
psi_plus = exp(- kappa * sqrt(x ** 2 + y ** 2 + (z  + R) ** 2))


    

In [20]:

    

psi = psi_plus  + psi_minus


    

In [29]:

    

lapl = simplify(diff(diff(psi, x), x) + diff(diff(psi, y), y) + diff(diff(psi, z), z))


    

In [27]:

    

S = exp(-kappa * R) * (1 + kappa*R + kappa**2 *R**2 / 3)
prefact = sqrt(kappa ** 3 / pi) / sqrt(2 * (1+ S))


    

In [28]:

    

prefact


    

    
Out[28]:





$$\frac{\sqrt{\kappa^{3}}}{\sqrt{\pi} \sqrt{2 \left(\frac{R^{2} \kappa^{2}}{3} + R \kappa + 1\right) e^{- R \kappa} + 2}}$$

In [30]:

    

ingra = (psi * lapl)


    

In [32]:

    

simplify(ingra)


    

    
Out[32]:





$$\frac{\kappa^{3}}{\pi \left(2 \left(\frac{R^{2} \kappa^{2}}{3} + R \kappa + 1\right) e^{- R \kappa} + 2\right)} \left(e^{- \kappa \sqrt{x^{2} + y^{2} + \left(R + z\right)^{2}}} + e^{- \kappa \sqrt{x^{2} + y^{2} + \left(R - z\right)^{2}}}\right) \left(\kappa^{2} e^{- \kappa \sqrt{x^{2} + y^{2} + \left(R + z\right)^{2}}} + \kappa^{2} e^{- \kappa \sqrt{x^{2} + y^{2} + \left(R - z\right)^{2}}} - \frac{2 \kappa e^{- \kappa \sqrt{x^{2} + y^{2} + \left(R + z\right)^{2}}}}{\sqrt{x^{2} + y^{2} + \left(R + z\right)^{2}}} - \frac{2 \kappa e^{- \kappa \sqrt{x^{2} + y^{2} + \left(R - z\right)^{2}}}}{\sqrt{x^{2} + y^{2} + \left(R - z\right)^{2}}}\right)$$

In [35]:

    

integrate(ingra, (x, -oo, +oo))


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-35-ae9b8df1c9e7> in <module>()
----> 1 integrate(ingra, (x, -oo, +oo))

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/integrals.py in integrate(*args, **kwargs)
   1278     if isinstance(integral, Integral):
   1279         return integral.doit(deep=False, meijerg=meijerg, conds=conds,
-> 1280                              risch=risch, manual=manual)
   1281     else:
   1282         return integral

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/integrals.py in doit(self, **hints)
    484                     function, xab[0],
    485                     meijerg=meijerg1, risch=risch, manual=manual,
--> 486                     conds=conds)
    487                 if antideriv is None and meijerg1 is True:
    488                     ret = try_meijerg(function, xab)

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/integrals.py in _eval_integral(self, f, x, meijerg, risch, manual, conds)
    885                 try:
    886                     if conds == 'piecewise':
--> 887                         h = heurisch_wrapper(g, x, hints=[])
    888                     else:
    889                         h = heurisch(g, x, hints=[])

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in heurisch_wrapper(f, x, rewrite, hints, mappings, retries, degree_offset, unnecessary_permutations)
    128 
    129     res = heurisch(f, x, rewrite, hints, mappings, retries, degree_offset,
--> 130                    unnecessary_permutations)
    131     if not isinstance(res, Basic):
    132         return res

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in heurisch(f, x, rewrite, hints, mappings, retries, degree_offset, unnecessary_permutations)
    660             solution = _integrate()
    661     else:
--> 662         solution = _integrate()
    663 
    664     if solution is not None:

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in _integrate(field)
    635 
    636         try:
--> 637             find_non_syms(raw_numer)
    638         except PolynomialError:
    639             return None

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:
--> 630                 list(map(find_non_syms, expr.args))
    631             else:
    632                 # TODO: Non-polynomial expression. This should have been

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/integrals/heurisch.py in find_non_syms(expr)
    625             elif expr in syms:
    626                 pass # ignore variables
--> 627             elif not expr.has(*syms):
    628                 non_syms.add(expr)
    629             elif expr.is_Add or expr.is_Mul or expr.is_Pow:

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/basic.py in has(self, *patterns)
   1134 
   1135         """
-> 1136         return any(self._has(pattern) for pattern in patterns)
   1137 
   1138     def _has(self, pattern):

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/basic.py in <genexpr>(.0)
   1134 
   1135         """
-> 1136         return any(self._has(pattern) for pattern in patterns)
   1137 
   1138     def _has(self, pattern):

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/basic.py in _has(self, pattern)
   1150         try:
   1151             match = pattern._has_matcher()
-> 1152             return any(match(arg) for arg in preorder_traversal(self))
   1153         except AttributeError:
   1154             return any(arg == pattern for arg in preorder_traversal(self))

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/basic.py in <genexpr>(.0)
   1150         try:
   1151             match = pattern._has_matcher()
-> 1152             return any(match(arg) for arg in preorder_traversal(self))
   1153         except AttributeError:
   1154             return any(arg == pattern for arg in preorder_traversal(self))

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/basic.py in __eq__(self, other)
    326                 return False
    327 
--> 328         return self._hashable_content() == other._hashable_content()
    329 
    330     def __ne__(self, other):

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/symbol.py in _hashable_content(self)
    219 
    220     def _hashable_content(self):
--> 221         return Symbol._hashable_content(self) + (self.dummy_index,)
    222 
    223 

/home/lorenzo/anaconda3/lib/python3.6/site-packages/sympy/core/symbol.py in _hashable_content(self)
    135         return {'_assumptions': self._assumptions}
    136 
--> 137     def _hashable_content(self):
    138         # Note: user-specified assumptions not hashed, just derived ones
    139         return (self.name,) + tuple(sorted(self.assumptions0.items()))

KeyboardInterrupt: 

In [ ]: