Open AccessArticle

A Note on a Meshless Method for Fractional Laplacian at Arbitrary Irregular Meshes

UNED, ETSII, 28040 Madrid, Spain

Instituto de Matemática Interdisciplinar, Departamento de Análisis Matemático y Matemática Aplicada, UCM, 28040 Madrid, Spain

Departamento de Análisis Matemático y Matemática Aplicada, UCM, 28040 Madrid, Spain

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Mathematics 2021, 9(22), 2843; https://doi.org/10.3390/math9222843

Submission received: 28 September 2021 / Revised: 3 November 2021 / Accepted: 8 November 2021 / Published: 10 November 2021

(This article belongs to the Special Issue Applications of Partial Differential Equations in Engineering)

Download

Browse Figures

Figure 1
Clouds of nodes with irregular distribution. "> Figure 2
Solutions given by the method for <math display="inline"><semantics> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> </semantics></math> (left) <math display="inline"><semantics> <mrow> <mi>m</mi> <mo>=</mo> <mn>3</mn> </mrow> </semantics></math> (centre) and <math display="inline"><semantics> <mrow> <mi>m</mi> <mo>=</mo> <mn>5</mn> </mrow> </semantics></math> (right) using the first cloud of <a href="#mathematics-09-02843-f001" class="html-fig">Figure 1</a>. "> Figure 3
Solutions given by the method for <math display="inline"><semantics> <mrow> <mi>m</mi> <mo>=</mo> <mn>2</mn> </mrow> </semantics></math> (left) <math display="inline"><semantics> <mrow> <mi>m</mi> <mo>=</mo> <mn>4</mn> </mrow> </semantics></math> (centre) and <math display="inline"><semantics> <mrow> <mi>m</mi> <mo>=</mo> <mn>10</mn> </mrow> </semantics></math> (right) using the second cloud of <a href="#mathematics-09-02843-f001" class="html-fig">Figure 1</a>. "> Figure 4
Solutions given by the method for <math display="inline"><semantics> <mrow> <mi>t</mi> <mo>=</mo> <mn>3</mn> <mo>,</mo> <mn>5</mn> </mrow> </semantics></math> and 10 in the second example. ">

Versions Notes

Abstract

The existence and uniqueness of the discrete solutions of a porous medium equation with diffusion are demonstrated. The Cauchy problem contains a fractional Laplacian and it is equivalent to the extension formulation in the sense of trace and harmonic extension operators. By using the generalized finite difference method, we obtain the convergence of the numerical solution to the classical/theoretical solution of the equation for nonnegative initial data sufficiently smooth and bounded. This procedure allows us to use meshes with complicated geometry (more realistic) or with an irregular distribution of nodes (providing more accurate solutions where needed). Some numerical results are presented in arbitrary irregular meshes to illustrate the potential of the method.

Keywords:

fractional Laplacian; generalized finite difference method; discrete maximum principle; convergence

1. Introduction

The well-known porous medium equation appears in the description of different phenomena. There are a number of physical applications where this simple model appear in a natural way, mainly to describe processes involving fluid flow, heat transfer or diffusion. Moreover, applications have been found in mathematical biology, spread of viscous fluids, boundary layer theory and other fields (see [1]). In this paper, we consider the porous medium equation with fractional Laplacian, which represents processes with anomalous diffusion. This kind of diffusion is observed in continuum mechanics, phase transition phenomena, population dynamics, optimal control, image processing, game theory, finance and others (see [2] and the references therein).

This paper deals with a numerical scheme given by the generalized finite difference method (GFDM) in order to solve the one dimensional porous medium equation with fractional diffusion

