[Wissenschaftlich-technische Berechnung durch Python] Lösung des Randwertproblems gewöhnlicher Differentialgleichungen im Matrixformat, numerische Berechnung

Einführung

** Die normale Differentialgleichung für $ y (x) $ ist homogen linear und die Parameter (proprietärer Wert lauten wie folgt. % 9C% 89% E5% 80% A4)) Wenn $ \ lambda $ ebenfalls linear ist, dann das Randwertproblem (https://ja.wikipedia.org/wiki/%E5%A2%83%E7) % 95% 8C% E5% 80% A4% E5% 95% 8F% E9% A1% 8C) wird auf die Eigenwertfrage der Matrix reduziert und durch Lösen der Matrixgleichung die Lösung (Eigenfunktion) und den Eigenwert der gewöhnlichen Differentialgleichung Es ist möglich, [1] zu bestimmen. ** ** **

** Es ist einfacher zu implementieren als die Methode zur Integration und Lösung gewöhnlicher Differentialgleichungen als Anfangswertproblem (Aufnahmemethode [2] usw.). ** ** **

Betrachten wir hier die gewöhnliche Differentialgleichung zweiter Ordnung.

p(x)y''(x)+q(x)y'(x)+r(x)y(x) = \lambda (u(x)y''(x)+v(x)y'(x)+w(x)y(x)) {\tag 1}

Als gleichzeitige Randbedingung

y(x_0)=0, y(x_N)=0 {\tag2}

Denk an.

Teilen Sie das Intervall $ x = [x_0, x_1] $ in N gleiche Teile

\delta x = (x_1-x_0)/N{\tag 3}

x_i=x_0+i\ \delta x{\tag 4}

y_i=y(x_i) {\tag 5}

Und.

Wenn nun die Ableitungen $ y '' (x) $ und $ y '(x) $ durch die Mitteldifferenznäherung ** ausgewertet werden

y'(x) = \frac{y_{i+1}-y_{i-1}}{2 \delta x} +O(\delta x^3){\tag 6}

y''(x) = \frac{y_{i+1}-2y_{i}+y_{i-1}}{\delta x ^2}+O(\delta x^3) {\tag 7}

Kann ausgedrückt werden als. Der Eigenwert / die Eigenfunktion, wenn (1) durch diese Näherungen gelöst wird, enthält einen Fehler von $ O (\ delta x ^ 2) $.

Durch Einsetzen von Gleichung (5,6) in Gleichung (1) und Umordnen dieser Gleichungen kann das Problem als ** allgemeines Eigenwertproblem ** im Matrixformat ausgedrückt werden.

$ A \ \mathbf{y} = \lambda \ B \ \mathbf{y} \ {\tag 8}$

Wobei $ \ mathbf {y} $ der Lösungsvektor ist. \mathbf{y}=(y_1, y_2, y_3, ...., y_{N-1}) {\tag 9}

$ A = (a_ {ij}) $ und $ B = (b_ {ij}) $ sind $ N-1 $ Dimensionsmatrizen und dreifach diagonale Matrizen. Die Matrixelemente sind $ i = 1,2, ..., N-1 $

a_{i,i-1}=\frac{p_i}{\delta x^2}+\frac{q_i}{2 \delta x}, \ a_{i,i}=\frac{-2p_i}{\delta x^2}+r_i, \ a_{i,i+1}= \frac{p_i}{\delta x^2}-\frac{q_i}{2 \delta x} {\tag {10}} a_{i,j} = 0 \ ( j\neq i-1, i, i+1){\tag {11}}

b_{i,i-1}=\frac{u_i}{\delta x^2}+\frac{v_i}{2 \delta x}, \ b_{i,i}=\frac{-2u_i}{\delta x^2}+w_i, \ b_{i,i+1}= \frac{u_i}{\delta x^2}-\frac{v_i}{2 \delta x} {\tag {12}} b_{i,j} = 0 \ ( j\neq i-1, i, i+1){\tag {13}}

Es wird ausgedrückt als.

** Gleichung (8) kann mit Python-Bibliotheken (numpy, scipy usw.) gelöst werden [3]. ** ** **

In diesem Artikel werde ich versuchen, ein einfaches Beispiel mit dieser Methode zu lösen.

Diese Methode unterscheidet sich von der Funktionserweiterungsmethode (** Spektrumsmethode **), bei der die Lösung durch ein geeignetes Funktionssystem erweitert wird und die Matrixgleichung für den Koeffizienten erhalten wird.

Inhalt

** Randwertproblem der homogenen normalen Differentialgleichung zweiter Ordnung (einfache Schwingung) **

y''(x)=-\lambda y, \ (y(0)=y(1)=0)

Löse ** mit der Matrixmethode. ** ** **

In diesem Fall ist jede Funktion in Gleichung (1)

p(x)=1, q(x)=0, r(x)=0, u(x)=0, v(x)=0, w(x)=1 {\tag {14}}

Und

