[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation

Introduction

In scientific and technical calculations, second-order ordinary differential equations that do not include first-order differentials,

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

Often appears (such as the one-dimensional Schrodinger equation).

** As a solution to this equation, there is a very simple and efficient algorithm called ** Numerov method **. This method is only first-order more accurate than the fourth-order Runge-Kutta method [1]. ** **

In this article, we will use Python to solve a simple example.

Numerov method

Uniform grid spacing

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

As, the solution by the Numerov method is

$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]. In particular, when $ S (x) = 0 $, the one-dimensional [Helmholtz equation](https://ja.wikipedia.org/wiki/%E3%83%98%E3%83%AB%E3%83] % A0% E3% 83% 9B% E3% 83% AB% E3% 83% 84% E6% 96% B9% E7% A8% 8B% E5% BC% 8F), and the solution by the Numerov method is

$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} $

Will be.

example

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

Is solved from x = 0 to 1.

The exact solution is

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

** The Numerov method requires a value of $ y [1] $ in addition to $ y [0] $. ** ** Here, from the forward difference representation of $ y'(0) $, Let $ y'(0) = (y [1] -y [0]) / h = 0 $ and $ y [1] = 1 $.

Also, set $ h = x_ {n} -x_ {n + 1} = 0.005 $

code

"""
Solving second-order ordinary differential equations by the Numerov method

"""
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])

#Initial conditions
y[0] = 1
y[1]=1  # y'(0) = (y[1]-y[0])/delta_Evaluate by forward difference with x

def Numerov (N,delta_x,k2,u):  #Development by the Numerov method
    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) #Numerov method solution

# 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') #x-axis label
plt.ylabel('Y') #y-axis label

plt.show()

result

The points represent the numerical solution by the Numerov method, and the red line represents the exact solution.

References

[1] Qiita article [Science and technology calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method

[2] Masanori Abe et al. (Translation), ["Kunin Computer Physics"](https://www.amazon.co.jp/%E8%A8%88%E7%AE%97%E6%A9%9F% E7% 89% A9% E7% 90% 86% E5% AD% A6-Steven-Koonin / dp / 4320032918)