NumPy・SciPyを用いた数値計算の高速化 : 応用その1 ref: http://qiita.com/jabberwocky0139/items/a9751d11caa64bc19226

file0.txt

\frac{d^2}{dx^2}\psi(x) = \frac{\psi(x + \Delta x) - 2\psi(x) + \psi(x - \Delta x)}{\Delta x^2} + {\cal O}(\Delta x^2)

file1.txt

\psi(x)\rightarrow \{\psi(x_0), \psi(x_1), ..., \psi(x_{n-1}), \psi(x_n)\} = \{\psi_0, \psi_1, ..., \psi_{n-1}, \psi_n\}\hspace{1cm}x_i = -L/2 + i\Delta x

file10.txt

>>> from scipy.integrate import quad
>>> quad(f, -np.inf, np.inf)
(1.7724538509055159, 1.4202636780944923e-08)

file11.txt

H = \frac{\hat{p}^2}{2} + \frac{\hat{q}^2}{2} = -\frac12\frac{d^2}{dq^2} + \frac{q^2}{2}

file12.txt

H\psi_\ell(q) = E_\ell\psi_\ell(q)

file13.txt

    \frac{q^2}{2}\psi(q) \rightarrow \frac12
    \begin{pmatrix}
    q_0^2&0&0&\cdots&0\\ 
    0&q_1^2 &0&\cdots&0\\
    0 & 0&q_2^2&& \vdots\\
    \vdots&&&\ddots&0\\
    0& \cdots& 0 &0& q_n^2
    \end{pmatrix}
    \begin{pmatrix}
    \psi_0\\
    \psi_1\\
    \vdots\\
    \psi_{n-1}\\
    \psi_n
    \end{pmatrix}=\frac12
    \begin{pmatrix}
    (-L/2)^2&0&0&\cdots&0\\ 
    0&(-L/2 + \Delta q)^2 &0&&0\\
    0 & 0&(-L/2 + 2\Delta q)^2&& \vdots\\
    \vdots&&&\ddots&0\\
    0& \cdots& 0 &0& (L/2)^2
    \end{pmatrix}
    \begin{pmatrix}
    \psi_0\\
    \psi_1\\
    \vdots\\
    \psi_{n-1}\\
    \psi_n
    \end{pmatrix} = V\psi

file14.txt

H\psi = (K + V)\psi = \frac12
\begin{pmatrix}
\frac{2}{\Delta q^2} + (-L/2)^2&\frac{-1}{\Delta q^2}&0&\cdots&0\\ 
\frac{-1}{\Delta q^2}&\frac{2}{\Delta q^2} + (-L/2 + \Delta q)^2 &\frac{-1}{\Delta q^2}&&0\\
0 & \frac{-1}{\Delta q^2}&\frac{2}{\Delta q^2} + (-L/2 + 2\Delta q)^2&& \vdots\\
\vdots&&&\ddots&\frac{-1}{\Delta q^2}\\
0& \cdots& 0 &\frac{-1}{\Delta q^2}& \frac{2}{\Delta q^2} + (L/2)^2
\end{pmatrix}
\begin{pmatrix}
\psi_0\\
\psi_1\\
\vdots\\
\psi_{n-1}\\
\psi_n
\end{pmatrix} = E_\ell\psi

file15.txt

# 系の設定
>>> L, N = 10, 200
>>> x, dx = np.linspace(-L/2, L/2, N), L / N

# 運動項K
>>> K = np.eye(N, N)
>>> K_sub = np.vstack((K[1:], np.array([0] * N)))
>>> K = dx**-2 * (2 * K - K_sub - K_sub.T)

# ポテンシャル項
>>> V = np.diag(np.linspace(-L/2, L/2, N)**2)

# エルミート行列の固有値方程式
# wが固有値, vが固有ベクトル
>>> H = (K + V) / 2
>>> w, v = np.linalg.eigh(H)

file16.txt

>>> np.diff(w)[:20]
[ 1.00470938  1.00439356  1.00407832  1.00376982  1.00351013  1.00351647
  1.00464172  1.00938694  1.02300462  1.05279364  1.10422331  1.17688369
  1.26494763  1.36132769  1.46086643  1.56082387  1.66005567  1.75820017
  1.85521207  1.95115074]

file17.txt

\frac{\partial}{\partial t}f(x, t) = \frac{\partial^2}{\partial x^2}f(x, t)

file18.txt

f(x, 0) = \exp(-x^2)\\
f_k \equiv f(k\Delta x, t), \hspace{0.7cm} \Delta x = L/N

file19.txt