\{\begin{matrix} w_{t} + {(- Δ)}^{1 / 2} w^{m} = 0, & x \in R^{N}, t \in (0, \infty), \\ w (x, 0) = f (x), & x \in R^{N}, \end{matrix}

(1)

for

m, N \geq 1

and nonnegative initial data in

L^{1} (R^{N}) \cap L^{\infty} (R^{N})

. The existence, uniqueness and regularity of solutions to (1) are well known [3]. In order to remove the nonlocality induced by the fractional Laplacian

{(- Δ)}^{1 / 2}

, Caffarelli and Silvestre showed that any power of the fractional Laplacian in

R^{N}

could be realized as an operator that mapped a Dirichlet boundary condition to a Neumann-type condition via an extension on the upper half-space

R^{N} \times R^{+}

. In [3], the authors stated that (1) is equivalent to the so-called extension formulation (see [4,5]):

\{\begin{matrix} Δ u (x, y, t) = 0, & x \in R^{N}, y > 0, t > 0, \\ \frac{\partial u^{1 / m}}{\partial t} (x, 0, t) = \frac{\partial u}{\partial y} (x, 0, t), & x \in R^{N}, y = 0, t > 0, \\ u (x, 0, 0) = f^{m} (x), & x \in R^{N} . \end{matrix}

(2)

The equivalence of both problems holds in the sense of trace and harmonic extension operators, that is,

w (x, t) = T r (u^{1 / m} (x, y, t)), u (x, y, t) = E (w^{m} (x, t)) .

See [3,6,7] for more details about the Silvestre–Caffarelli extension. The reader can also see [8,9] for more information about the method.

Numerical results for the discrete problem have been obtained using the finite difference method in [6]. In [7], authors used the standart finite difference method to solve a Caputo-type parabolic equation with a fractional laplacian. Furthermore, recently, Chen and Shen have solved numerically Poisson-type problems with the diagonalization matrix method and the enriched spectral method [10]. Padgett dealt with the same problem numerically in [11]. In addition, different numerical methods have been applied for approximating the fractional powers of the laplacian operator, see [12,13,14].

We approximate the solutions to (2) in the whole space by the solutions of the problem posed a bounded domain. In the sequel, we study the problem with

N = 1

, because the numerous indices in the formulas and to further simplify the notation but all the arguments are also valid for

N > 1

. Consider the problem in a bounded domain

Ω

with boundary

Γ = \partial Ω

The generalized finite difference method is a meshless procedure based on the Taylor series and moving least squares (see [15]). The main advantage of the meshless methods is the possibility of using completely irregular meshes. In this way, one can concentrate a higher number of nodes in certain parts of the domain and use the previous information of the problem. In this paper, we consider the porous medium phenomena, which is of a diffusive nature; therefore, the known information about the problem such as the Barenblatt profiles [16] can be used for placing a higher number of nodes in the centre of the profiles and less in the tails. Despite this fact, the explicit formulae remain in a simple form, even for nonlinear equations. Several applications of the GFDM have been implemented in the recent years such as modelling of chemotaxis models [17] and recently, for solving fractional differential equations [18]. For completeness, let us state the basics of the method, which can be found in the above cited references. Let

M = {x_{i} = (x_{i}, y_{i})}_{i = 1}^{P}

be a set of data points in the computational domain. For each

x_{i} \in M

, call

u (x_{i}, y_{i}, t_{j}) = u_{i}^{j}

the solution of (2) in the point of the mesh and

u (x_{i}, y_{i}, t_{j}) \approx U_{i}^{j}

the solution of the numerical method. Let

S_{0}

be a localized star centered at

x_{0} \in M

with s neighbors. Call

h_{i} = x_{i} - x_{0}

and

k_{i} = y_{i} - y_{0}

. Then, the Taylor’s expansions around

x_{i}

each data point of the

S_{0}

star are

u (x_{i}) = u (x_{0}) + h_{i} \frac{\partial u}{\partial x} |_{x_{0}} + k_{i} \frac{\partial u}{\partial y} |_{x_{0}} + \frac{1}{2} {[h_{i}^{2} \frac{\partial^{2} u}{\partial x^{2}}|}_{x_{0}} + 2 h_{i} k_{i} \frac{\partial^{2} u}{\partial x \partial y} |_{x_{0}} + k_{i}^{2} \frac{\partial^{2} u}{\partial y^{2}} |_{x_{0}}] + R^{i} .

Here,

R^{i}

is the remainder term. Consider the weighted sum of squares of these remainder terms based on Moving Least Squares (MLS) and then the residual function is defined as

R e (u) = \sum_{i = 1}^{s} {[{(u (x_{0}) - u (x_{i}) + h_{i} \frac{\partial u}{\partial x}|}_{x_{0}} + k_{i} \frac{\partial u}{\partial y}|}_{x_{0}} + \frac{1}{2} {[h_{i}^{2} \frac{\partial^{2} u}{\partial x^{2}}|}_{x_{0}} + 2 h_{i} k_{i} \frac{\partial^{2} u}{\partial x \partial y} |_{x_{0}} + k_{i}^{2} \frac{\partial^{2} u}{\partial y^{2}} |_{x_{0}}) δ_{0, i}]^{2},

(3)

where

δ_{i} : = δ (∥ x_{0} - x_{i} ∥)

is some weighting function defined in MLS theory (see [19]). Examples of

δ

can be found in [15]. By minimizing expression (3) with respect to the partial derivatives, the linear system

AD = b

is obtained (see [15]), where

A = (\begin{matrix} h_{1} & h_{2} & \dots & h_{s} \\ k_{1} & k_{2} & \dots & k_{s} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ h_{1} k_{1} & h_{2} k_{2} & \dots & h_{s} k_{s} \end{matrix}) (\begin{matrix} δ_{1}^{2} \\ δ_{2}^{2} \\ \dots \\ δ_{s}^{2} \end{matrix}) (\begin{matrix} h_{1} & k_{1} & \dots & h_{1} k_{1} \\ h_{2} & k_{2} & \dots & h_{2} k_{2} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ h_{s} & k_{s} & \dots & h_{s} k_{s} \end{matrix}),

D^{T} = (\frac{\partial u_{0}}{\partial x}, \frac{\partial u_{0}}{\partial y}, \frac{\partial^{2} u_{0}}{\partial x^{2}}, \frac{\partial^{2} u_{0}}{\partial y^{2}}, \frac{\partial^{2} u_{0}}{\partial x \partial y}) .

and

\begin{matrix} b^{T} = & (\sum_{i = 1}^{s} (- u_{0} + u_{i}) h_{i} δ_{i}^{2}, \sum_{i = 1}^{s} (- u_{0} + u_{i}) k_{i} δ_{i}^{2}, \sum_{i = 1}^{s} (- u_{0} + u_{i}) \frac{h_{i}^{2} δ_{i}^{2}}{2}, \\ \sum_{i = 1}^{s} (- u_{0} + u_{i}) \frac{k_{i}^{2} δ_{i}^{2}}{2}, \sum_{i = 1}^{s} (- u_{0} + u_{i}) h_{i} k_{i} δ_{i}^{2}) . \end{matrix}

By solving the above system, the GFD formulae are deduced

\begin{matrix} \frac{\partial u}{\partial x} |_{x_{0}} = - λ_{01} & u (x_{0}) + \sum_{i = 1}^{s} λ_{i 1} u (x_{i}), \frac{\partial u}{\partial y} |_{x_{0}} = - λ_{02} u (x_{0}) + \sum_{i = 1}^{s} λ_{i 2} u (x_{i}), \\ Δ u (x_{0}) = - λ_{00} u (x_{0}) + \sum_{i = 1}^{s} λ_{i 0} u (x_{i}) . \end{matrix}

(4)

Using the above GFDM formulae, the discretized version of (2) for each time step

j = 1, . . ., J

is as follows:

\{\begin{matrix} - λ_{00} U_{0}^{j} + \sum_{i = 1}^{s} λ_{i 0} U_{i}^{j} = 0, & (x_{0}, y_{0}) \in Ω, \\ \frac{{(U_{0}^{j})}^{1 / m} - {(U_{0}^{j - 1})}^{1 / m}}{Δ t} = - λ_{02} U_{0}^{j - 1} + \sum_{i = 1}^{s} λ_{i 2} U_{i}^{j - 1}, & (x_{0}, y_{0}) \in Γ_{1}, \\ U_{0}^{j} = 0, & (x_{0}, y_{0}) \in Γ_{2} . \end{matrix}

(5)

Here,

Γ_{1}

denotes the border

{y = 0, t > 0}

and

Γ_{2} = \partial Ω \ Γ_{1}

. For our computational examples, the domains of Figure 1 are considered, where the distribution of the nodes is irregular. More precisely,

Γ_{1} = [- 2, 2]

and

Ω = [- 2, 2] \times Γ_{2}

. The first mesh is generated in such a way that the number of nodes is higher in the centre, so the numerical solution is more accurate. The second mesh has a very irregular distribution of nodes near the border

Γ_{2} .

The second equation in system (5) is explicit in the time variable, i.e., it only depends on the solution of the numerical method in the previous time step. Therefore, to initialize the numerical method, one uses the solution of

\{\begin{matrix} - λ_{00} U_{0}^{0} + \sum_{i = 1}^{s} λ_{i 0} U_{i}^{0} = 0, & (x_{0}, y_{0}) \in Ω, \\ U_{0}^{0} = f^{m} (x_{0}), & (x_{0}, y_{0}) \in Γ_{1}, \\ U_{0}^{0} = 0, & (x_{0}, y_{0}) \in Γ_{2} . \end{matrix}

(6)

For the sake of simplicity, we introduce the following notation:

d_{i} = \sqrt{h_{i}^{2} + k_{i}^{2}}

(with

i \in {1, \dots, s}

) and

d = max_{i} {d_{i}}

for each star.

The paper is organized as follows: in Section 2 the order of the local truncation error and the consistency are proved. Using this together with a discrete maximum principle, the existence and uniqueness of solutions to the problem is shown. Section 3 is devoted to the proof of our convergence result. Finally, in Section 4, numerical examples are shown to validate the theoretical results.

2. Local Truncation Error and Consistency

Now, we substitute the solution u of the continuous problem (2) to be solved in the bounded domain

Ω

into the numerical scheme (5) and denote the local truncation error obtained by

τ_{i}^{j}

. Let us define

A = max_{i, j} | τ_{i}^{j} | .

The next theorem gives the order of the truncation error.

Theorem 1.

Assume that there exist two constants

C_{1}, C_{2} > 0

such that

C_{1} d^{2} < Δ t < C_{2} d^{2}

. The order of the local truncation error is

A = O (Δ t (d^{2} + Δ t))

Proof.

In all the following the fact that the theoretical solution is smooth enough is taken into consideration.

It is easy to observe that, in the boundary nodes on

Γ_{2}

, the local truncation error vanishes due to the imposed hypothesis on the solution. For the interior nodes, it is

τ_{0}^{j} = - λ_{00} u_{0}^{0} + \sum_{i = 1}^{s} λ_{i 0} u_{i}^{0} = Δ u (x_{0}, y_{0}, t_{j}) + O (d^{2}) = O (Δ t (d^{2} + Δ t))

and for the

Γ_{1}

nodes, it yields

\begin{matrix} τ_{0}^{j} & = Δ t (- λ_{02} u_{0}^{j} + \sum_{i = 1}^{s} λ_{i 2} u_{i}^{j}) + {(u_{0}^{j})}^{1 / m} - {(u_{0}^{j + 1})}^{1 / m} \\ = Δ t (\frac{\partial u}{\partial y} (x_{0}, y_{0}, t_{j}) + O (d^{2})) - Δ t (\frac{\partial u^{1 / m}}{\partial t} (x_{0}, y_{0}, t_{j}) + O (Δ t)) \\ = O (Δ t d^{2}) + O (Δ t^{2}) = O (Δ t (d^{2} + Δ t)) . \end{matrix}

(7)

□

Define the following variable

ρ = max_{x} {f^{m} (x)}

and denote by

ψ (x) = x^{m}

. In order to obtain the existence and uniqueness of the numerical solution

U_{i}^{j}

of (5) we need a discrete maximum principle given by the next theorem.

Theorem 2.

Let

U_{i}^{j}

be a solution of (5) with

m \geq 1

. Define

C (m, f) = {[m {(ρ)}^{m - 1}]}^{- 1}

and assume the following relation

Δ t \leq C (m, f) d^{2}

. Then, for every i and j, one has

0 \leq U_{i}^{j} \leq ρ .

Proof.

A maximum principle can be derived at the boundary nodes because at each time step we have a discrete harmonic extension problem. That is why the interior nodes are automatically smaller than they are. At each time step we proceed by induction. Observe that

0 \leq U_{0}^{0} \leq ρ .

By assuming

0 \leq U_{0}^{j} \leq ρ

, we obtain

{(U_{0}^{j + 1})}^{1 / m} = Δ t (- λ_{02} U_{0}^{j} + \sum_{i = 1}^{s} λ_{i 2} U_{i}^{j}) + {(U_{0}^{j})}^{1 / m},

(8)

and for the new variable

z_{i}^{j} = {(U_{0}^{j})}^{1 / m},

it can be obtained

z_{0}^{j + 1} = Δ t (- λ_{02} {(z_{0}^{j})}^{m} + \sum_{i = 1}^{s} λ_{i 2} {(z_{i}^{j})}^{m}) + z_{0}^{j} = Δ t (- \sum_{i = 1}^{s} λ_{i 2} {(z_{0}^{j})}^{m} + \sum_{i = 1}^{s} λ_{i 2} {(z_{i}^{j})}^{m}) + z_{0}^{j} .

(9)

By the Mean Value Theorem,

{(z_{i}^{j})}^{m} - {(z_{0}^{j})}^{m} = (z_{i}^{j} - z_{0}^{j}) ψ^{'} (ξ)

, for some

ξ \in (z_{i}^{j}, z_{0}^{j})

, rewrite (9) as

z_{0}^{j + 1} = ψ^{'} (ξ) Δ t \sum_{i = 1}^{s} λ_{i 2} z_{i}^{j} + (1 - ψ^{'} (ξ) Δ t \sum_{i = 1}^{s} λ_{i 2}) z_{0}^{j} .

(10)

Thanks to the induction hypothesis and the value of the constant

C (m, f)

, we obtain

ψ^{'} (ξ) Δ t λ_{02} \leq 1

and the following estimate holds

\begin{matrix} | z_{0}^{j + 1} | & \leq ψ^{'} (ξ) Δ t \sum_{i = 1}^{s} λ_{i 2} | z_{i}^{j} | + (1 - ψ^{'} (ξ) Δ t \sum_{i = 1}^{s} λ_{i 2}) | z_{0}^{j} | \\ \leq ψ^{'} (ξ) Δ t \sum_{i = 1}^{s} λ_{i 2} ρ^{1 / m} + (1 - ψ^{'} (ξ) Δ t \sum_{i = 1}^{s} λ_{i 2}) ρ^{1 / m} = ρ^{1 / m} . \end{matrix}

For

z_{0}^{j} \geq 0

the same result is obtained by performing similar steps. □

At this stage, the result of existence and uniqueness of solution can be presented.

Theorem 3.

Under the restriction

Δ t \leq C (m, f) d^{2},

the discrete scheme (5) has a unique solution.

Proof.

Proving the uniqueness of solutions is enough due to the linearity of the system of equations, the existence and the uniqueness are equivalent.

Suppose that there exist two different solutions

U_{0}^{j}

and

V_{0}^{j}

of (5) and we denote with

W_{0}^{j} = U_{0}^{j} - V_{0}^{j}

the difference between them, which, in principle, is non-zero. It is easy to see that

W_{0}^{0}

satisfies (5) with

f \equiv 0

and so, by the discrete maximum principle, we obtain

0 \leq W_{0}^{0} \leq 0

. Proceeding by induction we obtain

W_{0}^{j} = 0

; thus,

U_{0}^{j}

and

V_{0}^{j}

coincide and the uniqueness is obtain. □

3. Convergence of the Numerical Solution

With the above results we hereby provide the main result of the article, the convergence of the numerical solution of (5) to the theoretical solution of the Equation (2). The error of the numerical method is given by

e_{i}^{j} = u (x_{i}, y_{i}, t_{j}) - U_{i}^{j}, e^{j} = max_{i} | e_{i}^{j} | .

Theorem 4.

Let u be the solution to (2) and

U_{0}^{j}

be the solution to system (5) with

m \geq 1

. Assume that there exist two constants

C (m, f)

D > 0

such that

D d^{2} \leq Δ t \leq C (m, f) d^{2} .

Then,

e^{j} = O (d^{2} + Δ t) .

Proof.

The selection for the boundary conditions in

Γ_{2}

in the numerical method produces a zero error as in the local truncation error.

In addition,

e^{j} \leq max {e_{Γ_{1}}^{j}, O (d^{2})} \leq e_{Γ_{1}}^{j} + O (d^{2})

because of the approximation for the Laplacian of second order. To compute

e_{Γ_{1}}^{j}

, one has the Equations (7) (the first identity) and (8). Rewriting them in terms of z and

Z_{0}^{j}

and subtracting them the following is obtained

\begin{matrix} e_{0}^{j + 1} & = Δ t [- λ_{02} ({(z_{0}^{j})}^{m} - {(Z_{0}^{j})}^{m}) + \sum_{i = 1}^{s} λ_{i 2} ({(z_{i}^{j})}^{m} - {(Z_{i}^{j})}^{m})] + e_{0}^{j} - τ_{0}^{n} \\ = Δ t [- λ_{02} ψ^{'} (ξ_{0}) e_{0}^{j} + \sum_{i = 1}^{s} λ_{i 2} ψ^{'} (x i_{i}) e_{i}^{j}] + e_{0}^{j} - τ_{0}^{n}, \end{matrix}

(11)

for some

ξ_{0} \in (z_{0}^{j}, Z_{0}^{j})

and

ξ_{i} \in (z_{i}^{j}, Z_{i}^{j})

. Then,

e^{j + 1} \leq e^{j} (1 - Δ t λ_{02} ψ^{'} (ξ_{0})) + Δ t \sum_{i = 1}^{s} λ_{i 2} ψ^{'} (ξ_{i}) \cdot e^{j} + A .

(12)

There exists a constant (assuming enough regularity of the solution z)

K \geq 0

such that

| z_{0}^{j} - z_{i}^{j} | \leq K d

and

| Z_{0}^{j} - Z_{i}^{j} | \leq K d

, then since

\begin{matrix} ψ^{'} (ξ_{i}) & = \sum_{l = 1}^{m} {[z_{i}^{j}]}^{m - l} {[Z_{i}^{j}]}^{l - 1} \leq \sum_{l = 1}^{m} {[z_{0}^{j} + K d]}^{m - l} {[Z_{0}^{j} + K d]}^{l - 1} \\ \sum_{l = 1}^{m} {[z_{0}^{j}]}^{m - l} {[Z_{0}^{j}]}^{l - 1} + L d = ψ^{'} (ξ_{0}) + L d, \end{matrix}

(13)

we have

e^{j + 1} \leq e^{j} (1 - Δ t \sum_{i = 1}^{s} λ_{i 2} ψ^{'} (ξ_{0}) + \sum_{i = 1}^{s} λ_{i 2} ψ^{'} (ξ_{0}) + \sum_{i = 1}^{s} λ_{i 2} L d) + A = e^{j} [1 + λ_{02} L d] + A .

(14)

Remember that

(N + 1) Δ t = T

, where T was the final time and

N + 1

the number of elements in the time discretization. Then, the last expression can be rewritten, since

A = O (Δ t^{2})

, as

\begin{matrix} e^{N + 1} & \leq e^{N} [1 + C \frac{1}{N + 1}] + R Δ t^{2} \leq [1 + C \frac{1}{N + 1}] [e^{N} + R Δ t^{2}] \\ \leq [1 + C \frac{1}{N + 1}] [e^{N - 1} + R Δ t^{2} + R Δ t^{2}] \\ \leq {[1 + C \frac{1}{N + 1}]}^{2} [e^{N - 1} + 2 R Δ t^{2}] \\ \leq \dots \leq {[1 + C \frac{1}{N + 1}]}^{N + 1} [e^{0} + R Δ t^{2}] . \end{matrix}

(15)

However,

{(1 + C \frac{1}{N + 1})}^{N + 1} \leq e^{C}

and

e^{0} \leq D Δ t^{2},

for some

D > 0,

e^{N + 1} \leq e^{C} [D Δ t^{2} + T Δ t]

, and finally

e^{N + 1} = O (Δ t) = O (d^{2} + Δ t) .

□

4. Numerical Examples

Two numerical examples using the meshes of Figure 1 are performed in this section. The first one shows that the numerical solution given by the GFDM presents the well-known Barenblatt profiles. In the second, the error is computed using the two formulas:

l_{2} = \sqrt{\frac{\sum_{i = 1}^{P} {(W_{i} - w_{i})}^{2}}{P}}, l_{\infty} = max_{i} {| W_{i} - w_{i} |} .

Sensitivity analysis for the GFDM were performed in [19], where the authors considered different factors such as number of nodes of the stars, time step and weighting functions

δ

4.1. Example 1: Barenblatt Profiles

In this section we perform several examples showing the application of the described method over the domain given in Figure 1. For our computations we use the initial data

w (x, 0) = f (x) = e^{- \frac{1}{(1 - x) (1 + x)}} χ_{[- 1, 1]} (x),

where

χ_{[- 1, 1]} (x) = \{\begin{matrix} 1, & x \in [- 1, 1], \\ 0, & otherwise . \end{matrix}

Based on the cited reference [6] some Barenblatt profiles as in Vázquez [16] are expected.

In Figure 2 we show the numerical solutions for three different values of m (

m = 1, 3

and 5) for the rectangular domain of Figure 1. As expected, the slow diffusion occurs when m increases. Similarly, Figure 3 contains the cases

m = 2, 4

and 10 for the parabolic domain of Figure 1.

4.2. Example 2

For the following example we consider

m = 1

and

f (x) = \frac{1}{π} \frac{1}{{| x |}^{2} + 1}

[- 2, 2]

. The exact solution is

w (x, t) = \frac{1}{π} \frac{t + 1}{{| x |}^{2} + {(t + 1)}^{2}},

(see [6,20]). The errors of the numerical method are shown in Table 1 for times

t = 3, 5

and 10 s.

In order to obtain a general idea of the performance of the scheme, in [6] the author obtained an error of order

4.329 \times 10^{- 4}

in our second example. In [7], the authors solved a similar equation with a forcing term and a time Caputo fractional derivative obtaining

5.28 \times 10^{- 2}

and

1.50 \times 10^{- 2}

errors for small space steps

Δ x = Δ y

. The numerical solution given the GFDM are plotted in Figure 4 for times

t = 3, 5

and 10.

5. Conclusions and Future Work

The aim of this paper is to study the main properties of the discrete solution given by the generalized finite difference method of the porous media problem (1) with fractional laplacian. To do so, the Silvestre–Caffarelli extension is used. The numerical solution given by the GFDM allows us to discretize the new computational domain with an irregular distribution of nodes.

The order of the truncation error is obtained as well as a discrete maximum principle. Under some conditions in the time step, the convergence of the GFDM scheme is proved. Finally, some numerical examples are given, showing the application of the proposed scheme for solving this fractional porous media equation.

Some possible lines of future work are the following:

The extension of the proposed method for higher dimensional settings.
The application of the above procedure for solving the fractional laplacian equation.

Author Contributions

Conceptualization, M.N. and A.M.V.; methodology, F.U.; software, Á.G.; validation, F.U., A.M.V.; formal analysis, F.U.; investigation, M.N.; resources, Á.G.; data curation, M.N.; writing—original draft preparation, A.M.V.; writing—review and editing, M.N.; visualization, Á.G.; supervision, F.U.; funding acquisition, F.U. and M.N. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge the support of the Escuela Técnica Superior de Ingenieros Industriales (UNED) of Spain, project 2021-IFC02. This work is also partially support by the Project MTM2017-83391-P DGICT, Spain.

Conflicts of Interest

The authors declare no conflict of interest.

References

Vázquez, J.L. The Porous Medium Equation: Mathematical Theory; Oxford Mathematical Monographs, The Clarendon Press, Oxford University Press: Oxford, UK, 2007; ISBN 978-0-19-856903-9. [Google Scholar]
Vázquez, J.L. The Mathematical Theories of Diffusion: Nonlinear and Fractional Diffusion. In Nonlocal and Nonlinear Diffusions and Interactions: New Methods and Directions; Bonforte, M., Grillo, G., Eds.; Lecture Notes in Mathematics; Springer: Cham, Switzerland, 2017; Volume 2186. [Google Scholar]
De Pablo, A.; Quirós, F.; Rodríguez, A.; Vázquez, J.L. A fractional porous medium equation. Adv. Math. 2011, 226, 1378–1409. [Google Scholar] [CrossRef] [Green Version]
Valdinoci, E. From the long jump random walk to the fractional laplacian. SeMA J. Boletín Soc. Española Matemática Apl. 2009, 49, 33–44. [Google Scholar]
Landkof, N.S. Foundations of Modern Potential Theory; Translated from the Russian by Doohovskoy, A.P. and Die Grundlehren der mathematischen Wissenschaften Band; Springer: New York, NY, USA, 1972; Volume 180. (In Russian) [Google Scholar]
Del Teso, F. Finite difference method for a fractional porous medium equation. Calcolo 2014, 51, 615–638. [Google Scholar] [CrossRef] [Green Version]
Hu, Y.; Li, C.; Li, H. The finite difference method for Caputo-type parabolic equation with fractional Laplacian: One-dimension case. Chaos Solitons Fractals 2017, 102, 319–326. [Google Scholar] [CrossRef]
Arendt, W.; Elst, A.F.M.T.; Warma, M. Fractional powers of sectorial operators via the Dirichlet-to-Neumann operator. Commun. Partial Differ. Equ. 2018, 43, 1–24. [Google Scholar] [CrossRef]
Caffarelli, L.; Silvestre, L. An extension problem related to the fractional Laplacian. Commun. Partial Differ. Equ. 2007, 32, 1245–1260. [Google Scholar] [CrossRef] [Green Version]
Chen, S.; Shen, J. An Efficient and Accurate Numerical Method for the Spectral Fractional Laplacian Equation. J. Sci. Comput. 2020, 82, 17. [Google Scholar] [CrossRef]
Padgett, J.L. Analysis of an Approximation to a Fractional Extension Problem. BIT Numer. Math. 2020, 60, 715–739. [Google Scholar] [CrossRef] [Green Version]
D’Elia, M.; Gunzburger, M. The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator. Comput. Math. Appl. 2013, 66, 1245–1260. [Google Scholar] [CrossRef] [Green Version]
Huang, Y.; Oberman, A. Numerical methods for the fractional Laplacian: A finite difference-quadrature approach. SIAM J. Numer. Anal. 2014, 52, 3056–3084. [Google Scholar] [CrossRef]
Nochetto, R.H.; Otarola, E.; Salgado, A.J. A PDE approach to space-time fractional parabolic problems. SIAM J. Numer. Anal. 2016, 54, 848–873. [Google Scholar] [CrossRef] [Green Version]
Ureña, F.; Gavete, L.; Garcia, A.; Benito, J.J.; Vargas, A.M. Solving second order non-linear parabolic PDEs using generalized finite difference method (GFDM). J. Comput. Appl. Math. 2019, 354, 221–241. [Google Scholar] [CrossRef]
Vázquez, J.L. Barenblatt solutions and asymptotic behaviour for a nonlinear fractional heat equation of porous medium type. J. Eur. Math. Soc. 2014, 16, 769–803. [Google Scholar] [CrossRef] [Green Version]
Benito, J.J.; Garcia, A.; Gavete, L.; Negreanu, M.; Ureña, F.; Vargas, A.M. On the numerical solution to a parabolic-elliptic system with chemotactic and periodic terms using Generalized Finite Differences. Eng. Anal. Bound. Elem. 2020, 113, 181–190. [Google Scholar] [CrossRef]
Vargas, A.M. Finite difference method for solving fractional differential equations at irregular meshes. Math. Comput. Simul. 2021. [Google Scholar] [CrossRef]
Benito, J.J.; Ureña, F.; Gavete, L. Influence of several factors in the generalized finite difference method. Appl. Math. Model. 2001, 25, 1039–1053. [Google Scholar] [CrossRef]
Luchko, Y. Wave-diffusion dualism of the neutral-fractional processes. J. Comput. Phys. 2015, 293, 40–52. [Google Scholar] [CrossRef]

Figure 1. Clouds of nodes with irregular distribution.

Figure 2. Solutions given by the method for

m = 1

(left)

m = 3

(centre) and

m = 5

(right) using the first cloud of Figure 1.

Figure 2. Solutions given by the method for

m = 1

(left)

m = 3

(centre) and

m = 5

(right) using the first cloud of Figure 1.

Figure 3. Solutions given by the method for

m = 2

(left)

m = 4

(centre) and

m = 10

(right) using the second cloud of Figure 1.

Figure 3. Solutions given by the method for

m = 2

(left)

m = 4

(centre) and

m = 10

(right) using the second cloud of Figure 1.

Figure 4. Solutions given by the method for

t = 3, 5

and 10 in the second example.

Figure 4. Solutions given by the method for

t = 3, 5

and 10 in the second example.

Table 1.

l^{2}

and

l^{\infty}

norms of the errors in Example 2.

Table 1.

l^{2}

and

l^{\infty}

norms of the errors in Example 2.

$t (s)$	3	5	10
$l_{2}$	$3.7489 \times 10^{- 2}$	$3.9929 \times 10^{- 2}$	$3.7068 \times 10^{- 2}$
$l_{\infty}$	$6.2268 \times 10^{- 2}$	$6.2268 \times 10^{- 2}$	$6.2268 \times 10^{- 2}$

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

García, Á.; Negreanu, M.; Ureña, F.; Vargas, A.M. A Note on a Meshless Method for Fractional Laplacian at Arbitrary Irregular Meshes. Mathematics 2021, 9, 2843. https://doi.org/10.3390/math9222843

AMA Style

García Á, Negreanu M, Ureña F, Vargas AM. A Note on a Meshless Method for Fractional Laplacian at Arbitrary Irregular Meshes. Mathematics. 2021; 9(22):2843. https://doi.org/10.3390/math9222843

Chicago/Turabian Style

García, Ángel, Mihaela Negreanu, Francisco Ureña, and Antonio M. Vargas. 2021. "A Note on a Meshless Method for Fractional Laplacian at Arbitrary Irregular Meshes" Mathematics 9, no. 22: 2843. https://doi.org/10.3390/math9222843

APA Style

García, Á., Negreanu, M., Ureña, F., & Vargas, A. M. (2021). A Note on a Meshless Method for Fractional Laplacian at Arbitrary Irregular Meshes. Mathematics, 9(22), 2843. https://doi.org/10.3390/math9222843

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

A Note on a Meshless Method for Fractional Laplacian at Arbitrary Irregular Meshes

Abstract

1. Introduction

2. Local Truncation Error and Consistency

3. Convergence of the Numerical Solution

4. Numerical Examples

4.1. Example 1: Barenblatt Profiles

4.2. Example 2

5. Conclusions and Future Work

Author Contributions

Funding

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI