[PYTHON] Calcul numérique d'équations aux dérivées partielles avec singularité (en utilisant les équations thermiques de type Hardy-Hénon à titre d'exemple)
Préface
Ici, considérons l'équation suivante appelée équation thermique de type Hardy-Hénon.
\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*}
Cependant, $ 0 \ le \ gamma <1 $, $ \ alpha> 2 $.
La particularité de cette équation est que le terme non linéaire a une singularité à l'origine $ x = 0 $. Lorsque $ \ gamma = 0 $, il n'a pas de singularité et devient une équation thermique semi-linéaire normale.
Dans ce qui suit, une analyse numérique est effectuée dans le but d'étudier ** le comportement d'approche de la solution par rapport à une petite valeur initiale **.
Préparation
Réglez chaque paramètre comme suit:
\begin{align*}
& \alpha = 3, \\
& \gamma = 0, \, 0.25, \, 0.5, \, 0.75.
\end{align*}
La valeur initiale est la fonction suivante:
\begin{align*}
u_0(x) = \frac{1}{4} \mathrm{e}^{-10x^2}.
\end{align*}
De plus, en supposant que la solution se désintègre assez rapidement dans l'espace, nous résoudrons le problème de Diricre avec une valeur limite de $ 0 $ sur l'intervalle $ [-1, 1] $.
Après tout, pensez à résoudre numériquement l'équation suivante.
\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*}
La discrimination
Comme méthode de dispersion,Différence avant pour la différenciation temporelle,Nous avons décidé d'utiliser la différence de centre de second ordre pour le laplacien.,Le terme non linéaire utilise la valeur du point discret tel quel.pourtant,Parce que l'origine est unique,Seule l'origine doit être traitée différemment.ici,Utilisation de la moyenne intégrée uniquement autour de l'origine,Pour éviter la spécificité.Réellement,
\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*}
Est établi.
En résumé, les équations sont discrétisées comme suit:
\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*}
Résultats de la mise en œuvre et de l'analyse
Il a été implémenté en Python selon ce qui précède.Pour comparaison, une équation thermique linéaire a également été calculée en même temps.
%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")
#espace
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)
#temps
T = 3
dt = (dx ** 2) * 1/6
maxstep = int(T / dt)
#Calcul des équations différentielles
# gamma
gammas = np.linspace(0, 0.75, 4)
gammas = np.append(gammas, 0) #Pour les équations linéaires
#Valeur limite
L_BV = 0
R_BV = 0
#valeur initiale
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
#constant
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 #Pour les équations linéaires
c2 = dt * c2
a = 3 #Terme non linéaire
#Pour stocker les résultats
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" #Pour les équations linéaires
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")
Conclusion
Le comportement apocalyptique de la solution de l'équation thermique de type Hardy-Hénon par rapport à une petite valeur initiale devrait être similaire à l'équation thermique semi-linéaire et à l'équation thermique linéaire, au moins dans l'ordre inférieur.
Note) Cela ne signifie pas qu'ils sont exactement les mêmes, cela signifie qu'ils peuvent être considérés comme essentiellement les mêmes avec les modifications appropriées.
référence
・ K.Taniguchi (N.Chikami, M.Ikeda), Well-posedness and global dynamic for the critique Hardy-Hénon parabolic equation, Well-posedness Zoom Seminar.
Recommended Posts