[Calcul scientifique / technique par Python] Résolution de l'équation différentielle ordinaire du second ordre par la méthode Numerov, calcul numérique

introduction

Dans les calculs scientifiques et techniques, les équations différentielles ordinaires du second ordre qui n'incluent pas la différenciation du premier ordre,

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

Apparaît souvent (comme l'équation de Schrödinger unidimensionnelle).

** Pour résoudre cette équation, il existe un algorithme explicite très simple et efficace appelé ** méthode Numerov **. Cette méthode n'est que du premier ordre plus précise que la méthode Runge-Kutta du quatrième ordre [1]. ** **

Dans cet article, nous utiliserons Python pour résoudre un exemple simple.

Méthode Numerov

Espacement uniforme de la grille

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

Comme, la solution par la méthode Numerov est

$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]. En particulier, lorsque $ S (x) = 0 $, [l'équation de Helmholtz] unidimensionnelle (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), et la solution par la méthode Numerov est

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

Sera.

exemple

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

Est résolu de x = 0 à 1.

La solution exacte est

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

** La méthode Numerov nécessite une valeur de $ y [1] $ en plus de $ y [0] $. ** ** Ici, à partir de la représentation différentielle de $ y '(0) $ Soit $ y '(0) = (y [1] -y [0]) / h = 0 $ et $ y [1] = 1 $.

Définissez également $ h = x_ {n} -x_ {n + 1} = 0,005 $

code

"""
Résolution d'équations différentielles ordinaires du second ordre par la méthode Numerov

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

#Conditions initiales
y[0] = 1
y[1]=1  # y'(0) = (y[1]-y[0])/delta_Évaluer par différence avant avec x

def Numerov (N,delta_x,k2,u):  #Développement par la méthode Numerov
    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) #Solution de la méthode Numerov

# 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') #étiquette de l'axe des x
plt.ylabel('Y') #étiquette de l'axe y

plt.show()

résultat

Les points représentent la solution numérique par la méthode Numerov et la ligne rouge représente la solution exacte.

Les références

[1] Article Qiita [Calcul scientifique / technique par Python] Résolution de l'équation de Newton unidimensionnelle par la méthode Runge-Kutta du 4ème ordre

[2] Masanori Abe et al. (Traduction), ["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)