Résolution d'une équation d'onde unidimensionnelle par la méthode de la différence (Python)

Equation d'onde simple

 \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \\
0<x<l,0<t

Je veux résoudre ça

condition limite

u(0)=0,\frac{\partial u}{\partial x}\Big|_{x=l} = 0

L'extrémité gauche est la condition aux limites de Diricre qui spécifie la valeur de $ u $ à la limite, et l'extrémité droite est la condition aux limites de Neumann qui spécifie la valeur différentielle de la valeur de limite.

La condition aux limites est que la gauche est l'extrémité fixe et la droite est l'extrémité libre.

Conditions initiales

u(x,0) = f(x)\\
\frac{\partial u}{\partial t}=0

S'il s'agit de la vibration d'une corde ou de la vibration d'une barre verticale, la vitesse initiale est de 0.

Discret

x_i = i\Delta x\\
t^n = n\Delta t\\
u_{i}^n = u(x_i,t^n)\\
(i = 1\cdots m, \Delta x = \frac{l}{m})

Et met. Cependant, $ m $ est le nombre de divisions spatiales.

Condition initiale Condition aux limites

Tout d'abord, définissez la solution initiale.

u_{i}^0 = f(x_i)\\
u_{i}^1 = u_{i}^0

Puisque le différentiel de temps du premier ordre est 0, $ u_ {i} ^ 1 = u_ {i} ^ 0 $.

Ensuite, définissez les conditions aux limites. Cependant, le problème est que la condition aux limites de Neumann à l'extrémité droite est exprimée par différenciation spatiale, donc deux éléments doivent être utilisés pour l'exprimer. Par conséquent, ajoutez un autre élément faux à la limite la plus à droite. En d'autres termes

u_{m+1}^n = u_{m}^n

Si tel est le cas, la condition aux limites de Neumann à l'extrémité droite, à savoir que la différenciation spatiale du premier ordre est 0, est satisfaite. Ensuite, la condition aux limites de Diricle la plus à gauche peut être exprimée comme suit.

u_1^n = 0

Enfin, l'équation elle-même sera discrétisée.

Discrimination de la différenciation spatiale

\frac{\partial^2 u}{\partial x^2}\Big|_{x=x_i}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}

Cela peut être dérivé en développant sur mesure $ u_ {i + 1} $ et $ u_ {i-1} $ avec $ x = x_i $.

Séparation du développement du temps

Penser de la même manière que la différenciation spatiale

\frac{\partial^2 u}{\partial t^2} =\frac{u^{n+1}-2u^{n}+u^{n-1}}{\Delta t^2}

Parce que c'est

\begin{align}
u^{n+1}_i &= 2u^{n}_i-u^{n-1}_i + \Delta t^2 \frac{\partial^2 u}{\partial t^2}\Big|_{x=x_i}\\
&= 2u^{n}_i-u^{n-1}_i + c^2 \Delta t^2 \frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}
\end{align}

Écrivons le code

Exemple de putain de code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from copy import deepcopy


fig, ax = plt.subplots()

c = 1.
x_left = 0.
x_right = 1.
t = 0.
#Signification de se disperser en 101 points
num_points = 101
dx = (x_right-x_left)/(num_points-1)
#Réservez un supplément pour les limites de Neumann
x = np.linspace(x_left, x_right+dx, num_points+1)
# u = x/Initialiser à 1000
u = 1e-3*x
#Condition aux limites de Neumann
u[-1] = u[-2]

u_before = deepcopy(u)
u_next = deepcopy(u)

#Accumuler ici la valeur de u à intervalles réguliers
u_at_t = [deepcopy(u)]


dt = 0.1 * dx/c

count = 0


while t < 0.5:
    print(np.diff(u,2))
    u_next[1:num_points] = 2.*u[1:num_points] - u_before[1:num_points] + c*c*dt*dt/(dx*dx) * np.diff(u, 2)

    #Condition aux limites de Neumann
    u_next[-1] = u_next[-2]
    u_next[0] = 0.
    u_before = deepcopy(u)
    u = deepcopy(u_next)
    u_next = np.zeros(u.shape)
    t += dt
    count += 1
    #Une fois toutes les 16 étapes
    if count == 16:
        u_at_t.append(deepcopy(u))
        count = 0


def update(u, x):
    plt.clf()
    plt.xlim(x_left, x_right)
    plt.ylim(-1e-3, 1e-3)
    plt.plot(x, u)


#Fonction d'affichage de l'animation
anim = FuncAnimation(fig, update, fargs=(
    x,), frames=u_at_t, interval=10)

plt.show()

# anim.save("foo.mp4", writer="ffmpeg",fps=10)

Notez que l'enregistrement d'une vidéo peut être une tâche lourde

Continuer

référence

Techniques de réalisation de vidéos https://qiita.com/msrks/items/e264872efa062c7d6955 Formule de différence https://qiita.com/Ushio/items/0249fd7a5363ccd914dd

Recommended Posts

Résolution d'une équation d'onde unidimensionnelle par la méthode de la différence (Python)
[Calcul scientifique / technique par Python] Résolution de l'équation de Newton unidimensionnelle par la méthode Runge-Kutta du 4ème ordre
[Calcul scientifique / technique par Python] Résolution de l'équation de Schrödinger unidimensionnelle à l'état stationnaire par méthode de tir (1), potentiel de type puits, mécanique quantique
[Calcul scientifique / technique par Python] Résolution d'une équation de Schrödinger unidimensionnelle en régime permanent par méthode de tir (2), potentiel d'oscillateur harmonique, mécanique quantique
Python - Solution numérique d'équation différentielle Méthode d'Euler et méthode de différence centrale et méthode Rungekutta
[Méthode de calcul numérique, python] Résolution d'équations différentielles ordinaires par la méthode Eular
[Python] Différence entre fonction et méthode
[Python] Différence entre la méthode de classe et la méthode statique
Le début de la méthode des éléments finis (matrice de rigidité des éléments unidimensionnels)
Algorithme d'alignement par méthode d'insertion en Python
[Calcul scientifique / technique par Python] Résolution de l'équation différentielle ordinaire du second ordre par la méthode Numerov, calcul numérique
[Calcul scientifique / technique par Python] Solution numérique d'équations d'ondes unidimensionnelles et bidimensionnelles par méthode FTCS (méthode explicite), équations aux dérivées partielles bi-courbes
Résolution d'équations de mouvement en Python (odeint)
[Python] Calcul numérique d'une équation de diffusion thermique bidimensionnelle par méthode ADI (méthode implicite de direction alternative)