I need some help to minimize a numerical integration. It's about a classical problem in physics (hydrogen atom). It can be solved analytically but I need to solve it numerically in Python.

We have an electron's wave function $\psi(r,a)=e^{-r/a}$ and a Hamiltonian (in Ry units) $H=-\nabla^2-\frac{2}{r}$. The expected value of Energy is given by

$$\langle E(a)\rangle =\frac{\iiint \psi^{*}(r,a) H \psi(r,a) r^2 \sin(\theta) \mathrm dr \mathrm d\theta \mathrm d\phi}{\iiint |\psi(r,a)|^2 r^2 \sin(\theta) \mathrm dr \mathrm d\theta \mathrm d\phi}$$.

The integral it's over all space ($0\leq r \leq \infty$ ; $0\leq \theta \leq \pi$; $0\leq \phi \leq 2\pi$ )

Doing $\langle E \rangle$ numerically, how to find (also numerically) the value of $a$ which minimizes the value of the integral?

Edit 1- I know the value of '$a$' that minimizes the value of $<E(a)>$.

$a=1$ and $<E(1)>=-1$. But somehow I need to 'discovery' the value of $'a'$ that minimizes $<E(a)>$ numerically.

The code:

import numpy as np
from scipy.integrate import tplquad
import sympy as sp

# deffining the symbols
r, theta, phi, a = sp.symbols('r theta phi a')

# Deffining the function
f = sp.exp(-r/a)

# laplacian
laplacian = sp.diff(r**2 * sp.diff(f, r), r) / r**2 + \
            sp.diff(sp.sin(theta) * sp.diff(f, theta), theta) / 
(r**2 * sp.sin(theta)) + \
            sp.diff(f, phi, phi) / (r**2 * sp.sin(theta)**2)

# changing the functions into a lambda functions
laplacian_numeric = sp.lambdify((r, theta, phi, a), laplacian, 

f_numeric = sp.lambdify((r, theta, phi, a), f, 'numpy')

# deffining the kernel of integrations
def integrand(r, theta, phi, a):
    jac = r**2 * sp.sin(theta)
    return f_numeric(r, theta, phi, a) * laplacian_numeric(r, 
theta, phi, a)*jac

def integrand2(r,theta,phi,a):
    jac=r**2 * sp.sin(theta)
    return 2*jac*f_numeric(r, theta, phi, a)**2 /r

def integrand3(r,theta,phi,a):
    jac=r**2 * sp.sin(theta)
    return jac*f_numeric(r, theta, phi, a)**2

# integrations limits
r_limits = (0, np.inf)
theta_limits = (0, np.pi)
phi_limits = (0, 2*np.pi)

# value of the parameter a
a_value = 1

# numerical integration
h1, _ = tplquad(integrand, *phi_limits, *theta_limits,*r_limits, args=(a_value,))

h2,_ = tplquad(integrand2,*phi_limits, *theta_limits,*r_limits, args=(a_value,))

norm,_ = tplquad(integrand3,*phi_limits, *theta_limits,*r_limits, args=(a_value,))

print('result of numerical integration', (-h1-h2)/norm)
  • 4
    $\begingroup$ What have you tried? Do you have code that can compute the integrals correctly? If so, you can use one of the many optimization routines in, e.g., SciPy, such as scipy.optimize.minimize_scalar $\endgroup$
    – whpowell96
    Jan 22 at 17:53
  • 2
    $\begingroup$ If you want to minimize over only a single parameter (namely, $a$), why not just plot the quantity as a function of $a$? That will already give you intuition for what value of $a$ is approximately a minimum. $\endgroup$ Jan 22 at 18:34
  • $\begingroup$ I know the exact value of 'a' in this case. a=1. But I need to do this numerically. $\endgroup$ Jan 22 at 19:20
  • $\begingroup$ I'm doing a simple code now to explain better this problem. When I done, I will share with you. One moment, pls $\endgroup$ Jan 22 at 19:22


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.