Lösen einer eindimensionalen Wellengleichung mit der Differenzmethode (Python)

Einfache Wellengleichung

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

Ich möchte das lösen

Randbedingung

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

Das linke Ende ist die Diricre-Randbedingung, die den Wert von $ u $ an der Grenze angibt, und das rechte Ende ist die Neumann-Randbedingung, die den Differenzwert des Grenzwerts angibt.

Die Randbedingung ist, dass das linke das feste Ende und das rechte das freie Ende ist.

Anfangsbedingungen

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

Wenn es sich um die Vibration einer Saite oder die Vibration einer vertikalen Stange handelt, beträgt die Anfangsgeschwindigkeit 0.

Diskret

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

Und legen. $ M $ ist jedoch die Anzahl der räumlichen Unterteilungen.

Ausgangsbedingung Randbedingung

Stellen Sie zunächst die ursprüngliche Lösung ein.

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

Da das Zeitdifferential erster Ordnung 0 ist, ist $ u_ {i} ^ 1 = u_ {i} ^ 0 $.

Als nächstes legen Sie die Randbedingungen fest. Das Problem ist jedoch, dass die Neumann-Randbedingung am rechten Ende durch räumliche Differenzierung ausgedrückt wird, sodass zwei Elemente verwendet werden müssen, um sie auszudrücken. Fügen Sie daher der Grenze ganz rechts ein weiteres falsches Element hinzu. Mit anderen Worten

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

Wenn ja, ist die Neumann-Randbedingung am rechten Ende, dass die räumliche Differenzierung erster Ordnung 0 ist, erfüllt. Dann kann die Randbedingung ganz links in der Richtung wie folgt ausgedrückt werden.

u_1^n = 0

Schließlich wird die Gleichung selbst diskretisiert.

Diskriminierung der räumlichen Differenzierung

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

Dies kann durch maßgeschneiderte Erweiterung von $ u_ {i + 1} $ und $ u_ {i-1} $ mit $ x = x_i $ abgeleitet werden.

Trennung der Zeitentwicklung

Denken wie räumliche Differenzierung

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

Denn es ist

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

Schreiben wir den Code

Beispiel Fickcode

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.
#Das heißt, an 101 Punkten zu zerstreuen
num_points = 101
dx = (x_right-x_left)/(num_points-1)
#Reservieren Sie ein Extra für Neumann-Grenzen
x = np.linspace(x_left, x_right+dx, num_points+1)
# u = x/Bei 1000 initialisieren
u = 1e-3*x
#Neumann-Randbedingung
u[-1] = u[-2]

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

#Akkumulieren Sie hier in regelmäßigen Abständen den Wert von u
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)

    #Neumann-Randbedingung
    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
    #Einmal alle 16 Schritte
    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)


#Funktion zur Anzeige von Animationen
anim = FuncAnimation(fig, update, fargs=(
    x,), frames=u_at_t, interval=10)

plt.show()

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

Beachten Sie, dass das Speichern eines Videos eine schwere Aufgabe sein kann

Fortsetzen

Referenz

Techniken zum Erstellen von Videos https://qiita.com/msrks/items/e264872efa062c7d6955 Differenzformel https://qiita.com/Ushio/items/0249fd7a5363ccd914dd

Recommended Posts

Lösen einer eindimensionalen Wellengleichung mit der Differenzmethode (Python)
[Wissenschaftlich-technische Berechnung nach Python] Lösen der eindimensionalen Newton-Gleichung nach der Runge-Kutta-Methode 4. Ordnung
[Wissenschaftlich-technische Berechnung von Python] Lösen der eindimensionalen Schrödinger-Gleichung im stationären Zustand durch Schießmethode (1), Potential vom Well-Typ, Quantenmechanik
[Wissenschaftlich-technische Berechnung von Python] Lösen der eindimensionalen Schrödinger-Gleichung im stationären Zustand durch Aufnahmemethode (2), harmonisches Oszillatorpotential, Quantenmechanik
Python - Differentialgleichung Numerische Lösung Euler-Methode & Zentrale Differenzmethode & Rungekutta-Methode
[Numerische Berechnungsmethode, Python] Lösen gewöhnlicher Differentialgleichungen mit der Eular-Methode
[Python] Unterschied zwischen Funktion und Methode
[Python] Unterschied zwischen Klassenmethode und statischer Methode
Der Beginn der Finite-Elemente-Methode (eindimensionale Elementsteifigkeitsmatrix)
Ausrichtungsalgorithmus durch Einfügemethode in Python
[Wissenschaftlich-technische Berechnung mit Python] Lösen der gewöhnlichen Differentialgleichung zweiter Ordnung nach der Numerov-Methode, numerische Berechnung
[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
Lösen von Bewegungsgleichungen in Python (odeint)
[Python] Numerische Berechnung der zweidimensionalen Wärmediffusionsgleichung mit der ADI-Methode (implizite Methode mit alternativer Richtung)