QM 07.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二次元シュレーディンガー方程式の解：直交座標系で解く\n",
    "今回は、二次元シュレーディンガー方程式をもう少し力技で解いてみよう。\n",
    "二次元シュレーディンガー方程式は\n",
    "$$\n",
    "\\left( -\\frac{\\hbar^2}{2m} {\\bf \\nabla}^2 + V({\\bf r}) \\right)\\psi({\\bf r}) = \\epsilon \\psi({\\bf r})\n",
    "$$\n",
    "となる。ここで、ナブラ演算子${\\bf \\nabla}^2$は、直交座標系では\n",
    "$$\n",
    "{\\bf \\nabla}^2 = \\frac{\\partial^2}{\\partial x^2} + \\frac{\\partial^2}{\\partial y^2}\n",
    "$$\n",
    "である。\n",
    "\n",
    "これまでは系に回転対称性があるとして円筒座標系を用いて解いた。ここでは、回転対称性を課さずに、直交座標系で直接解いてみよう。\n",
    "これまでのノートと同じようにするために、半径$R=\\sqrt{x^2+y^2}$の円を考え、それより外側には非常に大きなポテンシャルがあるとする。それ以外にはポテンシャルはないとする。つまり、これまでのノートでのポテンシャルがない場合に相当する。\n",
    "計算する領域は幅$L$の正方形とする。\n",
    "このとき、\n",
    "二階微分をQM01と同様に\n",
    "$$\n",
    "\\frac{d^2}{dx^2}\\psi |_{x=x_i} \\sim \\frac{\\psi(x_i + a) -2\\psi(x_i) +\\psi(x_i - a)}{a^2}\n",
    "$$\n",
    "と差分化する。yも同様である。半径$R$よりも遠い領域ではポテンシャルを強くいれるため、正方形の境界条件はどのようにしてもよい。ここでは、周期境界条件としておく。\n",
    "\n",
    "行列は"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "make_matH (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "using SparseArrays \n",
    "function make_matH(Nx,Ny,V,L)\n",
    "        N = Nx*Ny\n",
    "        mat_H = spzeros(ComplexF64,N,N)\n",
    "        ax = L/(Nx-1)\n",
    "        ay = L/(Ny-1)\n",
    "        for ix=1:Nx\n",
    "            x = ax*(ix-1)-L/2\n",
    "            for iy=1:Ny\n",
    "                i = (iy-1)*Nx+ix\n",
    "                y = ay*(iy-1)-L/2\n",
    "                for dx = -1:1\n",
    "                    jx = ix + dx\n",
    "                    jx +=ifelse(jx > Nx,-Nx,0)\n",
    "                    jx +=ifelse(jx < 1,Nx,0)\n",
    "                    tx = 0.0\n",
    "                    ty = 0.0\n",
    "                \n",
    "                    if dx==0\n",
    "                        tx = -2.0\n",
    "                        for dy = -1:1\n",
    "                            jy = iy + dy\n",
    "                            jy +=ifelse(jy > Ny,-Ny,0)\n",
    "                            jy +=ifelse(jy < 1,Ny,0)\n",
    "                            j = (jy-1)*Nx + jx\n",
    "                            if dy == 0\n",
    "                                ty = -2.0\n",
    "                                mat_H[i,j] = -(tx/ax^2+ty/ay^2)+V(x,y)\n",
    "                            else\n",
    "                                ty = 1.0\n",
    "                                mat_H[i,j] = -(ty/ay^2)\n",
    "                            end\n",
    "                        end\n",
    "                    else\n",
    "                        jy = iy\n",
    "                        j = (jy-1)*Nx + jx\n",
    "                        tx = 1.0\n",
    "                        mat_H[i,j] = -(tx/ax^2)\n",
    "                    end\n",
    "                                        \n",
    "                    \n",
    "                end\n",
    "            end\n",
    "        end\n",
    "        return mat_H\n",
    "    end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで、二次元座標系を考えるために行列のサイズは$N_x \\times N_y$と大きなものになるため、疎行列を使うこととした。\n",
    "疎行列とは、行列のうちゼロではない要素のみを格納した行列のことであり、メモリと計算量を節約することができる。\n",
    "疎行列の対角化にはArpackを用いるので、"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Arpack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "としておく。\n",
    "\n",
    "ポテンシャルは、半径R0以上で非常に大きな値をとる関数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "V (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const R0 = 10.0\n",
    "const μ = 0.0\n",
    "function V(x,y)\n",
    "    r2=x^2+y^2\n",
    "    v = -μ\n",
    "    if r2 > R0^2\n",
    "        v += 10000.0        \n",
    "    end\n",
    "    return v\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "としておく。\n",
    "\n",
    "行列を作成し、対角化すると、"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  5.131422 seconds (6.84 M allocations: 424.800 MiB, 6.60% gc time)\n",
      "[0.0569595, 0.144572, 0.144572, 0.259354, 0.259906, 0.299925, 0.400561, 0.400561, 0.484221, 0.484221, 0.566293, 0.566474, 0.695771, 0.697519, 0.736271, 0.75639, 0.75639, 0.936223, 0.936223, 0.97003]\n"
     ]
    }
   ],
   "source": [
    "Nx = 100\n",
    "Ny = 100\n",
    "L = 20.1\n",
    "mat_H = make_matH(Nx,Ny,V,L)\n",
    "@time e,v = eigs(mat_H,nev=20,which=:SM)\n",
    "println(sort(real(e)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "と固有値が得られる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QM04で解いた$n=0$の時の固有値は"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.057774, 0.304407, 0.748117, 1.389, 2.22705, 3.26228, 4.49467, 5.92421, 7.55088, 9.37468]\n"
     ]
    }
   ],
   "source": [
    "println([0.057774, 0.304407, 0.748117, 1.389, 2.22705, 3.26228, 4.49467, 5.92421, 7.55088, 9.37468])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "であり、$n=1$の固有値は"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.146673, 0.491691, 1.03395, 1.77341, 2.71005, 3.84386, 5.17484, 6.70296, 8.42821, 10.3506]\n"
     ]
    }
   ],
   "source": [
    "println([0.146673, 0.491691, 1.03395, 1.77341, 2.71005, 3.84386, 5.17484, 6.70296, 8.42821, 10.3506])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "である。二次元直交座標で解くことによって、異なる$n$の解がまとめて出てきていることがわかる。なお、$n=0$の最低固有値以外が縮退しているのは、$n$の正負で縮退があるからである。\n",
    "\n",
    "なお、もう少し精度をあげると、"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  3.451619 seconds (17.68 k allocations: 360.574 MiB, 4.70% gc time)\n",
      "[0.0573951, 0.145703, 0.145703, 0.261462, 0.261977, 0.302363, 0.403898, 0.403898, 0.488333, 0.488333, 0.571103, 0.571474, 0.702129, 0.703579, 0.742888, 0.763228, 0.763228, 0.945023, 0.945023, 0.978886]\n"
     ]
    }
   ],
   "source": [
    "Nx = 200\n",
    "Ny = 200\n",
    "L = 20.1\n",
    "mat_H = make_matH(Nx,Ny,V,L)\n",
    "@time e,v = eigs(mat_H,nev=20,which=:SM)\n",
    "println(sort(real(e)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さらに値は近づく。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.7.0",
   "language": "julia",
   "name": "julia-0.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}