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, 0 \le \gamma<1Que, |x|^{- \gamma}Est intégrable autour de l'origine,

\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")

Hardy_Henon_a-3.gif

Hardy_Henon_maxval_a-3.png

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

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)
L'histoire du calcul numérique des équations différentielles avec TensorFlow 2.0
Analyse numérique des équations différentielles ordinaires avec l'odeint et l'ode de Scipy
[Calcul scientifique / technique par Python] Solution numérique de l'équation de Laplace-Poisson bidimensionnelle pour la position électrostatique par la méthode Jacobi, équation aux dérivées partielles elliptiques, problème des valeurs aux limites
[Calcul scientifique / technique par Python] Solution numérique d'une équation de conduction thermique non stationnaire unidimensionnelle par méthode Crank-Nicholson (méthode implicite) et méthode FTCS (méthode de solution positive)
Trouvez la solution numérique de l'équation différentielle ordinaire du second ordre avec scipy