{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Un problème en mécanique quantique numérique: Le potentiel de Lennard-Jones avec numpy et SciPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On va numériquement étudier l'équation de Schrödinger indpépendante du temps en une dimension avec un potentiel de Lennard-Jones. Ce potentiel décrit, en trois dimensions, les interactions entre des molecules de manière phénoménologique et approximative. Il est donné par $$V(x)=\\epsilon\\,\\left(\\left(\\frac{x_0}{x}\\right)^{12}-2\\,\\left(\\frac{x_0}{x}\\right)^6\\right)\\,.$$ On ne regarde que les $x$ positifs. Par inspection on déduit que\n", "* $V(x)\\,\\rightarrow\\,\\infty\\;$ lorsque $x\\,\\rightarrow\\,0$ et $V(x)\\,\\rightarrow\\,0\\;$ lorsque $x\\,\\rightarrow\\,\\infty$\n", "* il y a un minimum global à $x=x_0$ avec énergie $-\\epsilon$\n", "\n", "Traçons d'abord le potentiel pour plusieurs valeurs de $\\epsilon$ et $x_0=1$ entre $x=0.8$ et $x=2$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = 1.\n", "def V(x, epsilon):\n", " return epsilon * ((x0 / x)**12 - 2 * (x0 / x)**6)\n", "\n", "xpoints = np.linspace(.8, 2, 100)\n", "for epsilon in [1., 2., 10.]:\n", " plt.plot(xpoints, V(xpoints, epsilon), label=r'$\\epsilon$ = ' + str(epsilon))\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"V(x)\")\n", "plt.ylim(-5, 5)\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### La discrétisation de l'hamiltonien" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'équation de Schrödinger indépendante du temps est $$H\\psi(x)=E\\psi(x)\\,\\qquad\\qquad\\text{avec } H=\\left(-\\frac{\\hbar^2}{2m}\\frac{d^2}{dx^2}+V(x)\\right)\\,.$$ \n", "\n", "Notre premier objectif sera le calcul numérique des énergies $E$ des états liés, alors les valeurs propres négatives de l'opérateur $H$. Mais cet opérateur agît sur un espace vectoriel de dimension infinie (l'espace de Hilbert du problème). Pour le traitement numérique il faut *discrétiser* (prendre l'approximation d'un espace d'états de dimension finie), ce qui transforme l'opérateur différentiel linéaire $H$ en *matrice* $(N+1)\\times (N+1)$. Plus $N$ est grand, plus notre approximation sera précise.\n", "\n", "De plus, au lieu de considérer une fonction d'onde $\\psi(x)$ sur $\\mathbb{R}_+$ on va se limiter à un domaine de solutions fini. En fait on s'attend que $\\psi(x)\\,\\rightarrow\\,0$ pour $x$ suffisamment grand pour des modes normalisables, alors on résoudra le système pour $x\\in[0, x_{\\text{max}}]$ en choisissant $x_{\\text{max}}$ suffisamment grand. Plus $x_{\\text{max}}$ est grand, plus notre approximation sera précise.\n", "\n", "Mais l'espace de fonctions d'onde sur $[0, x_{\\text{max}}]$ est toujours de dimension infinie. Alors on regardera un ensemble fini de $N+1$ valeurs de fonction $\\psi(na)$ (avec $n=0\\ldots N$), où\n", "$$\n", "a = \\frac{x_{\\text{max}}}{N}\\,.\n", "$$\n", "\n", "Qu'est-ce que la discrétisation de l'opérateur différentiel $\\frac{d^2}{dx^2}$? On rappelle le développement limité d'une fonction $f(x)$,\n", "$$\\begin{split}f(x+a)=&\\;f(x)+a\\,f'(x)+\\frac{1}{2} a^2\\,f''(x)+\\frac{1}{6}a^3\\,f'''(x)+{\\cal O}(a^4),\\\\ f(x-a)=&\\;f(x)-a\\,f'(x)+\\frac{1}{2} a^2\\,f''(x)-\\frac{1}{6}a^3\\,f'''(x)+{\\cal O}(a^4)\\,.\\end{split}$$\n", "En prenant la somme de ces deux équations, les puissances impaires de $a$ se suppriment et on arrive à\n", "$$f(x+a)+f(x-a)=2\\,f(x)+a^2\\,f''(x)+{\\cal O}(a^4)\\;\\Rightarrow\\;f''(x)=\\frac{f(x+a)+f(x-a)-2\\,f(x)}{a^2}+{\\cal O}(a^2)\\,.$$\n", "On conclut que, pour $a$ suffisamment petit (c.à.d. $N$ suffisamment grand), on peut approximer\n", "$$\n", "\\frac{d^2}{dx^2}\\psi(x)\\approx \\frac{\\psi(x+a)+\\psi(x-a)-2\\,\\psi(x)}{a^2}\\,.\n", "$$\n", "Si on met tous les valeurs de la fonction d'onde $\\psi(na)$ dans un vecteur de colonne $\\Psi$,\n", "$$\n", "\\Psi=\\left(\\begin{array}{c}\\psi(0) \\\\ \\psi(a) \\\\ \\psi(2a) \\\\ \\vdots\\\\ \\psi(Na)\\end{array}\\right)\n", "$$\n", "alors l'action de l'opérateur $\\frac{d^2}{dx^2}$ devient une multiplication matricielle:\n", "$$\n", "\\frac{d^2}{dx^2}\\psi\\approx \\frac{1}{a^2}\\left(\\begin{array}{ccccc}\\ddots &&&&\\\\ 1 & -2 & 1 && \\\\ &1&-2&1& \\\\ &&1&-2&1 \\\\ &&&&\\ddots\\end{array}\\right)\\Psi\\,.\n", "$$\n", "D'où la partie cinétique du hamiltonien peut s'écrire approximativement \n", "$$\n", "\\frac{\\hbar^2}{2m}\\frac{d^2}{dx^2}\\approx \\frac{\\hbar^2}{2ma^2}\\,M\\,,\n", "%\\left(\\begin{array}{ccccc}\\ddots &&&&\\\\ 1 & -2 & 1 && \\\\ &1&-2&1& \\\\ &&1&-2&1 \\\\ &&&&\\ddots\\end{array}\\right)\\,.\n", "$$\n", "$M$ étant une matrice $(N+1)\\times(N+1)$ avec des $-2$ sur la diagonale principale, des $1$ aux côtés de la diagonale principale, et des $0$ ailleurs. Construisons cette matrice pour $N=10$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# numpy.identity(n) va construire la matrice d'identité n x n\n", "# numpy.diag(v, k=0) va construire une matrice avec le vecteur v sur la k-ème diagonale, avec k=0 la principale.\n", "def matriceM(N):\n", " M = -2 * np.identity(N + 1) # la matrice -2 I de taille (N+1)x(N+1)\n", " ones = np.ones(N) # le tableau [1, 1,... 1] de longueur N\n", " M += np.diag(ones, -1)\n", " M += np.diag(ones, 1)\n", " return M\n", "\n", "print(matriceM(10))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La représentation discrète de la partie $V(x)$ du hamiltonien est plus simple: C'est une matrice diagonale avec les entrées diagonales $V(0),\\,V(a),\\,V(2a) \\ldots V(Na)$. Le seul problème est que $V(0)$ diverge, alors on posera $V(0)=$ \"infini\" (= un très grand nombre) à la main. \n", "\n", "Nous sommes maintenant en mesure de construire le hamiltonien complet. En unités où $\\hbar=1$ et $m=1$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def hamiltonien(epsilon, xmax, N):\n", " a = xmax / N\n", " xn = np.linspace(0, xmax, N + 1)\n", " with np.errstate(divide='ignore', invalid='ignore'): # sinon Python va se plaindre car pour V(0) on divise par 0\n", " Vxn = V(xn, epsilon) \n", " Vxn[0] = 1.E8 # \"infini\" = une constante >> la plus grande autre échelle du problème\n", "\n", " H = -1 / (2 * a**2) * matriceM(N) + np.diag(Vxn)\n", " return H\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Détermier l'état fondamental d'un système avec un seul état lié" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour trouver les valeurs propres de cette matrice, on se servira de la méthode **scipy.linalg.eig(a)**. Cette méthode prend une matrice carrée comme argument et renvoie un tableau 1-dimensionnel qui contient les valeurs propres, ainsi qu'un tableau 2-dimensionnel dont les colonnes contiennent les composantes des vecteurs propres normalisés. Pour les valeurs des paramètres, nous choisirons $\\epsilon=10$, $x_{\\text{max}}=5$ et $N=100$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.linalg\n", "H = hamiltonien(10, 5, 100)\n", "E, vecteurspropres = scipy.linalg.eig(H) # fonction qui renvoie les valeurs propres et vecteurs propres normalisés\n", " # (ces derniers sont les colonnes de la matrice 'vecteurspropres')\n", "print(\"Énergies:\")\n", "print(np.sort(E))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il ne faut pas trop regarder les valeurs propres positives. En fait pour les états de diffusion on s'attend un spectre continu, la quantification des valeurs propres ici est un artefact numérique.\n", "\n", "Par contre il semble qu'il y a un état lié avec énergie $E=-0.13869\\ldots$ en nos unités. Le vecteur propre correspondant est la fonction d'onde $\\psi_0(x)$ de l'état fondamental, stocké dans la colonne du même indice de la matrice **vecteurspropres**. On trace $\\psi_0$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n0 = np.argmin(E) # donne l'indice de l'élément minimal de E = l'énergie de l'état fondamental\n", "psi0 = vecteurspropres[:, n0] # la fonction d'onde de l'état fondamental est la colonne qui correspond à cet indice\n", "x = np.linspace(0, 5, 101)\n", "plt.plot(x, np.abs(psi0))\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\psi_0$\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Est-ce qu'on peut faire confiance à ces résultats? Pour le savoir, on repète le calcul avec \n", "* un plus grand nombre $N$ de points\n", "* un intervalle $[0,x_{\\text{max}}]$ plus large\n", "* un plus grand nombre de points et un intervalle plus large.\n", "\n", "À noter que les vecteurs propres calculés par **scipy.linalg.eig** sont normalisés par la condition\n", "$$\n", "\\sum_{n=0}^N |\\psi(n\\,a)|^2=1\n", "$$\n", "tant que la bonne normalisation de la fonction d'onde nécessite\n", "$$\n", "1\\stackrel{!}{=} \\int |\\psi(x)|^2\\,dx\\approx a\\sum_{n=0}^N |\\psi(n\\,a)|^2\\,.\n", "$$\n", "Pour avoir des fonctions d'onde avec la bonne normalisation relative, il faut alors diviser par $\\sqrt{a}$.\n", "\n", "Regardons ce qui change si on augmente $N$ et $x_{\\text{max}}$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def etat_fond(epsilon, xmax, N): # renvoie l'énergie et la fonction d'onde de l'état fondamental\n", " H = hamiltonien(epsilon, xmax, N)\n", " E0, vecteurspropres = scipy.linalg.eig(H)\n", " n0 = np.argmin(E0)\n", " psi0 = vecteurspropres[:, n0]\n", " norm = np.sqrt(N / xmax)\n", " return E0, norm * np.abs(psi0)\n", "\n", "for xmax in [5, 10]:\n", " for N in [200, 500]:\n", " E, psi = etat_fond(10., xmax, N)\n", " print(\"Énergie de l'état fondamental pour xmax = \", xmax, \", N = \", N, \":\", E.min())\n", " x = np.linspace(0, xmax, N+1) \n", " plt.plot(x, psi, label = \"xmax = \" + str(xmax) + \", N = \" + str(N)) \n", "plt.xlim(right=5)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\psi_0$\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il semble que l'énergie est bien déterminée par notre procédure (avec une précision relative de $\\lesssim$ 1%). Les fonctions d'onde s'accordent bien également; en fait elles ne sont presque pas discernables par œil dans le graphique. Un peu plus précisément, on peut directement regarder la différence entre l'approximation la moins et la plus précise:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "E1, psi1 = etat_fond(10, 5, 200) # une approximation peu précise\n", "E2, psi2 = etat_fond(10, 10, 500) # une approximation plus précise\n", "x = np.linspace(0, 5, 51)\n", "plt.plot(x, psi1[::4] - psi2[:251:5])\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\Delta\\psi_0$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visiblement l'erreur est petite (de l'ordre de $10^{-3}$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plusieurs états liés" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On va maintenant augmenter la profondeur du potentiel pour trouver plus qu'un état lié:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "epsilon = 100.\n", "N = 500\n", "xmax = 10.\n", "H = hamiltonien(epsilon, xmax, N) \n", "norm = np.sqrt(N / xmax)\n", "E, vecteurspropres = scipy.linalg.eig(H)\n", "fonctionsdonde = [] # pour enregistrer les coefficients des vecteurs propres\n", "for e in E[E < 0]: # identifier les énergies des états liés = les valeurs propres négatives\n", " print(\"État lié trouvé avec énergie\", e)\n", " fonctionsdonde += [norm * vecteurspropres[:, np.where(E==e)].flatten()] \n", "x = np.linspace(0, 5, 251) # tracer les fonctions d'onde sur l'intervalle [0, 5] = la première moitié des points calculés\n", "for psi in fonctionsdonde: \n", " plt.plot(x, psi[:251])\n", "xx = np.linspace(0.8, 5, 300) # dans le même graphique, tracer aussi le potentiel entre x = 0.8 et x = 5\n", "plt.plot(xx, 0.1*V(xx, epsilon), 'r--', label=\"V(x)\") # (unités arbitraires, facteur 0.1 pour mieux voir la courbe)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\psi(x)$\")\n", "plt.ylim([-5, 3])\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On s'aperçoit de trois états liés avec les fonctions d'onde à $0$, $1$ et $2$ nœuds. \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparaison avec l'approximation harmonique" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Peut-on approximer le potentiel aux alentours de son minimum par un oscillateur harmonique? On note que, après un développement limité,\n", "$$\n", "V(x_0+x)=V(x_0)+x\\,V'(x_0)+\\frac{1}{2} x^2 V''(x_0)+{\\cal O}(x^3)\n", "$$\n", "où $V'(x_0)=0$ car il s'agît d'un minimum et $V(x_0)=-\\epsilon$ n'est qu'une constante. Le calcul de $V''(x_0)$ donne \n", "$$V''(x_0)=72\\,\\frac{\\epsilon}{x_0^2}\\,.$$\n", "Si les termes d'ordre $\\geq x^3$ sont négligéables, en comparant avec le potentiel habituel de l'oscillateur harmonique\n", "$$\n", "V_{\\text{OH}}=\\frac{m\\omega^2}{2}x^2\n", "$$\n", "on déduit qu'il faut identifier \n", "$$\\omega=\\frac{6}{x_0}\\,\\sqrt{\\frac{2\\,\\epsilon}{m}}\\,.$$\n", "Cependant, en juste regardant les courbes des deux potentiels on s'attend que l'approximation ne va pas être très bonne:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "omega = 6 / x0 * np.sqrt(2 * epsilon)\n", "plt.plot(xx, V(xx, epsilon), label=\"Potentiel de Lennard-Jones\")\n", "plt.plot(xx, omega**2 / 2 * (xx - x0)**2 - epsilon, label = \"Oscillateur harmonique\")\n", "plt.xlim(0.8, 2.0)\n", "plt.ylim(-100, 20)\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les énergies des états liés de l'oscillateur harmonique sont données par la formule bien connue\n", "$$\n", "E_n=\\hbar\\omega\\left(n+\\frac{1}{2}\\right)\\,,\\qquad n\\geq 0\n", "$$\n", "de laquelle il faut soustraire $\\epsilon$ pour arriver à la même normalisation du minimum que pour notre potentiel de Lennard-Jones. Alors:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Premiers niveaux d'énergie d'un oscillateur harmonique:\")\n", "for n in range(3):\n", " print(omega * (n + 0.5) - epsilon)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visiblement l'approximation est tout à fait assez mauxaise: l'énergie de l'état fondamental est réproduite à 10% près, mais les autres deux états liés sont complètement faux. \n", "\n", "On peut également comparer les fonctions d'onde. Les fonctions propres de l'oscillateur harmonique sont données par\n", "$$\n", "h_n(x)=N_n\\,\\exp\\left(-\\frac{m\\omega x^2}{2\\hbar}\\right)\\,H_n\\left(\\sqrt{\\frac{m\\omega}{\\hbar}}\\,x\\right)\n", "$$\n", "où $N_n=\\left(2^n\\,n!\\right)^{-1/2}\\left(\\frac{m\\omega}{\\pi\\hbar}\\right)^{1/4}$ est un facteur de normalisation et $H_n$ est le $n$-ème *polynôme d'Hermite*, défini dans la bibliothèque des fonctions mathématiques spéciales **scipy.special**. On définit alors" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.special\n", "# La fonction scipy.special.factorial donne la factorielle.\n", "# La fonction scipy.special.hermite(n) donne le n-ème polynôme de Hermite H_n.\n", "def h(x, n):\n", " Nn = (2**n * scipy.special.factorial(n))**(-1/2) * (omega / np.pi)**(1/4)\n", " return Nn * np.exp(- omega * x**2 / 2) * scipy.special.hermite(n)(np.sqrt(omega) * x)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Avant de les utiliser, vérifions numériquement que les fonctions $h_n(x)$ ainsi définies sont bien normalisées: il faut que $$\\int_{-\\infty}^\\infty |h_n(x)|^2\\,dx = 1\\,.$$ Avec la fonction **scipy.integrate.quad(f, a, b)** on peut calculer des intégrales numériques, même impropres. Elle renverra deux nombres: une estimation numérique de l'intégrale de la fonction $f$ entre $a$ et $b$, et une estimation de l'erreur numérique." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.integrate\n", "for n in range(3):\n", " def integrande(x):\n", " return h(x, n)**2\n", " print(\"Normalisation de h\" + str(n) + \":\", scipy.integrate.quad(integrande, -np.inf, np.inf))\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La normalisation est alors correcte dans les limites de la précision numérique dont nous avons besoin. \n", "\n", "Maintenant on peut comparer les deux fonctions d'onde pour les deux états fondamentaux, celui de l'oscillateur harmonique et celui du potentiel de Lennard-Jones. (On rappelle que le signe de la fonction d'onde n'est pas physique, alors il faut éventuellement multiplier par $-1$ pour mieux pouvoir comparer.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(x, fonctionsdonde[0][:251], label=r\"$\\psi_0(x)$\")\n", "h0 = h(x - x0, 0)\n", "plt.plot(x, h0 , label=r\"$h_0(x)$\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour l'état fondamental le profil est au moins qualitativement similaire, même si la fonction d'onde de l'oscillateur harmonique est légèrement décalé vers la gauche et moins large. Par contre les états excités ne se ressemblent pas beaucoup:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(x, fonctionsdonde[1][:251], label=r\"$\\psi_1(x)$\")\n", "h1 = h(x - x0, 1)\n", "plt.plot(x, h1, label=r\"$h_1(x)$\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(x, -fonctionsdonde[2][:251], label=r\"$\\psi_2(x)$\")\n", "h2 = h(x - x0, 2)\n", "plt.plot(x, h2, label=r\"$h_2(x)$\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### L'approximation WKB" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour les états liés d'énergie $E$ telle qu'il y a deux points tournants classiques $x_a$ et $x_b$, la méthode WKB donne lieu à une condition de quantification sur $E$: $$\\frac{1}{\\hbar}\\int_{x_a}^{x_b}{\\rm d}x'\\;\\sqrt{2m(E-V(x'))} = \\left(n+\\frac{1}{2}\\right)\\pi\\,,\\qquad\\qquad n\\in\\mathbb{N}$$. \n", "\n", "Notons que les points tournants $x$ pour le potentiel de Lennard-Jones vérifient $$E=\\epsilon\\,\\left(\\left(\\frac{x_0}{x}\\right)^{12}-2\\,\\left(\\frac{x_0}{x}\\right)^6\\right)\\,.$$ En posant $y=(x_0/x)^6$, on a alors $y^2-2y-\\epsilon/E=0$ et donc $$y=1\\pm\\sqrt{1+\\frac{E}{\\epsilon}}\\,,$$ soit $$x_{a,b}=x_0\\left(1\\pm\\sqrt{1+\\frac{E}{\\epsilon}}\\right)^{-1/6}\\,.$$\n", "Pour trouver les énergies qui vérifient la relation ci-dessus, on va définir une fonction $$f(E)=\\frac{1}{\\hbar}\\int_{x_a(E)}^{x_b(E)}{\\rm d}x'\\;\\sqrt{2m(E-V(x'))}-\\left(n+\\frac{1}{2}\\right)\\pi$$ et chercher un zéro de $f$ numériquement. On pourrait le faire, par exemple, avec la méthode de Newton, mais utilisons plutôt la bibliothèque **scipy** encore, plus précisément la méthode **scipy.optimize.root**. A cette dernière il faut fournir deux arguments: la fonction numérique dont on cherche le zéro, et la valeur de départ (qui peut être un seul nombre comme ici ou tout un tableau 1-dimensionnel, si la fonction est définie sur les vecteurs). Elle retourne un \"objet\", c.-à-d. une structure de données définie dans la même bibliothèque, dont le \"champ\" **x** contient la valeur du zéro." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def f(E, n):\n", " xa = x0 * (1 + np.sqrt(1+E/epsilon))**(-1/6)\n", " xb = x0 * (1 - np.sqrt(1+E/epsilon))**(-1/6)\n", " integrande = lambda x: np.sqrt(2 * (E - V(x, epsilon)))\n", " return scipy.integrate.quad(integrande, xa, xb)[0] - (n + 1/2)*np.pi\n", "\n", "import scipy.optimize\n", "\n", "EWKB = np.empty((3))\n", "\n", "for n in range(3):\n", " EWKB[n] = scipy.optimize.root(lambda e: f(e, n), -10).x\n", "print(EWKB)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pas trop mal: les valeurs avant étaient $-62.75$, $-17.99$ et $-1.82$, alors on a une précision rélative de l'ordre de 1%. Visiblement la méthode est tombée sur une valeur négative de $y$ à un moment donné pendant la recherche, d'où le message d'avertissement.\n", "\n", "On va enfin construire la fonction d'onde dans l'approximation WKB qui vérifie\n", "$$ \\psi(x)= \\frac{\\pm A}{(V(x)-E)^{1/4}}\\exp\\left(-\\frac{1}{\\hbar}\\int_x^{x_a}{\\rm d}x'\\,\\sqrt{2m(V(x')-E)}\\right)\\qquad (xx_b)$$\n", "Pour rappel, l'approximation WKB n'est pas valide si on est trop proche des points tournants, alors on ne s'attend pas que la fonction d'onde sera bien réproduite partout." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def psiWKB(n, x, A=1): # la normalisation A est laissée libre pour l'instant.\n", " E = EWKB[n]\n", " xa = x0 * (1 + np.sqrt(1+E/epsilon))**(-1/6)\n", " xb = x0 * (1 - np.sqrt(1+E/epsilon))**(-1/6)\n", " integrande = lambda x: np.sqrt(2 * np.abs(E - V(x, epsilon)))\n", " prefacteur = A / (np.abs(V(x, epsilon)))**(1/4)\n", " # les signes dans la suite sont choisis afin d'avoir la fonction d'onde psi[1, x] à peu près continue\n", " if x < xa:\n", " return -prefacteur * np.exp(-scipy.integrate.quad(integrande, x, xa)[0])\n", " elif x > xb:\n", " return prefacteur * np.exp(-scipy.integrate.quad(integrande, xb, x)[0])\n", " else:\n", " return -2 * prefacteur * np.cos(scipy.integrate.quad(integrande, xa, x)[0] - np.pi/4)\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = fonctionsdonde[1][50] / psiWKB(1, 1)\n", "X = np.linspace(.1, 3, 500)\n", "psi1 = np.array([psiWKB(1, x, A) for x in X])\n", "plt.plot(X, psi1, label=r\"$\\psi_1(x)$ (WKB)\")\n", "plt.plot(np.linspace(0, 3, 151), fonctionsdonde[1][:151],label=r\"$\\psi_1(x)$\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "C'est sûrement mieux que l'approximation harmonique, mais loin d'être parfait. On note en particulier le comportement singulier autour des points tournants à $x\\approx 0.9$ et $x\\approx 1.6$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 2 }