In $ A \ \ mathbf {y} = \ lambda \ B \ \ mathbf {y} $ ist $ A $

a_{i,i-1}=\frac{1}{\delta x^2}, \ a_{i,i}=\frac{-2}{\delta x^2}, \ a_{i,i+1}= \frac{1}{\delta x^2} {\tag {15}} Und

$ B $ ist ** Unit Matrix * * Um $ E $

B = -E{\tag {16}} Es wird ausgedrückt als.

** Dann für diese Matrizen ** $ A \ \mathbf{y_n} = - \lambda_n \ E \ \mathbf{y} {}\tag {17} $

Löse **, um den Eigenwert $ \ lambda_n $ und die entsprechende Eigenfunktion $ y_n (x) $ zu finden. ** ** **

Code

Lösen Sie das allgemeine Eigenwertproblem (Gleichung 17) mit ** scipy.linalg.eig ** [3].

"""
Lösen des Randwertproblems mit der Matrixmethode
"""
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt

delta_x = 0.025
x0, x1 = 0, 1 #Grenze x Koordinate
N=int((x1-x0)/delta_x) #Anzahl der Datengitter
print("N=",N)

y = np.zeros([N-1,N+1])

#Gleichzeitige Randbedingung
y[:,0] = 0
y[:,-1] = 0

A=np.zeros([N-1,N-1])
B=-np.identity(N-1) # B= -E

#Konstruktion der Matrix A.
for i in range(N-1):  #Achten Sie auf die Indexposition, da es sich um eine dreifach diagonale Matrix handelt.
    if i == 0:
        A[i,i] = -2/(delta_x**2) 
        A[i,i+1] = 1/(delta_x**2)
    elif i == N-2:
        A[i,i-1] = 1/(delta_x**2)
        A[i,i] = -2/(delta_x**2) 
    else:
        A[i,i-1] = 1/(delta_x**2)
        A[i,i] = -2/(delta_x**2)
        A[i,i+1] = 1/(delta_x**2)
#
        
eigen, vec=  scipy.linalg.eig(A,B) #Lösen Sie das allgemeine Eigenwertproblem und finden Sie den Eigenwert"eigen"， Spezifische Funktion"vec"Im Objekt speichern

print("eigen values=",eigen)
#print("eigen vectors=",vec)

for j in range(N-1):
    for i in range(1,N):
        y[j, i] = vec[i-1,j]

#
# for plot
X= np.linspace(x0,x1, N+1)
plt.plot(X, y[0,:],'-',markersize=5,label='y1')
plt.plot(X, y[1,:],'-',markersize=5,label='y2')
plt.plot(X, y[2,:],'-',markersize=5,label='y3')
plt.legend(loc='upper right')
plt.xlabel('X') #x-Achsenbeschriftung
plt.ylabel('Y') #y-Achsenbeschriftung

plt.show()

Ergebnis

Die ersten drei werden in aufsteigender Reihenfolge der eindeutigen Werte angezeigt. Von links $ \ lambda_1 $, $ \ lambda_2 $, $ \ lambda_3 $. eigen values= [ 9.86453205+0.j 39.39731010+0.j 88.41625473+0.j]

Die diesen Eigenwerten entsprechenden Eigenfunktionen $ y_1 (x) $, $ y_2 (x) $, $ y_3 (x) $ sind in der Abbildung dargestellt.

Verweise

[1] Ryo Natori, ["Numerische Analyse und ihre angewandte Computermathematik, Reihe 15"](https://www.amazon.co.jp/%E6%95%B0%E5%80%A4%E8%A7%A3 % E6% 9E% 90% E3% 81% A8% E3% 81% 9D% E3% 81% AE% E5% BF% 9C% E7% 94% A8% E3% 82% B3% E3% 83% B3% E3% 83% 94% E3% 83% A5% E3% 83% BC% E3% 82% BF% E6% 95% B0% E5% AD% A6% E3% 82% B7% E3% 83% AA% E3% 83% BC% E3% 82% BA-% E5% 90% 8D% E5% 8F% 96-% E4% BA% AE / dp / 4339025488), Corona, 1990.

[2] Beispiel für die Lösung gewöhnlicher Differentialgleichungen mit der Aufnahmemethode: [Wissenschaftliche / technische Berechnung mit Python] Lösen der eindimensionalen Schrödinger-Gleichung im stationären Zustand durch Aufnahmemethode (2), harmonisches Oszillatorpotential, Quantenmechanik //qiita.com/sci_Haru/items/edb475a6d2eb7e901905)

[3] So finden Sie die Eigenwerte / Eigenvektoren (Funktionen) einer Matrix mithilfe einer Bibliothek: [[Wissenschaftliche / technische Berechnung mit Python] Lösen (verallgemeinerter) Eigenwertprobleme mithilfe von numpy / scipy mithilfe der Bibliothek](http: / /qiita.com/sci_Haru/items/034c6f74d415c1c10d0b)