statio.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f50b8607",
   "metadata": {},
   "source": [
    "# Projet Numérique : Dynamique Quantique\n",
    "\n",
    "## Résolution en états stationnaire \n",
    "\n",
    "Numba : a voir vers la fin\n",
    "http://www.enit.rnu.tn/fr/Minds/quant/amphi6MQ.pdf pour explication sin et cos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84e2c304",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import special\n",
    "from scipy import constants\n",
    "import math\n",
    "\n",
    "#Nous faisons les calculs pour 1000 points\n",
    "N = 1000\n",
    "L = 20\n",
    "\n",
    "dx = L/N\n",
    "\n",
    "print(\"dx : \", dx)\n",
    "\n",
    "## On initialise x tel que \n",
    "\n",
    "x = np.linspace(-L/2,L/2,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e85605ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# On commence par former la liste de valeurs des vecteurs et V.\n",
    "\n",
    "V = np.zeros(N) #potentiel nul\n",
    "V2 = np.zeros(N) #potentiel harmonique\n",
    "\n",
    "omega =1\n",
    "masse = 1\n",
    "hbarre = 1\n",
    "\n",
    "for i in range(N):\n",
    "    V2[i] = masse*(omega*x[i])**2 / 2\n",
    "\n",
    "plt.plot(x,V)\n",
    "plt.plot(x,V2)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e4bf128",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Ensuite, on forme les matrices d et e diagonales principale et supérieure (resp) de H\n",
    "\n",
    "def get_matrice(V):\n",
    "\n",
    "    d = np.zeros(N)\n",
    "    e= np.zeros(N-1)\n",
    "\n",
    "    for i in range(V.size):\n",
    "        d[i] = 2/(dx**2) + V[i]\n",
    "\n",
    "    for i in range(e.size):\n",
    "        e[i] = -1/(dx**2)\n",
    "        \n",
    "    return d,e\n",
    "\n",
    "d,e = get_matrice(V)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c04c04d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import eigh_tridiagonal\n",
    "\n",
    "## Puis, on calcule w et v respectivement les valeurs propres et vecteurs propres de H\n",
    "\n",
    "w, v = eigh_tridiagonal(d,e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1b01961",
   "metadata": {},
   "outputs": [],
   "source": [
    "# On doit désormais normaliser les vecteurs propres :\n",
    "\n",
    "def normalize(m, dx):\n",
    "    m /= np.linalg.norm(m,axis=0)\n",
    "    m /= np.sqrt(dx)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1972a4a4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## On regarde la différence entre la courbe théorique et la courbe obtenu\n",
    "\n",
    "p = np.arange(0,N,1)\n",
    "\n",
    "Ep = np.zeros(N)\n",
    "for i in range(N):\n",
    "    Ep[i] = (np.pi*(i+1)/L)**2\n",
    "    \n",
    "#Pour voir plus distinctement la différence entre la courbe numérique et théorique, nous la mettons au carré\n",
    "R = (Ep-w)**2\n",
    "\n",
    "plt.plot(p,R,\"--\",label=\"$Erreur au carré$\")    \n",
    "plt.plot(p,w,label=\"$Numerical$\")\n",
    "plt.plot(p,Ep,label = \"$Theorical$\")\n",
    "plt.title('Courbe de la différence entre la fonction d\\'onde numérique et théorique')\n",
    "plt.xlabel('$p$')\n",
    "plt.ylabel('$Ep$')\n",
    "plt.xlim(0, 200)\n",
    "plt.ylim(0, 1000)\n",
    "plt.legend()\n",
    "plt.savefig('docs/statio/energie_num_vs_theoric_n%.i.pdf' %N,format='pdf')\n",
    "plt.savefig('docs/statio/energie_num_vs_theoric_n%.i' %N)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2af88075",
   "metadata": {},
   "outputs": [],
   "source": [
    "# On calcule le psi théorique pour le potentiel puit infini\n",
    "\n",
    "def psi_theo(p):\n",
    "    psi_theo = np.zeros(N)\n",
    "    for i in range(N):\n",
    "        if (p%2 == 0): #fonction pair\n",
    "            psi_theo[i] = np.sqrt(2/L)*np.cos(((p+1)*np.pi*x[i])/L+p/2*np.pi)\n",
    "        else: #fonction impair\n",
    "            psi_theo[i] = np.sqrt(2/L)*np.sin(((p+1)*np.pi*x[i])/L)\n",
    "    return psi_theo\n",
    "\n",
    "\n",
    "def energie_pot(p,w):\n",
    "    value_ep = np.zeros(N)\n",
    "    for i in range(value_ep.size):\n",
    "        value_ep[i] = w[p]\n",
    "    return value_ep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43a1db6b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# On crée la fonction pour former psi à partir de V\n",
    "\n",
    "def get_psi(V):\n",
    "    d,e = get_matrice(V)\n",
    "    \n",
    "    w, v = eigh_tridiagonal(d,e)\n",
    "    \n",
    "    for vec in v:\n",
    "        normalize(v,dx)\n",
    "    v0 = np.zeros((N,N))\n",
    "    for i in range(N):\n",
    "        for y in range(N):            #On échange ligne et colonnes\n",
    "            v0[i][y] = v[y][i]\n",
    "    \n",
    "    return w,v0\n",
    "\n",
    "# On crée la fonction pour afficher les graphiques\n",
    "\n",
    "def get_graph(w,v,name,theorical_psi,Ep):\n",
    "    for i in range(Ep):\n",
    "        plt.plot(x,v[i]+w[i],label=\"Numérique $p=%i$\" %i)\n",
    "        plt.plot(x,theorical_psi(i)+w[i],label = \"Théorique $p=%i$\" %i)\n",
    "        plt.plot(x,energie_pot(i,w),\"--\",label = \"Ep pour $p=%i$\" %i)\n",
    "    plt.xlabel('$x$')\n",
    "    plt.ylabel('$\\Psi(x)$')\n",
    "    plt.title('Fonction d\\'onde pour un potentiel '+name)\n",
    "    plt.legend(bbox_to_anchor=(1.05,1.0),loc='upper left')\n",
    "    plt.savefig('docs/statio/%s.pdf' %(name),bbox_inches='tight',format='pdf')\n",
    "    plt.savefig('docs/statio/%s' %(name),bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cf650e4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "w,v = get_psi(V)\n",
    "get_graph(w,v,\"puit_infini\",psi_theo,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ede9bfb8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#potentiel harmonique \n",
    "def factoriel(x):\n",
    "    if x <= 0:\n",
    "        return 1\n",
    "    else:\n",
    "        return x*factoriel(x-1)\n",
    "\n",
    "sigma = np.sqrt(np.sqrt(2.0))\n",
    "\n",
    "#On calcule la fonction d'onde théorique\n",
    "def psi_theo_harmo(p,y):\n",
    "    psi_theo = np.zeros(N)\n",
    "    p_hermite = special.hermite(p)\n",
    "    \n",
    "    for i in range(N):\n",
    "        psi_theo[i] = (1/(np.sqrt(2**p*factoriel(p)*np.sqrt(np.pi)*sigma)))*np.exp((-(x[i]+y)**2)/(2*sigma**2))*p_hermite((x[i]+y)/sigma)*(-1)**p\n",
    "    \n",
    "    return psi_theo\n",
    "\n",
    "w,v = get_psi(V2)\n",
    "\n",
    "#Fonction d'onde pour un potentiel harmonique centré en 0\n",
    "def psi_theo_harmonique(p):\n",
    "    return psi_theo_harmo(p,0)\n",
    "\n",
    "get_graph(w,v,\"harmonique\",psi_theo_harmonique,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41328000",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Double puit\n",
    "\n",
    "#Définition du potentiel\n",
    "V4 = np.zeros(N)\n",
    "a = 0.003\n",
    "for i in range(N):\n",
    "    V4[i] = a*(x[i]**4)-(x[i]**2/4)\n",
    "\n",
    "plt.plot(x,V4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29e8b405",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#psi_numerique\n",
    "w,v = get_psi(V4)\n",
    "\n",
    "#On calcule la fonction d'onde pour un double potentiel harmonique non couplé\n",
    "def psi_theo_doublepuits(p):\n",
    "    psi_theo = np.zeros(N)\n",
    "    b = 6.3\n",
    "    harmo_g = psi_theo_harmo(int(round(p/2.1,0)),b)\n",
    "    harmo_d = psi_theo_harmo(int(round(p/2.1,0)),-b)\n",
    "\n",
    "    for i in range(N):\n",
    "        psi_theo[i] = (harmo_g[i]+(-1)**p*harmo_d[i])\n",
    "        \n",
    "    \n",
    "    return normalize(psi_theo,dx)\n",
    "\n",
    "get_graph(w,v,\"double_puit\",psi_theo_doublepuits,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ccfb48a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Potentiel periodique\n",
    "\n",
    "V6 = np.zeros(N)\n",
    "Q = 50\n",
    "d = 1\n",
    "U = 3\n",
    "for i in range(N):\n",
    "    V6[i] = Q*np.cos(x[i]*np.pi/d)\n",
    "\n",
    "w,v = get_psi(V6)\n",
    "\n",
    "plt.plot(V6)\n",
    "plt.show()\n",
    "for i in range(U):\n",
    "    plt.plot(x,v[i]+w[i],label=\"Numérique $p=%i$\" %i)\n",
    "    plt.plot(x,energie_pot(i,w),\"--\",label = \"Ep pour $p=%i$\" %i)\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$\\Psi(x)$')\n",
    "plt.title('Fonction d\\'pour un potentiel harmonique periodique')\n",
    "plt.legend(bbox_to_anchor=(1.05,1.0),loc='upper left')\n",
    "plt.savefig('docs/statio/%s' %(\"periodique\"),bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b313eb1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}