[Calcul scientifique / technique par Python] Résolution du problème de la valeur aux limites des équations différentielles ordinaires au format matriciel, calcul numérique

introduction

** L'équation différentielle normale pour $ y (x) $ est linéaire homogène, et les paramètres ([valeur propriétaire](https://ja.wikipedia.org/wiki/%E5%9B%BA%E6] sont les suivants. % 9C% 89% E5% 80% A4)) Si $ \ lambda $ est également linéaire, alors le problème de la valeur limite (https://ja.wikipedia.org/wiki/%E5%A2%83%E7) % 95% 8C% E5% 80% A4% E5% 95% 8F% E9% A1% 8C) est réduit à la question des valeurs propres de la matrice, et en résolvant l'équation matricielle, la solution (fonction propre) et la valeur propre de l'équation différentielle ordinaire Il est possible de déterminer [1]. ** **

** Il est plus facile à mettre en œuvre que la méthode d'intégration et de résolution d'équations différentielles ordinaires comme problème de valeur initiale (méthode de prise de vue [2], etc.). ** **

Ici, considérons l'équation différentielle ordinaire du second ordre.

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}

En tant que condition aux limites simultanée

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

penser à.

Divisez l'intervalle $ x = [x_0, x_1] $ en N parties égales

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

Et.

Maintenant, lorsque les dérivées $ y '' (x) $ et $ y '(x) $ sont évaluées par l'approximation de la différence de centre **

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}

Peut être exprimé comme. La valeur propre / fonction propre lorsque (1) est résolu par ces approximations inclut une erreur de $ O (\ delta x ^ 2) $.

En substituant les équations (5,6) dans l'équation (1) et en les réorganisant, le problème peut être exprimé sous la forme d'un ** problème général de valeurs propres ** au format matriciel.

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

Où $ \ mathbf {y} $ est le vecteur de solution. \mathbf{y}=(y_1, y_2, y_3, ...., y_{N-1}) {\tag 9}

$ A = (a_ {ij}) $ et $ B = (b_ {ij}) $ sont des matrices de dimension $ N-1 $ et des matrices triples diagonales. Ces éléments de matrice sont $ 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}}

Il est exprimé comme.

** L'équation (8) peut être résolue en utilisant les bibliothèques Python (numpy, scipy, etc.) [3]. ** **

Dans cet article, je vais essayer de résoudre un exemple simple en utilisant cette méthode.

Cette méthode est différente de la méthode d'expansion des fonctions (** méthode du spectre **) dans laquelle la solution est développée par un système fonctionnel approprié et l'équation matricielle du coefficient est obtenue.

Contenu

** Problème de valeur aux limites de l'équation différentielle normale homogène du second ordre (vibration simple) **

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

Résolvez ** par la méthode matricielle. ** **

Dans ce cas, chaque fonction de l'équation (1) est

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

Et

Dans $ A \ \ mathbf {y} = \ lambda \ B \ \ mathbf {y} $, $ A $ est

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

$ B $ est ** Unit Matrix * * Par $ E $

B = -E{\tag {16}} Il est exprimé comme.

** Alors pour ces matrices ** $ A \ \mathbf{y_n} = - \lambda_n \ E \ \mathbf{y} {}\tag {17} $

Résolvez ** pour trouver la valeur propre $ \ lambda_n $ et la fonction propre correspondante $ y_n (x) $. ** **

code

Résolvez le problème général des valeurs propres (équation 17) en utilisant ** scipy.linalg.eig ** [3].

"""
Résolution du problème de la valeur limite par la méthode matricielle
"""
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt

delta_x = 0.025
x0, x1 = 0, 1 #Coordonnée frontière x
N=int((x1-x0)/delta_x) #Nombre de grilles de données
print("N=",N)

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

#Condition aux limites simultanée
y[:,0] = 0
y[:,-1] = 0

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

#Construction de la matrice A
for i in range(N-1):  #Faites attention à la position de l'index car il s'agit d'une matrice triple diagonale.
    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) #Résolvez le problème général des valeurs propres et trouvez la valeur propre"eigen"， Fonction spécifique"vec"Stocker dans l'objet

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

plt.show()

résultat

Les trois premiers sont affichés dans l'ordre croissant des valeurs uniques. À partir de la gauche, $ \ lambda_1 $, $ \ lambda_2 $, $ \ lambda_3 $. eigen values= [ 9.86453205+0.j 39.39731010+0.j 88.41625473+0.j]

Les fonctions propres $ y_1 (x) $, $ y_2 (x) $, $ y_3 (x) $ correspondant à ces valeurs propres sont représentées sur la figure.

Les références

[1] Ryo Natori, ["Numerical Analysis and Its Applied Computer Mathematics Series 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] Exemple de résolution d'équations différentielles ordinaires à l'aide de la méthode de prise de vue: [Calcul scientifique / technique par Python] Résolution d'une équation de Schrodinger unidimensionnelle en régime permanent par méthode de tir (2), potentiel d'oscillateur harmonique, mécanique quantique //qiita.com/sci_Haru/items/edb475a6d2eb7e901905)

[3] Comment trouver les valeurs propres / vecteurs propres (fonctions) d'une matrice à l'aide d'une bibliothèque: [[Calcul scientifique / technique par Python] Résolution de problèmes de valeurs propres (généralisés) en utilisant numpy / scipy, en utilisant la bibliothèque](http: / /qiita.com/sci_Haru/items/034c6f74d415c1c10d0b)