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

Hardy_Henon_a-3.gif

Hardy_Henon_maxval_a-3.png

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.

Recommended Posts

Numerical calculation of partial differential equations with singularity (for example, asymptotic behavior analysis of Hardy-Hénon type heat equation)
Numerical calculation of differential equations with TensorFlow 2.0
Numerical analysis of ordinary differential equations with Scipy's odeint and ode
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation
Find the numerical solution of the second-order ordinary differential equation with scipy
Play with numerical calculation of magnetohydrodynamics