introduction

La vibration transmise à travers la substance / l'espace lui-même (le milieu) qui se propage dans l'espace est généralement appelée mouvement d'onde ou d'onde [1]. Les ondes élastiques et le mouvement des fluides correspondent aux ondes mécaniques et les ondes électromagnétiques correspondent aux ondes électromagnétiques. L'équation satisfaite par la quantité physique $ u (x, y, z, t) $ qui vibre dans le milieu est l'équation d'onde.

** Dans cet article, Python est utilisé pour créer une équation d'onde bidimensionnelle (équation différentielle partielle bicurvée) ** $\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} $

Résolution ** par la méthode FTCS (Forward in Time and Centered Difference method) [Addendum], qui est un schéma de différence de la méthode explicite **

Notez qu'en dérivant la dynamique de cette équation, on suppose l'hypothèse d'une vibration infime et l'uniformité et la nature isotrope du milieu qui transmet l'onde (la vitesse de phase c est une constante).

Contenu

(1) Équation d'onde bidimensionnelle

Résolvez l'équation d'onde bidimensionnelle par la méthode FTCS, qui est un schéma de différence explicite. Les conditions de calcul sont les suivantes.

Zone de calcul: Carré avec Lx = 1, Ly = 1 Par incréments de 41 points dans la direction x (Nx) et la direction y (Ny) (Nx = Ny = 41)
La valeur maximale (t_max) du temps de calcul est de 0,5.
Pas de temps (delta_t): La valeur qui prend en compte la stabilité numérique est divisée par deux.
Vitesse de phase d'onde: c = 1
Condition aux limites: les quatre côtés sont des extrémités fixes (ne vibre pas et ne bouge pas)
Condition initiale 1: la fonction de distribution gaussienne est installée au centre du carré au temps 0. Cela «tire» la vague (voir la figure 1).
Condition initiale 2: $ \ frac {\ partial u} {\ partial t} = 0 \ (t = 0) $ (la vitesse initiale est 0)

De plus, le calcul nécessite des informations sur la vague $ u $ dans le pas t = 1. Dans cette simulation, la condition initiale 2 a été évaluée en utilisant le résultat obtenu par différence de centre [4].

Figure 1. Comme condition initiale, la gaussienne est installée au centre du carré.

(2) Équation d'onde unidimensionnelle

$\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} $ Résolvez avec la méthode FTCS.

La condition aux limites est $ u = 0 $ (extrémité fixe) au point final. La condition initiale est $ \ frac {\ partial u} {\ partial t} = 0 \ (t = 0) $ (la vitesse initiale est 0).

Les conditions initiales de la forme d'onde sont indiquées dans la figure ci-dessous.

Figure 2. Forme d'onde initiale

Code de calcul (1) Équation d'onde bidimensionnelle

