[PYTHON] Numerical calculation of partial differential equations with singularity (for example, asymptotic behavior analysis of Hardy-Hénon type heat equation)

Preface

Here, consider the following equation called the Hardy-Hénon type heat equation.

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

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

The feature of this equation is that the nonlinear term has singularity at the origin $ x = 0 $. When $ \ gamma = 0 $, it has no singularity and becomes a normal semi-linear heat equation.

In the following, numerical analysis is performed for the purpose of investigating ** the asymptotic behavior of the solution for a small initial value **.

Preparation

Set each parameter as follows:

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

The initial value is the following function:

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

Furthermore, assuming that the solution decays fast enough in space, we solve the Dirichlet problem with a boundary value of $ 0 $ on the interval $ [-1, 1] $.

After all, consider solving the following equation numerically.

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

Discretization

As a discretization method,Forward difference for time derivative,We decided to use the second-order center difference for the Laplacian.,Non-linear term uses the value of the discrete point as it is.However,Because the origin is unique,Only the origin needs to be treated differently.here,Using the integral mean only around the origin,To avoid specificity.Actually, 0 \le \gamma<1Than, |x|^{- \gamma}Is integrable around the origin,

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

Is established.

In summary, the equations are discretized as follows:

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

Implementation and analysis results

Implemented in Python according to the above. A linear heat equation was also calculated at the same time for comparison.

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

#space
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)

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

#Calculation of differential equations

# gamma
gammas = np.linspace(0, 0.75, 4)
gammas = np.append(gammas, 0) #For linear equations

#Boundary value
L_BV = 0
R_BV = 0

#initial value
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 #For linear equations
c2 = dt * c2
a = 3 #Non-linear term exponentiation

#For storing results
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" #For linear equations
    
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

The apocalyptic behavior of the solution of the Hardy-Hénon type heat equation with respect to a small initial value is expected to be similar to the semi-linear heat equation and the linear heat equation, at least in the lower order.

Note) It does not mean that they are exactly the same. It means that they can be regarded as essentially the same with appropriate modifications.

reference

・ K.Taniguchi (N.Chikami, M.Ikeda), Well-posedness and global dynamics for the critical Hardy-Hénon parabolic equation, Well-posedness Zoom Seminar.