[Wissenschaftlich-technische Berechnung mit Python] Lösen der gewöhnlichen Differentialgleichung zweiter Ordnung nach der Numerov-Methode, numerische Berechnung

Einführung

In wissenschaftlichen und technischen Berechnungen werden gewöhnliche Differentialgleichungen zweiter Ordnung, die keine Differenzierung erster Ordnung enthalten,

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

Erscheint oft (wie die eindimensionale Schrödinger-Gleichung).

** Als Lösung für diese Gleichung gibt es einen sehr einfachen und effizienten expliziten Algorithmus namens ** Numerov-Methode **. Diese Methode ist nur erster Ordnung genauer als die Runge-Kutta-Methode vierter Ordnung [1]. ** ** **

In diesem Artikel verwenden wir Python, um ein einfaches Beispiel zu lösen.

Numerov-Methode

Gleichmäßiger Gitterabstand

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

As ist die Lösung nach der Numerov-Methode

$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]. Insbesondere wenn $ S (x) = 0 $ ist, ist die eindimensionale [Helmholtz-Gleichung](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), und die Lösung nach der Numerov-Methode ist

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

Wird sein.

Beispiel

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

Wird von x = 0 bis 1 gelöst.

Die genaue Lösung ist

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

** Die Numerov-Methode erfordert zusätzlich zu $ y [0] $ den Wert $ y [1] $. ** ** ** Hier aus der Vorwärtsdifferenzdarstellung von $ y '(0) $ Sei $ y '(0) = (y [1] -y [0]) / h = 0 $ und $ y [1] = 1 $.

Setzen Sie außerdem $ h = x_ {n} -x_ {n + 1} = 0,005 $

Code

"""
Lösen gewöhnlicher Differentialgleichungen zweiter Ordnung nach der Numerov-Methode

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

#Anfangsbedingungen
y[0] = 1
y[1]=1  # y'(0) = (y[1]-y[0])/delta_Durch Vorwärtsdifferenz mit x auswerten

def Numerov (N,delta_x,k2,u):  #Entwicklung nach der Numerov-Methode
    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-Methodenlösung

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

plt.show()

Ergebnis

Die Punkte repräsentieren die numerische Lösung nach der Numerov-Methode, und die rote Linie repräsentiert die exakte Lösung.

Verweise

[1] Qiita-Artikel [Wissenschaftliche / technische Berechnung von Python] Lösen der eindimensionalen Newton-Gleichung nach der Runge-Kutta-Methode 4. Ordnung

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