#%matplotlib nbagg
"""
2D Wave Equation
Milieu uniforme:La vitesse de phase est constante dans le temps et dans l'espace
Méthode FTCS
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation #Importer des méthodes pour créer des animations


c=1
print ("c=",c)
Lx = 1 #longueur[m]
Ly = 1
t_max = 0.5

Nx = 41
Ny =41

delta_x = Lx/Nx
delta_y = Ly/Ny
delta_t=(1/(1/delta_x+1/delta_y))/c*0.5 #Largeur du pas de temps. Prenez-le petit, cependant.

print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
 

print("Nt=",Nt)

stx=delta_x/delta_t
sty=delta_y/delta_t



u =  np.empty((Nx,Ny,Nt),dtype=float)

#condition limite
u[0,:,:], u[-1,:,:] , u[:,0,:] , u[:,-1,:] = 0, 0, 0, 0

#Conditions initiales
for i in range (Nx) :
    for ii in range(Ny):
        u[i,ii,0]  = (np.exp(-((i*delta_x-Lx/2)**2+(ii*delta_y-Ly/2)**2)/0.01))  #Conditions initiales. Installation de la fonction de distribution gaussienne.

#Hado u[i,ii,1]Calcul à faire. Il est obtenu en centrant la condition initiale 2.
for i in range(1,Nx-1): 
    for ii in range(1,Ny-1):
        u[i,ii,1]= u[i,ii,0]+0.5*((delta_t/delta_x)**2)*(c**2)*(u[i+1,ii, 0]+u[i-1, ii,0]-2*u[i,ii,0])\
+0.5*((delta_t/delta_y)**2)*(c**2)*(u[i,ii+1, 0]+u[i,ii-1, 0]-2*u[i,ii,0])


# #Boucle pour le développement temporel du mouvement des vagues
for j in range(1,Nt-1):
    for i in range(1,Nx-1):
        for ii in range(1,Ny-1):
            u[i,ii, j+1] = 2*u[i,ii,j]-u[i,ii,j-1]+((c/stx)**2)* (u[i+1,ii,j]+u[i-1,ii,j]\
-2*u[i,ii,j])+((c/sty)**2)* (u[i,ii+1,j]+u[i,ii-1,j]-2*u[i,ii,j])

Le code suivant est utilisé pour visualiser la fonction d'onde obtenue u [x, y, t]. Il est tracé dans l'affichage filaire [3] et le résultat est animé.

#Programme de visualisation

%matplotlib nbagg 
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation #Importer des méthodes pour créer des animations

fig = plt.figure(figsize = (8, 6))
ax = fig.gca(projection='3d')

x=list(range(Nx))
y=list(range(Ny))
X, Y = np.meshgrid(x,y)

def update(i,u,fig_title):
    if i != 0:
        ax.clear()
    ax.view_init(elev=60., azim=60.) #Réglage de l'angle
    #ax.plot_surface(X,Y,u[X,Y,i],rstride=1,cstride=1,cmap='Greys',linewidth=0.3) #Tracé de surface
    ax.plot_wireframe(X,Y,u[X,Y,i],rstride=1,cstride=1,color='blue',linewidth=0.3) #Affichage de trame de fil
    
    ax.set_title(fig_title + 'i=' + str(i))
    ax.set_zlim(-0.4,0.4) #Spécification de la plage de dessin de z
    ax.set_xlabel('X') #étiquette
    ax.set_ylabel('Y')
    ax.set_zlabel('U')


    
    return ax


ani = animation.FuncAnimation(fig,update,fargs = (u,'Wave motion: time step='), interval =1, frames = Nt,blit=False)
fig.show()

ani.save("output.gif", writer="imagemagick")

Résultat (1)

Affiche le résultat animé. On peut voir que l'onde générée en position centrale se propage autour.

Code de calcul (2): équation d'onde unidimensionnelle


%matplotlib nbagg
"""
1D　Wave Equation
Méthode FTCS

"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation #Importer des méthodes pour créer des animations


T=40 #tension[N]
rho=0.01 #densité[Kg/m]
c = np.sqrt(T/rho) #Vitesse de phase

print ("c=",c)
Lx = 1 #longueur[m]
t_max = 0.1

Nx = 100

delta_x = Lx/Nx

delta_t=delta_x/c   #Pas de temps maximum qui satisfait la condition de stabilité
#print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
st=delta_x/delta_t 
print("c-st=",c-st) #Vérification des conditions de stabilité(Instable si négatif)

u =  np.zeros([Nx,Nt])

#condition limite
u[0,:] = 0
u[-1,:] = 0

#Conditions initiales
for i in range (Nx) :
    if i <= 0.8*Nx :
        u[i,0]=1.25*i/Nx
    
    else :  
        u[i,0] = 5.0*(1-i/Nx)

for i in range(1,Nx-1): # u[:,1]paramètres de
    u[i,1]  = u[i,0]+0.5*((delta_t/delta_x)**2)*(c**2)*(u[i+1,0]+u[i-1,0]-2*u[i,0])

for j in range(Nt-1):
    for i in range(1,Nx-1):
        u[i,j+1] = 2*u[i,j]-u[i,j-1]+(c**2/st**2)* (u[i+1,j]+u[i-1,j]-2*u[i,j])
        

