Einführung

Die durch die Substanz / den Raum selbst (Medium) übertragene Schwingung, die sich im Raum ausbreitet, wird allgemein als Welle oder Wellenbewegung bezeichnet [1]. Elastische Wellen und Flüssigkeitsbewegungen entsprechen mechanischen Wellen, und elektromagnetische Wellen entsprechen elektromagnetischen Wellen. Die Gleichung, die durch die physikalische Größe $ u (x, y, z, t) $ erfüllt wird, die im Medium schwingt, ist die Wellengleichung.

** In diesem Artikel wird Python verwendet, um eine zweidimensionale Wellengleichung (bikurvierte partielle Differentialgleichung) zu erstellen ** $\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} $

Lösen ** mit der FTCS-Methode (Forward in Time and Centered Difference-Methode) [Addendum], einem Differenzschema der expliziten Methode ** (Ergänzung: ** Einschließlich des Falls einer Dimension. Siehe Inhalt (2) **)

Es ist zu beachten, dass bei der Ableitung der Dynamik dieser Gleichung die Annahme einer winzigen Schwingung und die Gleichmäßigkeit und Isotropie des Mediums, das die Welle überträgt, angenommen werden (Phasengeschwindigkeit c ist eine Konstante).

Inhalt

(1) Zweidimensionale Wellengleichung

Lösen Sie die zweidimensionale Wellengleichung mit der FTCS-Methode, einem expliziten Differenzschema. Die Berechnungsbedingungen sind wie folgt.

Berechnungsbereich: Quadrat mit Lx = 1, Ly = 1 In 41-Punkt-Schritten sowohl in x-Richtung (Nx) als auch in y-Richtung (Ny) (Nx = Ny = 41)
Der Maximalwert (t_max) der Berechnungszeit beträgt 0,5.
Zeitschritt (delta_t): Der Wert, der die numerische Stabilität berücksichtigt, wird halbiert.
Wellenphasengeschwindigkeit: c = 1
Randbedingung: Alle vier Seiten sind feste Enden (vibriert nicht und bewegt sich nicht)
Anfangsbedingung 1: Die Gaußsche Verteilungsfunktion wird zum Zeitpunkt 0 in der Mitte des Quadrats installiert. Dies "zieht" die Welle (siehe Abbildung 1).
Anfangsbedingung 2: $ \ frac {\ partielles u} {\ partielles t} = 0 \ (t = 0) $ (anfängliche Geschwindigkeit ist 0)

Zusätzlich erfordert die Berechnung Informationen über die Welle $ u $ im Schritt t = 1. In dieser Simulation wurde die Anfangsbedingung 2 unter Verwendung des durch die Mittendifferenz erhaltenen Ergebnisses bewertet [4].

Abbildung 1. Als Ausgangsbedingung wird Gauß in der Mitte des Quadrats installiert.

(2) Eindimensionale Wellengleichung

$\frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} $ Lösen Sie mit der FTCS-Methode.

Die Randbedingung ist $ u = 0 $ (festes Ende) am Endpunkt. Die Anfangsbedingung ist $ \ frac {\ partielle u} {\ partielle t} = 0 \ (t = 0) $ (Anfangsgeschwindigkeit ist 0).

Die Anfangsbedingungen für die Wellenform sind in der folgenden Abbildung dargestellt.

Abbildung 2. Anfangswellenform

Berechnungscode (1) Zweidimensionale Wellengleichung

