소개

과학·기술 계산에서는, 1층 미분을 포함하지 않는 2층의 상미분 방정식,

$\frac{d^2 y}{d x^2} + k^2(x)y=S(x) {\tag 1} $

가 자주 나온다(1차원 슈레딩거 방정식 등).

이 방정식의 해법으로는 누메로프법이라고 불리는 매우 간단하고 효율적인 양해법의 알고리즘이 있다. 이 방법은 4차 룬게 쿠터법[1]보다 1차만 정밀도가 높다[2].

이 논문에서는 파이썬을 사용하여 간단한 예제를 풀 수 있습니다.

누메로프법

균일 한 격자 간격

$h = x_{n}-x_{n+1}{\tag 2}$

따라서 Numerov 방법에 의한 해결책은

$y_{n+1}=\frac{(1-5bh^2)y_n-(1+\frac{1}{2}k^2_{n-1} h^2)y_{n-1}+ b(S_{n+1}+10S_n+S_{n-1})}{1+bk^2_{n+1}h^2} +\mathcal{O(h^6)}{\tag 3} $

된다 [1]. 특히, $S(x)=0$일 때는 1차원의 헬름홀츠 방정식이 되고, 누메로프법에 의한 해는

$y_{n+1}=\frac{(1-5bh^2)y_n-(1+\frac{1}{2}k^2_{n-1} h^2)y_{n-1}} {1+bk^2_{n+1}h^2} +\mathcal{O(h^6)}{\tag 4} $

된다.

예제

$\frac{d^2 y}{d x^2} = -4\pi y, y(0)=1, y'(0)=0 {\tag 5} $

를 x=0에서 1까지 풉니다.

엄밀한 해는

$ y = cos(2\sqrt\pi x){\tag 6} $
이다.

누메로프법에서는 $y[0]$ 이외에 $y[1]$의 값이 필요하다.
여기서 $ y '(0) $의 전진 차분 표현에서,
$y'(0) = (y[1]-y[0])/h = 0 $로, $y[1] = 1$로 한다.

또, $h = x_{n}-x_{n+1} = 0.005$로 한다

코드

"""
ヌメロフ法による2階常微分方程式の解法

"""
import numpy as np
import matplotlib.pyplot as plt

delta_x=0.005
xL0, xR0  =0, 1
Nx =  int((xR0-xL0)/delta_x)

k2=np.zeros([Nx+1])
k2[:] = 4*np.pi

y=np.zeros([Nx])

#初期条件
y[0] = 1
y[1]=1  # y'(0) = (y[1]-y[0])/delta_xと前進差分して評価

def Numerov (N,delta_x,k2,u):  # ヌメロフ法による発展
    b = (delta_x**2)/12.0
    for i in range(1,N-1):
        u[i+1] = (2*u[i]*(1-5*b*k2[i])-(1+b*k2[i-1])*u[i-1])/(1+b*k2[i+1]) 



Numerov(Nx, delta_x, k2, y) # ヌメロフ法による解

# for plot
X= np.linspace(xL0,xR0, Nx)
y_exact = np.cos(2*np.sqrt(np.pi)*X)
plt.plot(X, y,'o',markersize=2,label='Numerov')
plt.plot(X, y_exact,'-',color='red',markersize=0.5,label='Exact')
plt.legend(loc='upper right')
plt.xlabel('X') # ｘ軸のラベル
plt.ylabel('Y') # y軸のラベル

plt.show()

결과

점이 누메로프법에 의한 수치해, 적선이 엄밀해를 나타낸다.

참고문헌

[1] Qiita 기사 [Python에 의한 과학·기술 계산] 4차의 룬게-쿠타법에 의한 1차원 뉴턴 방정식의 해법

[2] 아베 마사노리 외(역), 『쿠닌 계산기 물리학』