```python
from sympy import init_session
init_session()
```

    IPython console for SymPy 1.5.1 (Python 3.6.9-64-bit) (ground types: gmpy)
    
    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()
    
    Documentation can be found at https://docs.sympy.org/1.5.1/
    



```python
#init_printing?
```

### Parameters and two Gaussians


```python
a, b, c, a1, a2 = symbols('a b c a1 a2', positive=True, real=True)
```


```python
u1=x**2*exp(-a1*x**2)
u2=x**2*exp(-a2*x**2)
u1, u2
```

<br>

### Theory  (Fingers crossed)
We don't solve the radial equation for $R(r)$, but rather for $u(r) = r \cdot R(r)$.
So, if $R$ is a s-function, $u$ is a p-function, and if $R$ is a p-function, $u$ is a d-function, etc. <br>

$S = \int r^2dr R^2 = \int dr\, u^2$<br>
$T = \int r^2dr R_1 \hat{T} R_2 = \int dr\, u_1\, \frac{d^2}{dr^2} \, u_2 $
<br>

### Normalization constant


```python
N=integrate(u1*u1, (x, 0, oo))
N
```


```python
1/sqrt(N)
```


```python
printing.sstrrepr(1/sqrt(N))
```




    '4*2**(3/4)*sqrt(3)*a1**(5/4)/(3*pi**(1/4))'



### Overlap integral S


```python
S=integrate(u1*u2, (x, 0, oo))
S
```


```python
S.simplify()
```


```python
printing.sstrrepr(S.simplify())
```




    '3*sqrt(pi)/(8*(a1 + a2)**(5/2))'



### Kinetic energy $T = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} = \frac{1}{2m}\left(\frac{\hbar}{i}\frac{d}{dx} \right)^2$


```python
d1=diff(u1,x)
d2=diff(u2,x)
d1, d2
```


```python
Tsym = 1/2 * integrate(d1*d2, (x, 0, oo))
#T=T.simplify()
#T=T.factor()
Tsym.factor()
Tsym.simplify()
```


```python
dd1=diff(u1,x,2)
dd2=diff(u2,x,2)
dd1, dd2
```


```python
T1 = -1/2 * integrate(u1*dd2, (x, 0, oo))
T1.factor()
```


```python
T2 = -1/2 * integrate(u2*dd1, (x, 0, oo))
T2.factor()
```


```python
2/11
```


```python
printing.sstrrepr(Tsym.simplify())
```




    'sqrt(pi)*(1.875*a1*a2 - 0.25*(a1 + a2)**2)/(a1 + a2)**(7/2)'



### Potential $V(x) = (ax^2 - b)e^{-cx^2}$


```python
v=(a*x**2-b)*exp(-c*x**2)
v
```


```python
V = integrate(u1*v*u2, (x, 0, oo))
V
```


```python
V.factor()
```


```python
printing.sstrrepr(V.factor())
```




    '3*sqrt(pi)*(5*a - 2*a1*b - 2*a2*b - 2*b*c)/(16*(a1 + a2 + c)**(7/2))'



<br>

Angular momentum potential: $V = \frac{l(l+1)}{2\mu r^2 }$; for $l=1$: $V = \frac{1}{\mu r^2 }$


```python
vL=1/x**2
vL
```


```python
VL = integrate(u1*vL*u2, (x, 0, oo))
VL.simplify()
```


```python
printing.sstrrepr(VL.simplify())
```




    'sqrt(pi)/(4*(a1 + a2)**(3/2))'




```python

```

b9c75f8b7d6b9f70d793de8ded0ca23f3ac7f6cf
ipynb
Jupyter Notebook
notebooks/Stab_3D/.ipynb_checkpoints/Analytical_GTO_P-checkpoint.ipynb
tsommerfeld/L2-methods_for_resonances
acba48bfede415afd99c89ff2859346e1eb4f96c
MIT
notebooks/Stab_3D/.ipynb_checkpoints/Analytical_GTO_P-checkpoint.ipynb
tsommerfeld/L2-methods_for_resonances
acba48bfede415afd99c89ff2859346e1eb4f96c
MIT
notebooks/Stab_3D/.ipynb_checkpoints/Analytical_GTO_P-checkpoint.ipynb
tsommerfeld/L2-methods_for_resonances
acba48bfede415afd99c89ff2859346e1eb4f96c
MIT
Qwen/Qwen-72B
1. YES
2. YES