[PYTHON] Numerische Berechnung partieller Differentialgleichungen mit Singularität (am Beispiel der thermischen Gleichungen vom Hardy-Hénon-Typ)

Vorwort

Betrachten Sie hier die folgende Gleichung, die als thermische Gleichung vom Hardy-Hénon-Typ bezeichnet wird.

\begin{align*}
\left \{
\begin{aligned}
& \partial_t u - \Delta u = \frac{|u|^{\alpha - 2}u}{|x|^{\gamma}},
\quad (t,x) \in \mathbf{R}_+ \times \mathbf{R}, \\
& u(0,x)=u_0(x).
\end{aligned}
\right.
\end{align*}

$ 0 \ le \ gamma <1 $, $ \ alpha> 2 $.

Das Merkmal dieser Gleichung ist, dass der nichtlineare Term am Ursprung $ x = 0 $ Singularität aufweist. Wenn $ \ gamma = 0 $, hat er keine Singularität und wird zu einer normalen halblinearen thermischen Gleichung.

Im Folgenden wird eine numerische Analyse durchgeführt, um ** das Annäherungsverhalten der Lösung in Bezug auf einen kleinen Anfangswert ** zu untersuchen.

Vorbereitung

Stellen Sie jeden Parameter wie folgt ein:

\begin{align*}
& \alpha = 3, \\
& \gamma = 0, \, 0.25, \, 0.5, \, 0.75.
\end{align*}

Der Anfangswert ist die folgende Funktion:

\begin{align*}
u_0(x) = \frac{1}{4} \mathrm{e}^{-10x^2}.
\end{align*}

Unter der Annahme, dass die Lösung im Raum schnell genug zerfällt, lösen wir das Diricre-Problem mit einem Grenzwert von $ 0 $ für das Intervall $ [-1, 1] $.

Schließlich sollten Sie die folgende Gleichung numerisch lösen.

\begin{align*}
\left \{
\begin{aligned}
& \partial_t u - \Delta u = \frac{|u|u}{|x|^{\gamma}},
\quad (t,x) \in \mathbf{R}_+ \times (-1, 1), \\
& u(0,x)=\frac{1}{4} \mathrm{e}^{-10x^2}, \\
&u(t,-1)=u(t,1)=0.
\end{aligned}
\right.
\end{align*}

Diskriminierung

Als Ausbreitungsmethode,Vorwärtsdifferenz zur Zeitdifferenzierung,Wir haben uns entschieden, die Mittendifferenz zweiter Ordnung für Laplace zu verwenden.,Der nichtlineare Term verwendet den Wert des diskreten Punkts so wie er ist.jedoch,Weil der Ursprung einzigartig ist,Nur der Ursprung muss anders behandelt werden.Hier,Verwenden der integrierten Mittelwertbildung nur um den Ursprung,Um Spezifität zu vermeiden.Tatsächlich, 0 \le \gamma<1Als, |x|^{- \gamma}Ist um den Ursprung integrierbar,

\begin{align*}
\frac{1}{\Delta x} \int^{\Delta x / 2}_{- \Delta x / 2} |x|^{- \gamma} dx
= 2 \times \frac{1}{\Delta x} \int^{\Delta x / 2}_{0} |x|^{- \gamma} dx
= \frac{2^{\gamma}}{1 - \gamma} {\Delta x}^{- \gamma}
\end{align*}

Ist festgelegt.

Zusammenfassend sind die Gleichungen wie folgt diskretisiert:

\begin{align*}
u^{n+1}_{i} = 
\left \{
\begin{aligned}
& u^{n}_{i} + 
\frac{\Delta t}{{\Delta x}^2} (u^{n}_{i+1} - 2 u^{n}_{i} + u^{n}_{i-1})
+ \Delta t \frac{|u^{n}_{i}|u^{n}_{i}}{|x_{i}|^{\gamma}}, \quad (x_{i} \neq 0), \\
& u^{n}_{i} + 
\frac{\Delta t}{{\Delta x}^2} (u^{n}_{i+1} - 2 u^{n}_{i} + u^{n}_{i-1}) 
+ \Delta t \frac{2^{\gamma}}{1 - \gamma} \frac{|u^{n}_{i}|u^{n}_{i}}{{\Delta x}^{\gamma}}, \quad (x_{i} = 0).
\end{aligned}
\right.
\end{align*}

Implementierungs- und Analyseergebnisse

Es wurde in Python wie oben implementiert implementiert. Zum Vergleich wurde gleichzeitig auch eine lineare thermische Gleichung berechnet.

%matplotlib nbagg
import numpy as np
import matplotlib as mpl
import matplotlib.animation as animation
import matplotlib.pyplot as plt
from numpy import pi, sin, cos, exp, real, imag
from scipy.fftpack import fft, ifft, fftfreq
import warnings
warnings.filterwarnings("ignore")

#Raum
L_end = -1
R_end = 1
N = 2**8 + 1
x = np.linspace(L_end, R_end, N)
dx = (R_end - L_end)/N
origin_idx = int((len(x)-1)/2)

#Zeit
T = 3
dt = (dx ** 2) * 1/6
maxstep = int(T / dt)

#Berechnung von Differentialgleichungen

# gamma
gammas = np.linspace(0, 0.75, 4)
gammas = np.append(gammas, 0) #Für lineare Gleichungen

#Grenzwert
L_BV = 0
R_BV = 0

#Ursprünglicher Wert
alpha = 1/4
beta = 10
u = alpha * exp(-beta * (x**2))
u = np.tile(u, (len(gammas), 1))
u[:, 0] = L_BV
u[:, -1] = R_BV

#Konstante
c1 = dt/(dx**2)
c2 = np.zeros_like(u)
for i, gamma in enumerate(gammas):
    c2[i, 0:origin_idx] = abs(x[0:origin_idx]) ** (-gamma)
    c2[i, origin_idx+1:] = abs(x[origin_idx+1:]) ** (-gamma)
    c2[i, origin_idx] = (2**gamma) / (1 - gamma) * (dx**(-gamma))
c2[-1] = 0 #Für lineare Gleichungen
c2 = dt * c2
a = 3 #Nichtlinearer Term

#Zum Speichern von Ergebnissen
fig = plt.figure(figsize=(6,3))
fig.subplots_adjust(bottom=0.2)
ims = []
cmap = plt.get_cmap("tab10")
max_vals = []
labels = []
for gamma in gammas:
    labels.append("gamma: " + str(gamma))
labels[-1] = "linear" #Für lineare Gleichungen
    
for n in range(maxstep):
    u[:, 1:-1] += c1 * (u[:, 2:] - 2*u[:, 1:-1] + u[:, 0:-2]) + c2[:, 1:-1] * abs(u[:, 1:-1])**(a-2) * u[:, 1:-1]
    u[:, 0] = L_BV
    u[:, -1] = R_BV
    
    if n % 1000 < 1e-10: 
        im = []
        for i in range(len(gammas)):
            im += plt.plot(x, u[i,:], color=cmap(i))
        ims.append(im)
        max_val = []
        for i in range(len(gammas)):
            max_val.append(np.max(u[i, :]))
        max_vals.append(max_val)

plt.title("Solutions(a=3)")
plt.xlabel("x")
plt.ylabel("u")
plt.legend(labels)
ani = animation.ArtistAnimation(fig, ims, interval=1)

fig2 = plt.figure(figsize=(6,3))
fig2.subplots_adjust(bottom=0.2)
for i in range(len(gammas)):
    plt.plot(dt * 1000 * np.array(range(len(max_vals))), np.array(max_vals)[:, i], color=cmap(i), label = labels[i])
plt.title("Max Values(a=3)")
plt.xlabel("t")
plt.ylabel("max(u)")
plt.legend()

plt.show()
ani.save(f"results\Hardy_Henon_a-3.gif", savefig_kwargs={"bbox_inches":"tight"})
plt.savefig(f"results\Hardy_Henon_maxval_a-3.png ", bbox_inches="tight")

Fazit

Es wird erwartet, dass das apokalyptische Verhalten der Lösung der thermischen Gleichung vom Hardy-Hénon-Typ in Bezug auf einen kleinen Anfangswert der semilinearen thermischen Gleichung und der linearen thermischen Gleichung zumindest in der unteren Ordnung ähnlich ist.

Hinweis) Dies bedeutet nicht, dass sie genau gleich sind. Dies bedeutet, dass sie mit geeigneten Modifikationen als im Wesentlichen gleich angesehen werden können.

Referenz

・ K. Taniguchi (N.Chikami, M.Ikeda), Gutmütigkeit und globale Dynamik für die kritische Hardy-Hénon-Parabolgleichung, Gutmütigkeit-Zoomseminar.