\frac{d}{dt}f_0 = \frac{f_1 - 2f_0}{\Delta x^2}\\
\frac{d}{dt}f_1 = \frac{f_2 - 2f_1 + f_0}{\Delta x^2}
...\\
\frac{d}{dt}f_k = \frac{f_{k+1} - 2f_k + f_{k-1}}{\Delta x^2}\\
...\\
\frac{d}{dt}f_{n-1} = \frac{f_{n} - 2f_{n-1} + f_{n-2}}{\Delta x^2}\\
\frac{d}{dt}f_n = \frac{-2f_n + f_{n-1}}{\Delta x^2}

file2.txt

\frac{d^2}{dx^2}\psi(x) = \frac{1}{\Delta x^2}
\begin{pmatrix}
\psi_{-1} -2\psi_0 + \psi_1\\
\psi_0 -2\psi_1 + \psi_2\\
\vdots\\
\psi_{n-2} -2\psi_{n-1} + \psi_n\\
\psi_{n-1} -2\psi_n + \psi_{n+1}
\end{pmatrix}
\simeq \frac{1}{\Delta x^2}
\begin{pmatrix}
-2&1&0&\cdots&0\\ 
1&-2 &1&&0\\
0 & 1&-2&& \vdots\\
\vdots&&&\ddots&1\\
0& \cdots& 0 &1& -2
\end{pmatrix}
\begin{pmatrix}
\psi_0\\
\psi_1\\
\vdots\\
\psi_{n-1}\\
\psi_n
\end{pmatrix}
=K\psi

file20.txt

from scipy.integrate import odeint

def equation(f, t=0, N, L):
    dx = L / N
    gamma = dx**-2
    i = np.arange(1, N, 1)

    # f0
    arr = np.array(gamma * (f[1] - 2 * f[0]))
    # f1 to fN-1
    arr = np.append(arr, gamma * (f[i+1] - 2 * f[i] + f[i-1]))
    # fN
    arr = np.append(arr, gamma * (-2 * f[N] + f[N-1]))

    return arr

file21.txt

from scipy.integrate import odeint

def equation(f, t=0, N, L):
    dx = L / N
    arr = np.gradient(np.gradient(f, dx), dx)

    return arr

file22.txt

def f(x):
    return np.exp(-x**2)

file23.txt

# initial parameter(optional)
N, L = 200, 40

# coordinate
q = np.linspace(-L/2, L/2, N)

# initial value for each fk
fk_0 = f(q)

# time
t_max, t_div = 10, 5
t = np.linspace(0, t_max, t_div)

trajectories = odeint(equation, fk_0, t, args=(N, L))

file3.txt

>>> import numpy as np
>>> def f(x):
...     return np.exp(-x**2)
... 
>>> L, N = 7, 100
>>> x = np.linspace(-L/2, L/2, N)
>>> psi = f(x)
>>> psi
array([  4.78511739e-06,   7.81043424e-06,   1.26216247e-05,
         2.01935578e-05,   3.19865883e-05,   5.01626530e-05,
         7.78844169e-05,   1.19723153e-04,   1.82206228e-04,
         ...,
         2.74540100e-04,   1.82206228e-04,   1.19723153e-04,
         7.78844169e-05,   5.01626530e-05,   3.19865883e-05,
         2.01935578e-05,   1.26216247e-05,   7.81043424e-06,
         4.78511739e-06])

file4.txt

>>> K = np.eye(N, N)
>>> K
array([[ 1.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  1.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  1., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  1.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  1.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  1.]])
>>> K_sub = np.vstack((K[1:], np.array([0] * N)))
>>> K_sub
array([[ 0.,  1.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  1., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  1.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  1.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
>>> K = -2 * K + K_sub + K_sub.T
>>> K
array([[-2.,  1.,  0., ...,  0.,  0.,  0.],
       [ 1., -2.,  1., ...,  0.,  0.,  0.],
       [ 0.,  1., -2., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ..., -2.,  1.,  0.],
       [ 0.,  0.,  0., ...,  1., -2.,  1.],
       [ 0.,  0.,  0., ...,  0.,  1., -2.]])

file5.txt

dx = L/N
psi_2dot = dx**-2 * np.dot(K, psi)

file6.txt

psi_dot_np = np.gradient(psi, dx)
psi_2dot_np = np.gradient(psi_dot_np, dx)

file7.txt

\int_a^b dx\; f(x) = \sum_{i=0}^{n-1}\frac{h}{6}\left[f(x_i) + 4f\left(x_i + \frac{h}{2}\right) + f(x_i + h)\right]\hspace{1cm}x_i = a + ih

file8.txt

>>> simp_arr = f(x) + 4 * f(x + dx / 2) + f(x + dx)
>>> dx / 6 * np.sum(simp_arr)
1.75472827313

file9.txt

>>> from scipy.integrate import simps
>>> simps(psi, x)
1.7724525416390826