x=list(range(Nx))
y=list(range(Nt))

X, Y = np.meshgrid(x,y)

def functz(u):
    z=u[X,Y]
    return z

Z = functz(u)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z, color='r')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('U')

plt.show()

Résultat (2) Équation d'onde unidimensionnelle

Figure 3. Solution d'une équation d'onde unidimensionnelle

Ensuite, le résultat est illustré par une animation.

%matplotlib nbagg 
from matplotlib.animation import ArtistAnimation #Importer des méthodes pour créer des animations

fig = plt.figure()

anim = [] #Une liste pour stocker les données du diagramme para-para dessiné pour l'animation

for i in range(Nt):
    U=list(u[:,i])
    x=list(range(Nx))
    if i % int(Nt*0.01) ==0: 
        im=plt.plot(x,U, '-', color='red',markersize=10, linewidth = 2, aa=True)
        anim.append(im)

anim = ArtistAnimation(fig, anim) #Création d'animation
plt.xlabel('x')
plt.ylabel('U')

fig.show() 
anim.save("t.gif", writer='imagemagick')   #Animation.Enregistrez-le en tant que gif et créez un fichier d'animation gif.

Animation de la solution de l'équation d'onde unidimensionnelle

Addenda

(1) Méthode FTCS

La ** méthode FTCS (méthode Forward in Time and Centered Difference) [2] ** dans l'équation différentielle partielle est un schéma simple qui prend la différence directe pour le différentiel temporel et la différence centrale pour le différentiel spatial. ** Facile à mettre en œuvre, mais doit être minuté pour éviter l'instabilité numérique **

(2) La simulation traitée ici ne prend pas en compte l'effet d'amortissement. De plus, la dépendance temps / espace de la vitesse de phase n'est pas prise en compte. Les références [4] en fournissent une explication rudimentaire.

Les références

[1] Yosuke Nagaoka, ["Vibration and Waves"](https://www.amazon.co.jp/%E6%8C%AF%E5%8B%95%E3%81%A8%E6%B3%A2 -% E9% 95% B7% E5% B2% A1-% E6% B4% 8B% E4% BB% 8B / dp / 4785320451) [2] William H. Press, ["Recette numérique en version japonaise de la mer-Recette pour le calcul numérique en langage C"](https://www.amazon.co.jp/%E3%83%8B% E3% 83% A5% E3% 83% BC% E3% 83% A1% E3% 83% AA% E3% 82% AB% E3% 83% AB% E3% 83% AC% E3% 82% B7% E3% 83% 94-% E3% 82% A4% E3% 83% B3-% E3% 82% B7% E3% 83% BC-% E6% 97% A5% E6% 9C% AC% E8% AA% 9E% E7 % 89% 88-C% E8% A8% 80% E8% AA% 9E% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E6% 95% B0% E5% 80% A4% E8 % A8% 88% E7% AE% 97% E3% 81% AE% E3% 83% AC% E3% 82% B7% E3% 83% 94-Press-William-H / dp / 4874085601 / ref = sr_1_1? S = livres & ie = UTF8 & qid = 15039997535 & sr = 1-1 & mots-clés =% E3% 83% 8B% E3% 83% A5% E3% 83% BC% E3% 83% A1% E3% 83% AA% E3% 82% AB% E3% 83% AB% E3% 83% AC% E3% 82% B7% E3% 83% 94) [3] Article Qiita: [Calcul scientifique / technique par Python] Dessin et visualisation de lignes de contour 2D (couleur), etc., matplotlib [4] Rubin H. Landau, ["Computational Physics Application" (traduction)](https://www.amazon.co.jp/%E8%A8%88%E7%AE%97%E7%89% A9% E7% 90% 86% E5% AD% A6-% E5% BF% 9C% E7% 94% A8% E7% B7% A8-Landau-Rubin-H / dp / 4254130872 / ref = sr_1_2? S = livres & ie = UTF8 & qid = 15039997611 & sr = 1-2 & keywords = Landau + calcul + physique)

[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