#%matplotlib nbagg
"""
2D Wave Equation
Einheitliches Medium:Die Phasengeschwindigkeit ist zeitlich und räumlich konstant
FTCS-Methode
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation #Importieren Sie Methoden zum Erstellen von Animationen


c=1
print ("c=",c)
Lx = 1 #Länge[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 #Zeitschrittbreite. Nimm es aber klein.

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)

#Randbedingung
u[0,:,:], u[-1,:,:] , u[:,0,:] , u[:,-1,:] = 0, 0, 0, 0

#Anfangsbedingungen
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))  #Anfangsbedingungen. Installation der Gaußschen Verteilungsfunktion.

#Hado u[i,ii,1]Berechnung zu machen. Es wird durch Zentrieren der Anfangsbedingung 2 erhalten.
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])


# #Schleife zur zeitlichen Entwicklung der Wellenbewegung
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])

Der folgende Code wird verwendet, um die erhaltene Wellenfunktion u [x, y, t] zu visualisieren. Es wird in der Drahtrahmenanzeige [3] aufgezeichnet und das Ergebnis wird animiert.

#Programm zur Visualisierung

%matplotlib nbagg 
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation #Importieren Sie Methoden zum Erstellen von Animationen

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.) #Winkeleinstellung
    #ax.plot_surface(X,Y,u[X,Y,i],rstride=1,cstride=1,cmap='Greys',linewidth=0.3) #Flächendiagramm
    ax.plot_wireframe(X,Y,u[X,Y,i],rstride=1,cstride=1,color='blue',linewidth=0.3) #Drahtrahmenanzeige
    
    ax.set_title(fig_title + 'i=' + str(i))
    ax.set_zlim(-0.4,0.4) #Angeben des Zeichenbereichs von z
    ax.set_xlabel('X') #Etikette
    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")

Ergebnis (1)

Zeigt das animierte Ergebnis an. Es ist zu sehen, dass sich die an der Mittelposition erzeugte Welle herum ausbreitet.

Berechnungscode (2): Eindimensionale Wellengleichung


%matplotlib nbagg
"""
1D　Wave Equation
FTCS-Methode

"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import ArtistAnimation #Importieren Sie Methoden zum Erstellen von Animationen


T=40 #Spannung[N]
rho=0.01 #Dichte[Kg/m]
c = np.sqrt(T/rho) #Phasengeschwindigkeit

print ("c=",c)
Lx = 1 #Länge[m]
t_max = 0.1

Nx = 100

delta_x = Lx/Nx

delta_t=delta_x/c   #Maximaler Zeitschritt, der die Stabilitätsbedingung erfüllt
#print("delta_t=",delta_t)

Nt = int(t_max/delta_t)
st=delta_x/delta_t 
print("c-st=",c-st) #Stabilitätsbedingungen prüfen(Instabil wenn negativ)

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

#Randbedingung
u[0,:] = 0
u[-1,:] = 0

#Anfangsbedingungen
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]Einstellungen von
    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()

Ergebnis (2) Eindimensionale Wellengleichung

Abbildung 3. Lösung der eindimensionalen Wellengleichung

Als nächstes wird das Ergebnis durch Animation veranschaulicht.

%matplotlib nbagg 
from matplotlib.animation import ArtistAnimation #Importieren Sie Methoden zum Erstellen von Animationen

fig = plt.figure()

anim = [] #Eine Liste zum Speichern der Daten des für die Animation gezeichneten Para-Para-Diagramms

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) #Animationserstellung
plt.xlabel('x')
plt.ylabel('U')

fig.show() 
anim.save("t.gif", writer='imagemagick')   #Animation.Speichern Sie es als GIF und erstellen Sie eine GIF-Animationsdatei.

Animation der Lösung der eindimensionalen Wellengleichung

Nachtrag

(1) FTCS-Methode

Die ** FTCS-Methode (Forward in Time and Centered Difference-Methode) [2] ** in der partiellen Differentialgleichung ist ein einfaches Schema, bei dem die Vorwärtsdifferenz für die Zeitdifferenz und die zentrale Differenz für die räumliche Differenz verwendet werden. ** Einfach zu implementieren, muss jedoch genau abgestimmt werden, um numerische Instabilität zu vermeiden **

(2) Die hier behandelte Simulation berücksichtigt nicht den Dämpfungseffekt. Darüber hinaus wird die zeitliche / räumliche Abhängigkeit der Phasengeschwindigkeit nicht berücksichtigt. Referenzen [4] liefern eine rudimentäre Erklärung für diese.

Verweise

[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, ["Numerisches Rezept in der japanischen Seeversion - Rezept für die numerische Berechnung in C-Sprache"](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. = Bücher & dh = UTF8 & qid = 15039997535 & sr = 1-1 & Schlüsselwörter =% 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] Qiita-Artikel: [Wissenschaftliche / technische Berechnung von Python] Zeichnen und Visualisieren von 2D- (Farb-) Konturlinien usw., matplotlib [4] Rubin H. Landau, ["Computational Physics Application" (Übersetzung)](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 = Bücher & dh = UTF8 & qid = 15039997611 & sr = 1-2 & keywords = Landau + Computer + Physik)

[Wissenschaftlich-technische Berechnung nach Python] Numerische Lösung von 1-dimensionalen und 2-dimensionalen Wellengleichungen nach der FTCS-Methode (explizite Methode), doppelt gekrümmte partielle Differentialgleichungen