# Numerical Solution to Coupled Burgers’ Equations by Gaussian-Based Hermite Collocation Scheme

Authors: Nissaya Chuathong; Sayan Kaennakham
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

One of the most challenging PDE forms in fluid dynamics namely Burgers equations is solved numerically in this work. Its transient, nonlinear, and coupling structure are carefully treated. The Hermite type of collocation mesh-free method is applied to the spatial terms and the 4th-order Runge Kutta is adopted to discretize the governing equations in time. The method is applied in conjunction with the Gaussian radial basis function. The effect of viscous force at high Reynolds number up to 1,300 is investigated using the method. For the purpose of validation, a conventional global collocation scheme (also known as “Kansa” method) is applied parallelly. Solutions obtained are validated against the exact solution and also with some other numerical works available in literature when possible.


---

## Body

## 1. Introduction

Fluid is known to dominate most part of the planet and can be modelled mathematically by the well-known Navier-Stoke equation. To understand its behavior under a certain set of conditions, one seeks solutions (if any or if possible) of the equation and this is not an easy task. The combination of challenging factors, nonlinearity, coupling, and transient nature, of the equation makes it one of the most complex systems both analytically and numerically.Acknowledged as a simple case of the Navier-Stoke equation, the famous “Burgers” equations, a system of equations that describes the interaction between two crucial physical manners of nature, convection, and diffusion, are used to model variety of applications. This includes flows through a shock wave traveling in viscous fluid, the phenomena of turbulence, sedimentation of two kinds of particles in fluid suspensions under the effect of gravity, and also dispersion of pollutants in rivers (see [1–3]).The standard form of this type of equation is expressed as(1)
∂
u
∂
t
+
u
∂
u
∂
x
+
v
∂
u
∂
y
=
1
R
e
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
(2)
and
∂
v
∂
t
+
u
∂
v
∂
x
+
v
∂
v
∂
y
=
1
R
e
∂
2
v
∂
x
2
+
∂
2
v
∂
y
2The ratio of inertial forces to viscous forces is represented by the Reynolds number, Re, and the x- and y-velocity components are noted by u
x
,
y
,
t and v
x
,
y
,
t, respectively.For a domainΩ
=
x
,
y
:
a
≤
x
,
y
≤
b with its boundary Γ, the above system is subject to the initial conditions(3)
u
x
,
y
,
0
=
β
1
x
,
y
,
x
,
y
∈
Ω
(4)
v
x
,
y
,
0
=
β
2
x
,
y
,
x
,
y
∈
Ωand the boundary conditions(5)
u
x
,
y
,
t
=
γ
1
x
,
y
,
t
,
x
,
y
∈
Γ
,
t
>
0
,
(6)
v
x
,
y
,
t
=
γ
2
x
,
y
,
t
,
x
,
y
∈
Γ
,
t
>
0
,Here β
1
,
β
2
,
γ
1, and γ
2 are known functions.In 1983, Fletcher [4] applied the Hopf-Cole transformation and successfully solved the equation analytically. Ever since, the coupled two-dimensional Burgers’ Equations have been attracting a number of investigations both numerically and analytically. Some of the recent computational researches include the application of a fully implicit finite-difference (see [5]), by Adomain-Pade technique (see [6]), by adopting variational iteration method (see [7]), and by applying a mesh-free technique (see [8, 9]). Even more recently, more numerical works have been successfully carried out and nicely documented in literature: the application of the dual reciprocity boundary element method (see [10, 11]), POD/discrete empirical interpolation method- (DEIM-) reduced-order model (ROM) (see [12]), and a combination between the homotopy analysis method and finite differences (see [13]).Amongst the well-known numerical scheme, finite volume, finite difference, and finite element method that have been invented, developed, and applied in a wide range of science and engineering problems, a rather young idea was discovered and has been categorized as “meshless/mesh-free” methods. The methods under this category have recently become promising alternative tools for numerically solving variety of science and engineering problems. Generally, these meshless schemes can be grouped into two main classes:(i)
Strong forms: the finite point method (see [14]), The hp-meshless cloud method (see [15]), the collocation method (see [16]), and references therein.
(ii)
Weak forms: the diffuse element method (see [17]), the element-free Galerkin method (EFGM; see [18]), the point interpolation method (see [19]), and references therein.This depends on the requirement for integration. In this work, particularly, the main attention is paid to the development of the so-called “Collocation Method” characterized in the strong form family. The method is truly mesh-free meaning that no mesh is at all required at any point of the whole computing process. Amongst several versions of its further improvements, this work focuses on an application of the Hermite form of the conventional collocation where the involved differential operator is applied twice, explained in more detail in Section2.The paper is organized as follows. Section2 provides the necessary mathematical background concerning the Gaussian radial basis function and also the fundamental idea of collocation mesh-free schemes in general. The implementation of the schemes to the couple-Burgers equations is then detailed in Section 3. Section 4 demonstrates the numerical experiments and discusses general aspect of the application before some interesting findings are listed in Section 5.


## 2. The Gaussian-Based Collocation Meshless Method

Lettingx
c be a point in a d-dimensional Euclidean space R
d, for any point x
∈
R
d that depends only on the distance, r, from the fixed x
c, then the function φ
:
R
d
→
R that is written as (7)
φ
r
,
ε
=
φ
x
-
x
cis called a Radial Basis Function (RBF), with x
c being called the center and some nonnegative parameter ε.Among different types of RBFs invented, improved, and applied nowadays, in this work, one of the most widely used forms of radial functions, known as “Gaussian,” is focused on. The standard form of this type of RBF is expressed as(8)
φ
r
=
exp
⁡
-
ε
r
2where r
∈
R and ε is called “shape parameter,” that is related to the variance of the normal distribution function by ε
2
=
1
/
2
δ
2. Under the context of collocation schemes, this relation is linked to the Euclidean distance function, .
2, via the following multivariate function:(9)
Ξ
k
x
=
e
-
χ
with
χ
=
ε
2
x
-
x
k
2
2for a fixed center x
∈
R
d and x
k
∈
R
d. It can then be seen that the connection between and is as follows:(10)
Ξ
k
x
=
φ
x
-
x
k
2The collocation scheme starts with considering the following elliptical partial differential equation defined on a bounded and connected domainΩ:(11)
Φ
u
x
=
f
x
for
x
∈
Ω
⊂
R
n
(12)
B
1
u
x
=
g
x
for
x
∈
Γ
1
(13)
B
1
u
x
=
h
x
for
x
∈
Γ
1where ∂
Ω is the domain boundary containing two nonoverlap sections: Γ
1 and Γ
2, with Γ
1
∩
Γ
2
=
ϕ. These differential operators, Φ and B
1
,
B
2 are applied on the domain and the two boundary sections, respectively. Three known functions f
(
x
)
,
g
(
x
)
,
a
n
d
h
(
x
) can well be dependent of space and/or time. Let X
c
=
x
j
j
=
1
N be a set of randomly selected points, known as “collocation” or “centers”, on the domain, where x
j
j
=
1
N
i are those contained within, where x
j
j
=
N
i
+
1
N
i
+
N
1 and x
j
j
=
N
i
+
N
1
+
1
N are those on the boundaries Γ
1 and Γ
2, respectively.


### 2.1. The Conventional Kansa Collocation Method

Initially, a collocation scheme writes the approximate solution,u
~
(
x
), as a linear combination of the basis function φ
:
-
j
N, shown in the following form:(14)
u
x
≃
u
~
x
=
∑
j
=
1
N
α
j
φ
x
-
x
j
2where α
j are coefficients and .
2 is the Euclidean norm. The basis function used now is the radial type as defined previously. Applying Φ and B
1
,
B
2 to both domain and boundary sections, satisfying the governing system of equations, allows the system to arrive at(15)
A
α
=
bwhere α
=
α
1
α
2
.
.
.
α
N, and the known b vector expressed is as follows:(16)
b
=
f
x
1
f
x
2
⋯
f
x
N
i
g
x
N
i
+
1
⋯
g
x
N
i
+
N
1
h
x
N
i
+
N
1
+
1
⋯
h
x
N
Tand by setting φ to be a matrix with entries φ
i
j
=
φ
x
i
-
x
j
2 for i
,
j
=
1,2
,
.
.
.
,
N we have (17)
A
=
Φ
φ
B
1
φ
B
2
φOnce α
j are obtained by (15), the approximate solution is straightforward yielded. This method is known as “Kansa,” in honor of a great mathematician E. J. Kansa [22] who discovered the idea in 1990. The method has been applied to a wide range of problems ever since [23, 24]. It, nevertheless, has shortcomings as it is known to suffer the problem of unsymmetric interpolation matrix, A, and the rigorous mathematical proof of its solvability is still not available [25], and very often produces low quality results particularly in boundary-adjacent region [26, 27].


### 2.2. The Hermite Collocation Method

In order to alleviate the difficulty mentioned earlier, in 1997, Fasshauer [28] proposed a new way of interpolating by applying the self-adjoint operators Φ
∗ and B
1
∗
,
B
2
∗ to the governing system of equations and rewriting the approximate solution (14) as(18)
u
x
≃
u
~
x
=
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
i
+
N
1
α
j
B
1
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
N
1
+
1
N
α
j
B
2
∗
φ
x
-
x
j
2This leads to a new interpolation matrix A, shown as follows:(19)
A
=
Φ
Φ
∗
φ
Φ
B
1
∗
φ
Φ
B
2
∗
φ
B
1
Φ
∗
φ
B
1
B
1
∗
φ
B
1
B
2
∗
φ
B
2
Φ
∗
φ
B
2
B
1
∗
φ
B
2
B
2
∗
φIn this study, nevertheless, the problem domain has only one continuous boundary with no differential operator. An application of this Hermite type of collocation method to elastostatic problem was successfully investigated by Leitao [29]. Additionally, some interesting implementations of the scheme to the transient and nonlinear plate problems can be found in [30–32].


## 2.1. The Conventional Kansa Collocation Method

Initially, a collocation scheme writes the approximate solution,u
~
(
x
), as a linear combination of the basis function φ
:
-
j
N, shown in the following form:(14)
u
x
≃
u
~
x
=
∑
j
=
1
N
α
j
φ
x
-
x
j
2where α
j are coefficients and .
2 is the Euclidean norm. The basis function used now is the radial type as defined previously. Applying Φ and B
1
,
B
2 to both domain and boundary sections, satisfying the governing system of equations, allows the system to arrive at(15)
A
α
=
bwhere α
=
α
1
α
2
.
.
.
α
N, and the known b vector expressed is as follows:(16)
b
=
f
x
1
f
x
2
⋯
f
x
N
i
g
x
N
i
+
1
⋯
g
x
N
i
+
N
1
h
x
N
i
+
N
1
+
1
⋯
h
x
N
Tand by setting φ to be a matrix with entries φ
i
j
=
φ
x
i
-
x
j
2 for i
,
j
=
1,2
,
.
.
.
,
N we have (17)
A
=
Φ
φ
B
1
φ
B
2
φOnce α
j are obtained by (15), the approximate solution is straightforward yielded. This method is known as “Kansa,” in honor of a great mathematician E. J. Kansa [22] who discovered the idea in 1990. The method has been applied to a wide range of problems ever since [23, 24]. It, nevertheless, has shortcomings as it is known to suffer the problem of unsymmetric interpolation matrix, A, and the rigorous mathematical proof of its solvability is still not available [25], and very often produces low quality results particularly in boundary-adjacent region [26, 27].


## 2.2. The Hermite Collocation Method

In order to alleviate the difficulty mentioned earlier, in 1997, Fasshauer [28] proposed a new way of interpolating by applying the self-adjoint operators Φ
∗ and B
1
∗
,
B
2
∗ to the governing system of equations and rewriting the approximate solution (14) as(18)
u
x
≃
u
~
x
=
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
i
+
N
1
α
j
B
1
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
N
1
+
1
N
α
j
B
2
∗
φ
x
-
x
j
2This leads to a new interpolation matrix A, shown as follows:(19)
A
=
Φ
Φ
∗
φ
Φ
B
1
∗
φ
Φ
B
2
∗
φ
B
1
Φ
∗
φ
B
1
B
1
∗
φ
B
1
B
2
∗
φ
B
2
Φ
∗
φ
B
2
B
1
∗
φ
B
2
B
2
∗
φIn this study, nevertheless, the problem domain has only one continuous boundary with no differential operator. An application of this Hermite type of collocation method to elastostatic problem was successfully investigated by Leitao [29]. Additionally, some interesting implementations of the scheme to the transient and nonlinear plate problems can be found in [30–32].


## 3. Implementation to Burgers’ Equations

As previously mentioned, the coupled Burgers’ equations can be seen as a simple model that mimics the coupling of fluid flow to thermal dynamics, as well as some other scientific fields. High and low Reynolds numbers play important roles in both theoretically modeling and numerical simulation. In this section, the Hermite collocation is now applied based on the Gaussian type of radial basis function.Since the principle implementation of the conventional collocation method to nonlinear and time dependent problem has already been documented in some of our previous studies [33, 34], this section is therefore dedicated to the application of the Hermite version of collocation as explained in details in the following subsections.


### 3.1. Space Discretization with the Hermite Scheme

With Hermite interpolation technique, it begins with writing the approximation of solution foru
~
(
x
) and v
~
(
x
), respectively, with the same radial basis function φ
x
-
x
j
2, as follows: (20)
u
x
,
t
≃
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
t
B
∗
φ
x
-
x
j
2and (21)
v
x
,
t
≃
v
~
x
,
t
=
∑
j
=
1
N
i
β
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
β
j
t
B
∗
φ
x
-
x
j
2Note that it is now assumed that there is only one section of boundary Γ
1
=
Γ
2
=
Γ. Therefore, at n
t
h
- derivative, this leads to(22)
∂
n
u
x
,
t
∂
x
n
≃
∑
j
=
1
N
i
α
j
t
∂
n
Φ
∗
φ
x
-
x
j
2
∂
x
n
+
∑
j
=
N
i
+
1
N
α
j
t
∂
n
B
∗
φ
x
-
x
j
2
∂
x
nand (23)
∂
n
v
x
,
t
∂
y
n
≃
∑
j
=
1
N
i
β
j
t
∂
n
Φ
∗
φ
x
-
x
j
2
∂
y
n
+
∑
j
=
N
i
+
1
N
β
j
t
∂
n
B
∗
φ
x
-
x
j
2
∂
y
nTherefore, the governing equations are now rewritten, at a i
t
h
-center node, as(24)
∂
u
∂
t
i
=
1
R
e
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
x
x
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
x
x
+
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
y
y
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
y
y
-
u
i
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
x
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
x
-
v
i
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
y
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
y
(25)
∂
v
∂
t
i
=
1
R
e
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
x
x
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
x
x
+
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
y
y
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
y
y
-
u
i
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
x
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
x
-
v
i
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
y
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
ywhere (26)
H
i
j
x
x
=
∂
2
φ
x
i
-
x
j
∂
x
2
,
H
i
j
y
y
=
∂
2
φ
x
i
-
x
j
∂
y
2
,
H
i
j
x
=
∂
φ
x
i
-
x
j
∂
x
,
and
H
i
j
y
=
∂
φ
x
i
-
x
j
∂
yleading their matrix forms expresses as below:(27)
u
1
t
u
2
t
⋮
u
N
t
=
ψ
11
ψ
12
⋯
ψ
1
N
ψ
21
ψ
22
⋯
ψ
1
N
⋮
⋮
⋱
⋮
ψ
N
1
ψ
N
2
⋯
ψ
N
N
-
u
1
0
⋯
0
0
u
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
u
N
ζ
11
ζ
12
⋯
ζ
1
N
ζ
21
22
⋯
ζ
2
N
⋮
⋮
⋱
⋮
ζ
N
1
ζ
N
2
⋯
ζ
N
N
-
v
1
0
⋯
0
0
v
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
v
N
ϑ
11
ϑ
12
⋯
ϑ
1
N
ϑ
21
ϑ
22
⋯
ϑ
2
N
⋮
⋮
⋱
⋮
ϑ
N
1
ϑ
N
2
⋯
ϑ
N
N
y
α
1
α
2
⋮
α
NSimilarly,(28)
v
1
t
v
2
t
⋮
v
N
t
=
ψ
11
ψ
12
⋯
ψ
1
N
ψ
21
ψ
22
⋯
ψ
1
N
⋮
⋮
⋱
⋮
ψ
N
1
ψ
N
2
⋯
ψ
N
N
-
u
1
0
⋯
0
0
u
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
u
N
ζ
11
ζ
12
⋯
ζ
1
N
ζ
21
22
⋯
ζ
2
N
⋮
⋮
⋱
⋮
ζ
N
1
ζ
N
2
⋯
ζ
N
N
-
v
1
0
⋯
0
0
v
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
v
N
ϑ
11
ϑ
12
⋯
ϑ
1
N
ϑ
21
ϑ
22
⋯
ϑ
2
N
⋮
⋮
⋱
⋮
ϑ
N
1
ϑ
N
2
⋯
ϑ
N
N
y
β
1
β
2
⋮
β
Nwhere (29)
ψ
i
j
=
1
R
e
Φ
∗
H
i
j
x
x
+
Φ
∗
H
i
j
y
y
;
1
≤
j
≤
N
i
1
R
e
B
∗
H
i
j
x
x
+
B
∗
H
i
j
y
y
;
N
i
+
1
≤
j
≤
N
,
and
ζ
i
j
=
Φ
∗
H
i
j
x
+
B
∗
H
i
j
x
;
1
≤
j
≤
N
i
Φ
∗
H
i
j
x
+
B
∗
H
i
j
x
;
N
i
+
1
≤
j
≤
N
,
and
ϑ
i
j
=
Φ
∗
H
i
j
y
+
B
∗
H
i
j
y
;
1
≤
j
≤
N
i
Φ
∗
H
i
j
y
+
B
∗
H
i
j
y
;
N
i
+
1
≤
j
≤
NConsequently, the systematics of matrices involved in this coupling process is(30)
U
t
=
ψ
-
d
i
a
g
U
ζ
-
d
i
a
g
V
ϑ
α
=
Α
~
α
(31)
V
t
=
ψ
-
d
i
a
g
U
ζ
-
d
i
a
g
V
ϑ
β
=
Β
~
β


### 3.2. Time Discretization

In order to carry out the time-dependent terms, the4
t
h-order Runge Kutta is applied to both terms appearing in the above equations, that is, by setting(32)
u
~
=
u
~
x
1
,
u
~
x
2
,
…
,
u
~
x
N
T
and
v
~
=
v
~
x
1
,
v
~
x
2
,
…
,
v
~
x
N
TThat is,(33)
d
u
~
d
t
=
U
t
=
Α
~
α
=
F
u
~
(34)
d
v
~
d
t
=
V
t
=
Β
~
β
=
G
v
~The time progressing of variables involved can be estimated as follows:(35)
u
t
+
1
=
u
t
+
1
6
k
1
+
2
k
2
+
2
k
3
+
k
4where(36)
k
1
=
Δ
t
F
u
i
,
t
i
,
k
2
=
Δ
t
F
u
i
+
k
1
2
,
t
i
+
Δ
t
2
,
k
3
=
Δ
t
F
u
i
+
k
2
2
,
t
i
+
Δ
t
2
,
k
4
=
Δ
t
F
u
i
+
k
3
,
t
i
+
Δ
twhere the similar manner is applied for the other case, d
v
~
/
d
t, and the systems are treated and advanced in time with an ODE method. Then the approximate solution obtained by this final process is expressed in terms of the radial basis function at N collocation nodes as follows:(37)
u
~
x
=
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
i
+
N
α
j
B
∗
φ
x
-
x
j
2
(38)
and
v
~
x
=
∑
j
=
1
N
i
β
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
i
+
N
β
j
B
∗
φ
x
-
x
j
2


### 3.3. Gaussian-Based Hermite Collocation

Lettingφ
r be a radial basis function depending on a distance function, r
i
j
=
x
i
-
x
j
2, where x
i
,
x
j
∈
R
2, it can be seen that(39)
∂
φ
∂
x
=
x
i
-
x
j
r
d
φ
d
r
,
∂
φ
∂
y
=
y
i
-
y
j
r
d
φ
d
r
,
∂
2
φ
∂
x
2
=
x
i
-
x
j
2
r
2
d
2
φ
d
r
2
+
y
i
-
y
j
3
r
3
d
φ
d
r
,
and
∂
2
φ
∂
y
2
=
y
i
-
y
j
2
r
2
d
2
φ
d
r
2
+
x
i
-
x
j
3
r
3
d
φ
d
r
,This leads to the Laplacian expressed as follows:(40)
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
r
=
d
2
d
r
2
φ
r
+
1
r
d
d
r
φ
rFurthermore, when the Hermite concept of collocation is employed, it is necessary that the Lapacian be applied twice resulting in the following so-called fourth-order biharmonic form; (41)
∂
4
∂
x
4
+
2
∂
4
∂
x
2
∂
y
2
+
∂
4
∂
y
4
φ
r
=
d
4
d
r
4
φ
r
+
2
r
d
3
d
r
3
φ
r
-
1
r
2
d
2
d
r
2
φ
r
+
1
r
3
d
d
r
φ
rFor this work, Gaussian type of radial basis function, is used and its first four orders of derivatives can, be respectively, expressed as(42)
φ
x
-
x
j
2
=
exp
⁡
-
ε
r
2
=
exp
⁡
-
ε
x
-
x
j
2
2
(43)
Thus,
d
d
r
φ
r
=
-
2
r
ε
2
exp
⁡
-
ε
r
2
,
(44)
d
2
d
r
2
φ
r
=
2
ε
2
exp
⁡
-
ε
r
2
2
ε
r
2
-
1
,
(45)
d
3
d
r
3
φ
r
=
2
ε
2
exp
⁡
-
ε
2
r
2
6
ε
2
r
-
4
ε
4
r
3
,
(46)
and
d
4
d
r
4
φ
r
=
2
ε
2
exp
⁡
-
ε
2
r
2
8
ε
6
r
4
-
24
ε
4
r
2
+
6
ε
2for each pair of center nodes: x
i
,
x
j.What follows is the brief demonstration of how the self-adjoint operator works.Hermite starts with the following solution approximation:(47)
u
x
,
t
≃
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
t
B
∗
φ
x
-
x
j
2From (24),(25), (27), and (28), when considering for 1
≤
j
≤
N
i, we have ψ
i
j
=
1
/
R
e
Φ
∗
H
i
j
x
x
+
Φ
∗
H
i
j
y
y. That is, (48)
∂
2
∂
x
2
+
∂
2
∂
y
2
u
x
,
t
≃
∂
2
∂
x
2
+
∂
2
∂
y
2
u
~
x
,
t
=
∂
2
∂
x
2
+
∂
2
∂
y
2
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
Φ
∗
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2For this particular example we have Φ
∗
=
∂
2
/
∂
x
2
+
∂
2
/
∂
y
2 and there is no derivative operation taking place on the boundary. Therefore, the above equation becomes(49)
∂
2
∂
x
2
+
∂
2
∂
y
2
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
∂
4
∂
x
4
+
2
∂
4
∂
x
2
∂
y
2
+
∂
4
∂
y
4
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
d
2
d
r
2
φ
r
+
1
r
d
d
r
φ
r
=
∑
j
=
1
N
i
α
j
d
4
d
r
4
φ
x
-
x
j
2
+
2
r
d
3
d
r
3
φ
x
-
x
j
2
-
1
r
2
d
2
d
r
2
φ
x
-
x
j
2
+
1
r
3
d
d
r
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
d
2
d
r
2
φ
x
-
x
j
2
+
1
r
d
d
r
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
2
ε
2
exp
⁡
-
ε
2
x
-
x
j
2
2
8
ε
6
x
-
x
j
2
4
-
24
ε
4
x
-
x
j
2
2
+
6
ε
2
+
2
2
ε
2
exp
⁡
-
ε
2
x
-
x
j
2
2
6
ε
2
-
4
ε
4
x
-
x
j
2
2
-
1
x
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
1
x
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2For 1
≤
j
≤
N
i, we have ζ
i
j=Φ
∗
H
i
j
x
+
B
∗
H
i
j
x
(50)
u
x
,
t
≃
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
t
B
∗
φ
x
-
x
j
2
∂
∂
x
u
x
,
t
≃
∂
∂
x
u
~
x
,
t
=
∂
∂
x
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
∂
∂
x
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
∂
∂
x
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
Φ
∗
∂
∂
x
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
∂
∂
x
φ
x
-
x
j
2Now, we have Φ
∗
=
∂
/
∂
x and with the absence of the derivative operation on the boundary, B
∗, the remaining term of the above equation can be written as(51)
∂
∂
x
u
x
,
t
≃
∂
∂
x
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
∂
∂
x
∂
φ
x
-
x
j
2
∂
x
+
∑
j
=
N
i
+
1
N
α
j
∂
∂
x
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
x
i
-
x
j
2
x
-
x
j
2
2
d
2
φ
d
r
2
+
y
i
-
y
j
3
x
-
x
j
2
3
d
φ
d
r
+
∑
j
=
N
i
+
1
N
α
j
x
i
-
x
j
r
d
φ
d
r
=
∑
j
=
1
N
i
α
j
x
i
-
x
j
2
x
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
y
i
-
y
j
3
x
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
x
i
-
x
j
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2Similarly, (52)
∂
∂
y
u
x
,
t
≃
∂
∂
y
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
y
i
-
y
j
2
x
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
x
i
-
x
j
3
x
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
y
i
-
y
j
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2By substituting (49), (51), and (52) back into the governing equations(53)
∂
u
∂
t
=
1
R
e
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
-
u
∂
u
∂
x
-
v
∂
u
∂
yTherefore, for each center x
k, the following is reached:(54)
∂
u
∂
t
k
=
1
R
e
∑
j
=
1
N
i
α
j
2
ε
2
exp
⁡
-
ε
2
x
k
-
x
j
2
2
8
ε
6
x
k
-
x
j
2
4
-
24
ε
4
x
k
-
x
j
2
2
+
6
ε
2
+
2
2
ε
2
exp
⁡
-
ε
2
x
k
-
x
j
2
2
6
ε
2
-
4
ε
4
x
k
-
x
j
2
2
-
1
x
k
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
2
ε
x
k
-
x
j
2
2
-
1
+
1
x
k
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
-
u
k
∑
j
=
1
N
i
α
j
x
i
-
x
j
2
x
k
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
2
ε
x
k
-
x
j
2
2
-
1
+
y
i
-
y
j
3
x
k
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
x
i
-
x
j
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
-
v
k
∑
j
=
1
N
i
α
j
y
i
-
y
j
2
x
k
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
2
ε
x
k
-
x
j
2
2
-
1
+
x
i
-
x
j
3
x
k
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
y
i
-
y
j
-
2
ε
2
e
x
p
-
ε
x
-
x
j
2
2The whole process explained above is applied to the other equation of the coupled governing equations. Note that all the terms with1
/
x
k
-
x
j
2 shall be cancelling each other out so there is no zero-division problem in computing process.


### 3.4. Computing Setup and Algorithm

Each computation under this work was completed following the computing steps detailed below.Step 1.
ChooseN collocation or center nodes X
c
=
x
j
j
=
1
N, on the domain Ω.Step 2.
Specify the desired values of(2.1)
The Reynolds number, Re.
(2.2)
The final timet.
(2.3)
The time stepΔ
t.
(2.4)
The Gaussian shape parameterε.Step 3.
Compute the initial collocation matricesΑ
~ and Β
~ using the initial Hermite approximation equations (20)-(21).Step 4.
Apply the initial conditions,U
0
,
V
0, to get α and β from (30), (31), (33), (34), (35), (36), (37), (38), (40), (41), (42), and (43).Step 5.
Compute the solutions for the next time step via the time-discretization explained above.Step 6.
Construct the collocation matricesΑ
~ and Β
~ in the forms expressed in (30), (31), (33), (34), (35), (36), (37), (38), (40), (41), (42), and (43) using the solutions previously obtained in Step 5.Step 7.
Carry on Steps5 and 6 until reaching the final time step t.All computing experiments carried out in this study were executed on a computer notebook: Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz with 8.00 GB of RAM and 64-bit operating system.


## 3.1. Space Discretization with the Hermite Scheme

With Hermite interpolation technique, it begins with writing the approximation of solution foru
~
(
x
) and v
~
(
x
), respectively, with the same radial basis function φ
x
-
x
j
2, as follows: (20)
u
x
,
t
≃
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
t
B
∗
φ
x
-
x
j
2and (21)
v
x
,
t
≃
v
~
x
,
t
=
∑
j
=
1
N
i
β
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
β
j
t
B
∗
φ
x
-
x
j
2Note that it is now assumed that there is only one section of boundary Γ
1
=
Γ
2
=
Γ. Therefore, at n
t
h
- derivative, this leads to(22)
∂
n
u
x
,
t
∂
x
n
≃
∑
j
=
1
N
i
α
j
t
∂
n
Φ
∗
φ
x
-
x
j
2
∂
x
n
+
∑
j
=
N
i
+
1
N
α
j
t
∂
n
B
∗
φ
x
-
x
j
2
∂
x
nand (23)
∂
n
v
x
,
t
∂
y
n
≃
∑
j
=
1
N
i
β
j
t
∂
n
Φ
∗
φ
x
-
x
j
2
∂
y
n
+
∑
j
=
N
i
+
1
N
β
j
t
∂
n
B
∗
φ
x
-
x
j
2
∂
y
nTherefore, the governing equations are now rewritten, at a i
t
h
-center node, as(24)
∂
u
∂
t
i
=
1
R
e
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
x
x
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
x
x
+
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
y
y
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
y
y
-
u
i
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
x
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
x
-
v
i
∑
j
=
1
N
i
α
j
Φ
∗
H
i
j
y
+
∑
j
=
N
i
+
1
N
α
j
B
∗
H
i
j
y
(25)
∂
v
∂
t
i
=
1
R
e
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
x
x
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
x
x
+
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
y
y
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
y
y
-
u
i
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
x
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
x
-
v
i
∑
j
=
1
N
i
β
j
Φ
∗
H
i
j
y
+
∑
j
=
N
i
+
1
N
β
j
B
∗
H
i
j
ywhere (26)
H
i
j
x
x
=
∂
2
φ
x
i
-
x
j
∂
x
2
,
H
i
j
y
y
=
∂
2
φ
x
i
-
x
j
∂
y
2
,
H
i
j
x
=
∂
φ
x
i
-
x
j
∂
x
,
and
H
i
j
y
=
∂
φ
x
i
-
x
j
∂
yleading their matrix forms expresses as below:(27)
u
1
t
u
2
t
⋮
u
N
t
=
ψ
11
ψ
12
⋯
ψ
1
N
ψ
21
ψ
22
⋯
ψ
1
N
⋮
⋮
⋱
⋮
ψ
N
1
ψ
N
2
⋯
ψ
N
N
-
u
1
0
⋯
0
0
u
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
u
N
ζ
11
ζ
12
⋯
ζ
1
N
ζ
21
22
⋯
ζ
2
N
⋮
⋮
⋱
⋮
ζ
N
1
ζ
N
2
⋯
ζ
N
N
-
v
1
0
⋯
0
0
v
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
v
N
ϑ
11
ϑ
12
⋯
ϑ
1
N
ϑ
21
ϑ
22
⋯
ϑ
2
N
⋮
⋮
⋱
⋮
ϑ
N
1
ϑ
N
2
⋯
ϑ
N
N
y
α
1
α
2
⋮
α
NSimilarly,(28)
v
1
t
v
2
t
⋮
v
N
t
=
ψ
11
ψ
12
⋯
ψ
1
N
ψ
21
ψ
22
⋯
ψ
1
N
⋮
⋮
⋱
⋮
ψ
N
1
ψ
N
2
⋯
ψ
N
N
-
u
1
0
⋯
0
0
u
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
u
N
ζ
11
ζ
12
⋯
ζ
1
N
ζ
21
22
⋯
ζ
2
N
⋮
⋮
⋱
⋮
ζ
N
1
ζ
N
2
⋯
ζ
N
N
-
v
1
0
⋯
0
0
v
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
v
N
ϑ
11
ϑ
12
⋯
ϑ
1
N
ϑ
21
ϑ
22
⋯
ϑ
2
N
⋮
⋮
⋱
⋮
ϑ
N
1
ϑ
N
2
⋯
ϑ
N
N
y
β
1
β
2
⋮
β
Nwhere (29)
ψ
i
j
=
1
R
e
Φ
∗
H
i
j
x
x
+
Φ
∗
H
i
j
y
y
;
1
≤
j
≤
N
i
1
R
e
B
∗
H
i
j
x
x
+
B
∗
H
i
j
y
y
;
N
i
+
1
≤
j
≤
N
,
and
ζ
i
j
=
Φ
∗
H
i
j
x
+
B
∗
H
i
j
x
;
1
≤
j
≤
N
i
Φ
∗
H
i
j
x
+
B
∗
H
i
j
x
;
N
i
+
1
≤
j
≤
N
,
and
ϑ
i
j
=
Φ
∗
H
i
j
y
+
B
∗
H
i
j
y
;
1
≤
j
≤
N
i
Φ
∗
H
i
j
y
+
B
∗
H
i
j
y
;
N
i
+
1
≤
j
≤
NConsequently, the systematics of matrices involved in this coupling process is(30)
U
t
=
ψ
-
d
i
a
g
U
ζ
-
d
i
a
g
V
ϑ
α
=
Α
~
α
(31)
V
t
=
ψ
-
d
i
a
g
U
ζ
-
d
i
a
g
V
ϑ
β
=
Β
~
β


## 3.2. Time Discretization

In order to carry out the time-dependent terms, the4
t
h-order Runge Kutta is applied to both terms appearing in the above equations, that is, by setting(32)
u
~
=
u
~
x
1
,
u
~
x
2
,
…
,
u
~
x
N
T
and
v
~
=
v
~
x
1
,
v
~
x
2
,
…
,
v
~
x
N
TThat is,(33)
d
u
~
d
t
=
U
t
=
Α
~
α
=
F
u
~
(34)
d
v
~
d
t
=
V
t
=
Β
~
β
=
G
v
~The time progressing of variables involved can be estimated as follows:(35)
u
t
+
1
=
u
t
+
1
6
k
1
+
2
k
2
+
2
k
3
+
k
4where(36)
k
1
=
Δ
t
F
u
i
,
t
i
,
k
2
=
Δ
t
F
u
i
+
k
1
2
,
t
i
+
Δ
t
2
,
k
3
=
Δ
t
F
u
i
+
k
2
2
,
t
i
+
Δ
t
2
,
k
4
=
Δ
t
F
u
i
+
k
3
,
t
i
+
Δ
twhere the similar manner is applied for the other case, d
v
~
/
d
t, and the systems are treated and advanced in time with an ODE method. Then the approximate solution obtained by this final process is expressed in terms of the radial basis function at N collocation nodes as follows:(37)
u
~
x
=
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
i
+
N
α
j
B
∗
φ
x
-
x
j
2
(38)
and
v
~
x
=
∑
j
=
1
N
i
β
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
i
+
N
β
j
B
∗
φ
x
-
x
j
2


## 3.3. Gaussian-Based Hermite Collocation

Lettingφ
r be a radial basis function depending on a distance function, r
i
j
=
x
i
-
x
j
2, where x
i
,
x
j
∈
R
2, it can be seen that(39)
∂
φ
∂
x
=
x
i
-
x
j
r
d
φ
d
r
,
∂
φ
∂
y
=
y
i
-
y
j
r
d
φ
d
r
,
∂
2
φ
∂
x
2
=
x
i
-
x
j
2
r
2
d
2
φ
d
r
2
+
y
i
-
y
j
3
r
3
d
φ
d
r
,
and
∂
2
φ
∂
y
2
=
y
i
-
y
j
2
r
2
d
2
φ
d
r
2
+
x
i
-
x
j
3
r
3
d
φ
d
r
,This leads to the Laplacian expressed as follows:(40)
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
r
=
d
2
d
r
2
φ
r
+
1
r
d
d
r
φ
rFurthermore, when the Hermite concept of collocation is employed, it is necessary that the Lapacian be applied twice resulting in the following so-called fourth-order biharmonic form; (41)
∂
4
∂
x
4
+
2
∂
4
∂
x
2
∂
y
2
+
∂
4
∂
y
4
φ
r
=
d
4
d
r
4
φ
r
+
2
r
d
3
d
r
3
φ
r
-
1
r
2
d
2
d
r
2
φ
r
+
1
r
3
d
d
r
φ
rFor this work, Gaussian type of radial basis function, is used and its first four orders of derivatives can, be respectively, expressed as(42)
φ
x
-
x
j
2
=
exp
⁡
-
ε
r
2
=
exp
⁡
-
ε
x
-
x
j
2
2
(43)
Thus,
d
d
r
φ
r
=
-
2
r
ε
2
exp
⁡
-
ε
r
2
,
(44)
d
2
d
r
2
φ
r
=
2
ε
2
exp
⁡
-
ε
r
2
2
ε
r
2
-
1
,
(45)
d
3
d
r
3
φ
r
=
2
ε
2
exp
⁡
-
ε
2
r
2
6
ε
2
r
-
4
ε
4
r
3
,
(46)
and
d
4
d
r
4
φ
r
=
2
ε
2
exp
⁡
-
ε
2
r
2
8
ε
6
r
4
-
24
ε
4
r
2
+
6
ε
2for each pair of center nodes: x
i
,
x
j.What follows is the brief demonstration of how the self-adjoint operator works.Hermite starts with the following solution approximation:(47)
u
x
,
t
≃
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
t
B
∗
φ
x
-
x
j
2From (24),(25), (27), and (28), when considering for 1
≤
j
≤
N
i, we have ψ
i
j
=
1
/
R
e
Φ
∗
H
i
j
x
x
+
Φ
∗
H
i
j
y
y. That is, (48)
∂
2
∂
x
2
+
∂
2
∂
y
2
u
x
,
t
≃
∂
2
∂
x
2
+
∂
2
∂
y
2
u
~
x
,
t
=
∂
2
∂
x
2
+
∂
2
∂
y
2
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
Φ
∗
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2For this particular example we have Φ
∗
=
∂
2
/
∂
x
2
+
∂
2
/
∂
y
2 and there is no derivative operation taking place on the boundary. Therefore, the above equation becomes(49)
∂
2
∂
x
2
+
∂
2
∂
y
2
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
∂
2
∂
x
2
+
∂
2
∂
y
2
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
∂
4
∂
x
4
+
2
∂
4
∂
x
2
∂
y
2
+
∂
4
∂
y
4
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
d
2
d
r
2
φ
r
+
1
r
d
d
r
φ
r
=
∑
j
=
1
N
i
α
j
d
4
d
r
4
φ
x
-
x
j
2
+
2
r
d
3
d
r
3
φ
x
-
x
j
2
-
1
r
2
d
2
d
r
2
φ
x
-
x
j
2
+
1
r
3
d
d
r
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
d
2
d
r
2
φ
x
-
x
j
2
+
1
r
d
d
r
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
2
ε
2
exp
⁡
-
ε
2
x
-
x
j
2
2
8
ε
6
x
-
x
j
2
4
-
24
ε
4
x
-
x
j
2
2
+
6
ε
2
+
2
2
ε
2
exp
⁡
-
ε
2
x
-
x
j
2
2
6
ε
2
-
4
ε
4
x
-
x
j
2
2
-
1
x
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
1
x
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2For 1
≤
j
≤
N
i, we have ζ
i
j=Φ
∗
H
i
j
x
+
B
∗
H
i
j
x
(50)
u
x
,
t
≃
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
t
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
t
B
∗
φ
x
-
x
j
2
∂
∂
x
u
x
,
t
≃
∂
∂
x
u
~
x
,
t
=
∂
∂
x
∑
j
=
1
N
i
α
j
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
∂
∂
x
Φ
∗
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
∂
∂
x
B
∗
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
Φ
∗
∂
∂
x
φ
x
-
x
j
2
+
∑
j
=
N
i
+
1
N
α
j
B
∗
∂
∂
x
φ
x
-
x
j
2Now, we have Φ
∗
=
∂
/
∂
x and with the absence of the derivative operation on the boundary, B
∗, the remaining term of the above equation can be written as(51)
∂
∂
x
u
x
,
t
≃
∂
∂
x
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
∂
∂
x
∂
φ
x
-
x
j
2
∂
x
+
∑
j
=
N
i
+
1
N
α
j
∂
∂
x
φ
x
-
x
j
2
=
∑
j
=
1
N
i
α
j
x
i
-
x
j
2
x
-
x
j
2
2
d
2
φ
d
r
2
+
y
i
-
y
j
3
x
-
x
j
2
3
d
φ
d
r
+
∑
j
=
N
i
+
1
N
α
j
x
i
-
x
j
r
d
φ
d
r
=
∑
j
=
1
N
i
α
j
x
i
-
x
j
2
x
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
y
i
-
y
j
3
x
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
x
i
-
x
j
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2Similarly, (52)
∂
∂
y
u
x
,
t
≃
∂
∂
y
u
~
x
,
t
=
∑
j
=
1
N
i
α
j
y
i
-
y
j
2
x
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
x
i
-
x
j
3
x
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
y
i
-
y
j
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2By substituting (49), (51), and (52) back into the governing equations(53)
∂
u
∂
t
=
1
R
e
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
-
u
∂
u
∂
x
-
v
∂
u
∂
yTherefore, for each center x
k, the following is reached:(54)
∂
u
∂
t
k
=
1
R
e
∑
j
=
1
N
i
α
j
2
ε
2
exp
⁡
-
ε
2
x
k
-
x
j
2
2
8
ε
6
x
k
-
x
j
2
4
-
24
ε
4
x
k
-
x
j
2
2
+
6
ε
2
+
2
2
ε
2
exp
⁡
-
ε
2
x
k
-
x
j
2
2
6
ε
2
-
4
ε
4
x
k
-
x
j
2
2
-
1
x
k
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
2
ε
x
k
-
x
j
2
2
-
1
+
1
x
k
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
2
ε
x
-
x
j
2
2
-
1
+
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
-
u
k
∑
j
=
1
N
i
α
j
x
i
-
x
j
2
x
k
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
2
ε
x
k
-
x
j
2
2
-
1
+
y
i
-
y
j
3
x
k
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
x
i
-
x
j
-
2
ε
2
exp
⁡
-
ε
x
-
x
j
2
2
-
v
k
∑
j
=
1
N
i
α
j
y
i
-
y
j
2
x
k
-
x
j
2
2
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
2
ε
x
k
-
x
j
2
2
-
1
+
x
i
-
x
j
3
x
k
-
x
j
2
2
-
2
ε
2
exp
⁡
-
ε
x
k
-
x
j
2
2
+
∑
j
=
N
i
+
1
N
α
j
y
i
-
y
j
-
2
ε
2
e
x
p
-
ε
x
-
x
j
2
2The whole process explained above is applied to the other equation of the coupled governing equations. Note that all the terms with1
/
x
k
-
x
j
2 shall be cancelling each other out so there is no zero-division problem in computing process.


## 3.4. Computing Setup and Algorithm

Each computation under this work was completed following the computing steps detailed below.Step 1.
ChooseN collocation or center nodes X
c
=
x
j
j
=
1
N, on the domain Ω.Step 2.
Specify the desired values of(2.1)
The Reynolds number, Re.
(2.2)
The final timet.
(2.3)
The time stepΔ
t.
(2.4)
The Gaussian shape parameterε.Step 3.
Compute the initial collocation matricesΑ
~ and Β
~ using the initial Hermite approximation equations (20)-(21).Step 4.
Apply the initial conditions,U
0
,
V
0, to get α and β from (30), (31), (33), (34), (35), (36), (37), (38), (40), (41), (42), and (43).Step 5.
Compute the solutions for the next time step via the time-discretization explained above.Step 6.
Construct the collocation matricesΑ
~ and Β
~ in the forms expressed in (30), (31), (33), (34), (35), (36), (37), (38), (40), (41), (42), and (43) using the solutions previously obtained in Step 5.Step 7.
Carry on Steps5 and 6 until reaching the final time step t.All computing experiments carried out in this study were executed on a computer notebook: Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz with 8.00 GB of RAM and 64-bit operating system.


## 4. Numerical Experiments and General Discussion

The test case chosen for numerical investigation carried out in this study is the most popular form of Burgers’ equations where the corresponding exact solutions for validation are provided by using a Hopf-Cole transformation nicely carried out by Fletcher [4] in 1983 and are expressed as follows:(55)
u
x
,
y
,
t
=
3
4
-
1
4
1
+
exp
⁡
-
4
x
+
4
y
-
t
R
e
32
-
1
(56)
v
x
,
y
,
t
=
3
4
+
1
4
1
+
exp
⁡
-
4
x
+
4
y
-
t
R
e
32
-
1where both the initial and the boundary conditions to be imposed to the equation system are generated directly from the above exact forms over the domain Ω
=
(
x
,
y
)
:
0
≤
x
≤
1,0
≤
y
≤
1. With their rich of challenging features, the equations have received interest from several researchers and been treated with variety of numerical techniques as mentioned in Section 1. Nevertheless, it can be clearly seen that, with the reasons previously mentioned, most studies are concerned only with only cases of relatively low Reynolds numbers, i.e., 1
≤
R
e
≤
500. In this work, the investigation covers a wide range of Reynolds number from low to moderate R
e
∈
1,800 and, moreover, goes above R
e
=
1,000, where no numerical works have ever reached before.Solutions obtained from applying the Hermite collocation method are validated by comparing to both those documented in literatures and those produced by the exact formula. Error measurements are carried out using the following error norms:(1)
Maximum error (L
∞);(57)
L
∞
=
m
a
x
1
≤
i
≤
N
u
~
x
i
-
u
x
i
(2)
Root-mean-square error (L
r
m
s);(58)
L
r
m
s
=
1
N
∑
j
=
1
N
u
~
x
j
-
u
x
j
2
1
/
2In addition to these norms, the results presented below are sometimes referred to as “Average-Error” and it is to be understood as(59)
A
v
e
r
a
g
e
-
L
∞
=
L
∞
U
+
L
∞
V
2where L
∞
U and L
∞
V are of U- and V-velocity component, respectively.The whole study is split into two main parts, one concerned with low- to moderate-Reynolds number cases,1
≤
R
e
≤
500, and that with relatively high value of Reynolds number, 500
<
R
e
≤
1,300. The so-called Reynolds number plays crucial rules in determining the ratio of forces exerted in the system and this is remained as one of the most challenging tasks to simulate and mimic both numerically and experimentally.In order to demonstrate the overall effectiveness of the Hermite type of collocation scheme, investigations using the conventional global collocation known as “Kansa” method, (14), are also studied parallelly. The results obtained from using the Hermite type and Kansa are referred to hereafter as “Hermite” and “Conventional Kansa/Kansa,” respectively.


### 4.1. Low- to Moderate-Reynolds Number Cases

As it is known to be one of the most challenging tasks for any collocation method, finding a suitable shape parameter,ε, is to be firstly presented. In the past, some attempts to pinpoint the optimal value of ε involve the classic work of Hardy [35] where it was shown that, by fixing the shape at ε
=
1
/
0.815
d, where d
=
1
/
N
∑
i
=
1
N
d
i and d
i is the distance from the node to its nearest neighbor, good results should be anticipated. Also, in the work of Franke and Schaback [36] where the choice of a fixed shape of the form ε
=
0.8
N
/
D with D is the diameter of the smallest circle containing all data nodes, can also be a good alternative. In 2000, Zhang et al. [37] demonstrated and concluded that the optimal shape parameter is highly problem dependent. In 2002, Wang and Lui [38] pointed out that, by analyzing the condition number of the collocation matrix, a suitable range of derivable values of ε can be found. Later in 2003, Lee et al. [39] suggested that the final numerical solutions obtained are found to be less affected by the method when the approximation is applied locally rather than globally. A rather recent work is the selection of an interval for variable shape parameter by Biazar and Hosami in 2016 [40], where a novel algorithm for determining and interval was proposed.It is to be noted, however, that even though this topic has long been widely studied (see also in [41–44]), the majority of studies are done on “Multiquadric RBF” or “Inverse-Multiquadric RBF,” whereas works on Gaussian type are rather rare. A confirmation on problem-dependent nature of Gaussian RBF shape value was confirmed in the work of Davydov and Oanh [45]. A more recent work is the application of the traditional Kansa collocation to linear PDE based on a variable shape parameter called symmetric variable shape parameter (SVSP) carried out by Ranjbar [46].Even though Gaussian RBF was focused on in these two later works, the scheme does not contain any self-adjoint operator. Moreover, the complications introduced by the structure of the governing equation at hand are known to greatly affect the final choice of the shape. Taking these facts into consideration, with the transient, nonlinear, and coupling nature of Burgers’ equations, this work chose to directly measure the overall errors generated in the system. To complete this task, a large number of numerical experiments were carried out in this work aiming to reach the optimal value ofε. The impact of shape parameter on the quality of the computed solution is demonstrated in Figures 1–4. Figures 1 and 2 show that, within the range of 0
<
ε
≤
1.5, both methods reveal slight differences in the trends of L
r
m
s, with some strong fluctuations being found for the case of Kansa. For this particular set of conditions, at Re = 10 with t
=
0.5
,
Δ
t
=
0.005
,
N
=
10
×
10, it can be noted that a more suitable shape parameter may be found at ε
<
0.1 for the Hermite scheme. This observation can also be useful for Kansa case as indicated in the Figure. Figures 3 and 4 indicate clearly that both of the methods become less sensitive to the shape value when it goes beyond ε
≈
300. Above this point, it should be observed that L
r
m
s remains under 10
E
-
01 when using Hermite scheme whereas L
r
m
s stays constantly above 10
E
-
01.Figure 1
Effect of the Gaussian shape parameter (ε
≤
1.5) on the result accuracy for the Hermite collocation scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 2
Effect of the Gaussian shape parameter (ε
≤
1.5) on the result accuracy for the conventional Kansa scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 3
Effect of the Gaussian shape parameter (ε
≫
1.5) on the result accuracy for the Hermite collocation scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 4
Effect of the Gaussian shape parameter (ε
≫
1.5) on the result accuracy for the conventional Kansa scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Information provided previously on the effect of shape parameter leads to decision made for the case with Re = 50. At randomly chosen points, Table1 provides the solutions of both U- and V-velocity components compared against the corresponding exact values. It can be seen that the solutions are in good agreement for both choices of shape parameter; ε
=
0.001 and ε
=
0.1. A bigger picture of results obtained in this work is shown in Table 2 and Table 3, with cases of Reynolds number starting from 1 up to 400.Table 1
Numerical solutions for Re = 50,Δ
t
=
0.001, t
=
0.5, and N
=
12
×
12.
x
y
ε
=
0.001
ε
=
0.1
U-velocity
V-velocity
U-velocity
V-velocity
Hermite
Exact
Hermite
Exact
Hermite
Exact
Hermite
Exact
0
0.3
0.6873
0.6872
0.8127
0.8128
0.6873
0.6872
0.8127
0.8128
0.1
0.2
0.6653
0.6628
0.8347
0.8372
0.6553
0.6631
0.8247
0.8367
0.3
0.8
0.7381
0.7395
0.7719
0.7605
0.7281
0.7398
0.7719
0.7601
0.4
0.6
0.6938
0.6943
0.8162
0.8057
0.6838
0.6946
0.8162
0.8052
0.6
0.7
0.6553
0.6628
0.8347
0.8372
0.6653
0.6631
0.8247
0.8367
0.7
0.2
0.5149
0.5105
0.9951
0.9895
0.5049
0.5109
0.9951
0.9890
0.7
0.6
0.5892
0.5872
0.9008
0.9128
0.5892
0.5875
0.9108
0.9124
0.8
0.1
0.5014
0.5031
0.9986
0.9969
0.5014
0.5034
0.9986
0.9964
1
0.4
0.5027
0.5027
0.9973
0.9973
0.5027
0.5027
0.9973
0.9973
1
0.7
0.5164
0.5164
0.9836
0.9836
0.5164
0.5164
0.9836
0.9836Table 2
AveragedL
r
m
s-error obtained with N
=
100 and ε
=
0.001 for low to moderate Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
1
100
5E-04
1.11E-06
1.56E-05
3.55E-05
3.12E-03
10
100
0.0050
1.08E-03
1.58E-04
5.65E+26
2.91E-02
200
900
0.0010
8.97E-03
1.19E-03
3.66E-01
1.32E-01
400
900
0.01
1.55E-02
1.77E-03
4.48E-01
1.39E-01Table 3
AveragedL
∞-error obtained with N
=
100, ε
=
0.001 for low to moderate Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
1
100
5E-04
2.10E-06
1.95E-05
7.51E-05
3.91E-03
10
100
0.0050
2.20E-03
1.95E-04
1.10E+27
3.87E-02
200
900
0.0010
2.38E-02
3.90E-03
8.36E-01
2.49E-01
400
900
0.01
4.11E-02
7.80E-03
1.03E+00
2.50E-01From these two Tables, it can be clearly seen that, with different values of time stepΔ
t and the number of collocation nodes N, both methods produce approximately the same quality of results at both points of measured times: t
=
0.01 and t
=
2.0.For a validation purpose, a comparison of computed solution obtained from using Hermite to other numerical works available in literature is done and presented in Tables4 and 5. In this case, the number of nodes is 441 and the measurement is taken at Re = 100, Δ
t
=
0.0001, and t
=
2.0. The comparison clearly shows that good agreement in terms of solution quality is reached using the method under investigation.Table 4
Computed value of U-velocity at for Re = 100,Δ
t
=
0.0001, t
=
2.0, and N
=
21
×
21
.
x
y
Exact
Hermite
[20]
[8]
[21]
[5]
0.1
0.1
0.500482
0.502321
0.500470
0.50035
0.50012
0.49983
0.3
0.3
0.500482
0.504311
0.500441
0.50042
0.50042
0.49977
0.5
0.5
0.500482
0.502553
0.500414
0.50046
0.50041
0.49973
0.3
0.7
0.555675
0.552258
0.554805
0.55609
0.55587
0.55429
0.1
0.9
0.744256
0.741302
0.744197
0.74409
0.74416
0.74340
0.5
0.9
0.555675
0.556811
0.554489
0.55604
0.55637
0.55413Table 5
Computed value of V-velocity at Re = 100 withΔ
t
=
0.0001, t
=
2.0, and N
=
21
×
21
.
x
y
Exact
Hermite
[20]
[8]
[21]
[5]
0.1
0.1
0.999518
0.99695
0.999530
0.99936
0.99946
0.99826
0.3
0.3
0.999518
0.99771
0.999559
0.99951
0.99938
0.99861
0.5
0.5
0.999518
0.99882
0.999586
0.99958
0.99941
0.99821
0.3
0.7
0.944325
0.94502
0.945195
0.94387
0.94387
0.94409
0.1
0.9
0.755744
0.75499
0.755803
0.75592
0.75558
0.75500
0.5
0.9
0.944325
0.94102
0.945511
0.94392
0.94345
0.94345At Reynolds number of 500, solutions produced witht
=
0.5
,
Δ
t
=
0.001
,
N
=
30
×
30 are depicted in Figures 5 and 6. Once again, reasonably good quality of computed solutions is obtained.Figure 5
U-velocity profile produced by (left) the Hermite collocation scheme and (right) the exact, at Re = 500 witht
=
0.5
,
Δ
t
=
0.001, and N
=
30
×
30.Figure 6
V-velocity profile Average produced by (left) the Hermite collocation scheme and (right) The exact, at Re = 500 witht
=
0.5
,
Δ
t
=
0.001, and N
=
30
×
30.For all the information presented so far, it strongly indicates that Hermite is now successfully applied to this type of PDEs. This gives some confidence to move to more challenging cases with higher Reynolds number and the results are to be presented in the following sections.


### 4.2. Moderate- to High-Reynolds Number Cases

It is well known that the problem becomes more difficult when the Reynolds number increases as the instability in the numerical process can completely ruin the final solution [47]. This problem is generally encountered by any study that concerns Reynolds number above 800. In this work, the numerical investigation starts from Re = 600 and moves to the highest of 1,300. With using only 900 collocation nodes and a fixed value of shape ε
=
0.001, Tables 6 and 7, respectively, show L
r
m
s and L
∞ of solutions computed by Hermite and Kansa at two different times t
=
0.01 and t
=
2.0. At the smaller time t
=
0.01, both tables reveal that both Kansa and Hermite produce approximately the same results quality with those obtained from Hermite being only slightly better by an order of magnitude. This is no longer the case when the computation reaches the time of t
=
2.0 where it can clearly been seen from both tables that Kansa has a strong growth in errors with L
r
m
s
≈ 1.0E-02 growing to L
r
m
s
≈ 1.0E+04 and L
∞
≈ 1.0E-02 to L
∞
≈ 1.0E+04. However, the Hermite loses only slightly its accuracy with L
r
m
s
≈ 1.0E-03 to L
r
m
s
≈ 1.0E-01 and L
∞
≈ 1.0E-02 to L
∞
≈ 1.0E-01.Table 6
AveragedL
r
m
s- error obtained with ε
=
0.001 for moderate to high Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
600
900
0.001
1.52E-02
2.27E-03
4.23E-01
1.41E-01
800
900
0.001
4.94E-02
2.87E-03
1.04E+04
1.42E-01
1000
900
0.001
4.88E-02
3.54E-03
1.11E+04
1.42E-01
1200
900
0.001
4.87E-02
4.21E-03
1.20E+04
1.42E-01Table 7
AveragedL
∞- error obtained with ε
=
0.001 for moderate to high Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
600
900
0.001
4.37E-02
1.17E-02
9.55E-01
2.50E-01
800
900
0.001
1.50E-01
1.55E-02
2.87E+04
2.50E-01
1000
900
0.001
1.52E-01
1.94E-02
3.09E+04
2.50E-01
1200
900
0.001
1.55E-01
2.32E-02
3.34E+04
2.50E-01The influence of the shape parameter to the solution accuracy is further studied also for these cases of high Reynolds number and the general observation can be made based on information given in Tables8 and 9. At a very small number of collocation nodes, N
=
100, and a fixed time step Δ
t
=
0.001, both tables reveal the same conclusion as previously shown. Kansa is capable of reproducing good results quality only for the case with small time, i.e., t
=
0.01. At t
=
2.0, on the other hand, while Kansa is seen to fail to provide good solutions, Hermite clearly manages to keep the quality of numerical solutions around 1.00E-01 for both error norms and for all values of shapes under investigation, ε
∈
0.01,0.1,1.0. At a particular Re = 1,000, a wider range of shape parameter, 0
<
ε
≤
3, is further investigated and the measurement of L
r
m
s is shown in Figures 7 and 8 where a clear observation on the effect of shape parameter is illustrated. The figures indicate clearly that the Hermite scheme is much less sensitive to shape parameter when the Reynolds number increases when compared with the conventional Kansa method. In addition, it is also observed that the Hermite scheme studied and applied in this work is comparatively less sensitive to the progression of time even when the Reynolds number is taken up to 1,300. This observation is illustrated in Figures 9 and 10 with surface plots of solutions obtained from Kansa and Hermite being depicted in Figure 11 together with the exact ones. It must be noted here that results of U-velocity component are not shown here as they are in similar trend to those for V-velocity component.Table 8
AveragedL
r
m
s- error obtained with N
=
100, Δ
t
=
0.001 for high Reynolds number.
ε
Re
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
0.01
1000
6.70E-03
5.49E-03
4.22E+01
1.22E-01
0.1
1100
7.55E-03
6.02E-03
2.90E+02
1.22E-01
1
1200
3.91E-02
6.58E-03
2.27E+08
1.85E-01Table 9
AveragedL
∞- error obtained with N
=
100
,
Δ
t
=
0.001 for high Reynolds number.
ε
Re
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
0.01
1000
2.49E-02
1.94E-02
1.31E+02
2.50E-01
0.1
1100
2.91E-02
2.13E-02
1.03E+03
2.50E-01
1
1200
1.35E-01
2.38E-02
8.82E+08
5.86E-01Figure 7
AverageL
r
m
s-error obtained from using Kansa collocation scheme at Re = 1,000 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 8
AverageL
r
m
s- error obtained from using Hermite collocation scheme at Re = 1,000 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 9
L
∞-error progression in time 0.00
≤
t
≤
2.00of the U-maximum error comparison computed at Re = 1,300 with ε
=
0.005
,
Δ
t
=
0.0005, and N
=
10
×
10.Figure 10
L
∞-error progression in time 0.00
≤
t
≤
2.00of the V-maximum error comparison computed at Re = 1,300 with ε
=
0.005
,
Δ
t
=
0.0005, and N
=
10
×
10.Figure 11
V-velocity profile at Re = 1,300 withε
=
0.005
,
t
=
0.5
,
Δ
t
=
0.0005
,
N
=
30
×
30; (a) Conventional Kansa, (b) Hermite scheme, and (c) exact solution.
(a)
(b)
(c)


## 4.1. Low- to Moderate-Reynolds Number Cases

As it is known to be one of the most challenging tasks for any collocation method, finding a suitable shape parameter,ε, is to be firstly presented. In the past, some attempts to pinpoint the optimal value of ε involve the classic work of Hardy [35] where it was shown that, by fixing the shape at ε
=
1
/
0.815
d, where d
=
1
/
N
∑
i
=
1
N
d
i and d
i is the distance from the node to its nearest neighbor, good results should be anticipated. Also, in the work of Franke and Schaback [36] where the choice of a fixed shape of the form ε
=
0.8
N
/
D with D is the diameter of the smallest circle containing all data nodes, can also be a good alternative. In 2000, Zhang et al. [37] demonstrated and concluded that the optimal shape parameter is highly problem dependent. In 2002, Wang and Lui [38] pointed out that, by analyzing the condition number of the collocation matrix, a suitable range of derivable values of ε can be found. Later in 2003, Lee et al. [39] suggested that the final numerical solutions obtained are found to be less affected by the method when the approximation is applied locally rather than globally. A rather recent work is the selection of an interval for variable shape parameter by Biazar and Hosami in 2016 [40], where a novel algorithm for determining and interval was proposed.It is to be noted, however, that even though this topic has long been widely studied (see also in [41–44]), the majority of studies are done on “Multiquadric RBF” or “Inverse-Multiquadric RBF,” whereas works on Gaussian type are rather rare. A confirmation on problem-dependent nature of Gaussian RBF shape value was confirmed in the work of Davydov and Oanh [45]. A more recent work is the application of the traditional Kansa collocation to linear PDE based on a variable shape parameter called symmetric variable shape parameter (SVSP) carried out by Ranjbar [46].Even though Gaussian RBF was focused on in these two later works, the scheme does not contain any self-adjoint operator. Moreover, the complications introduced by the structure of the governing equation at hand are known to greatly affect the final choice of the shape. Taking these facts into consideration, with the transient, nonlinear, and coupling nature of Burgers’ equations, this work chose to directly measure the overall errors generated in the system. To complete this task, a large number of numerical experiments were carried out in this work aiming to reach the optimal value ofε. The impact of shape parameter on the quality of the computed solution is demonstrated in Figures 1–4. Figures 1 and 2 show that, within the range of 0
<
ε
≤
1.5, both methods reveal slight differences in the trends of L
r
m
s, with some strong fluctuations being found for the case of Kansa. For this particular set of conditions, at Re = 10 with t
=
0.5
,
Δ
t
=
0.005
,
N
=
10
×
10, it can be noted that a more suitable shape parameter may be found at ε
<
0.1 for the Hermite scheme. This observation can also be useful for Kansa case as indicated in the Figure. Figures 3 and 4 indicate clearly that both of the methods become less sensitive to the shape value when it goes beyond ε
≈
300. Above this point, it should be observed that L
r
m
s remains under 10
E
-
01 when using Hermite scheme whereas L
r
m
s stays constantly above 10
E
-
01.Figure 1
Effect of the Gaussian shape parameter (ε
≤
1.5) on the result accuracy for the Hermite collocation scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 2
Effect of the Gaussian shape parameter (ε
≤
1.5) on the result accuracy for the conventional Kansa scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 3
Effect of the Gaussian shape parameter (ε
≫
1.5) on the result accuracy for the Hermite collocation scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 4
Effect of the Gaussian shape parameter (ε
≫
1.5) on the result accuracy for the conventional Kansa scheme computed at Re = 10 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Information provided previously on the effect of shape parameter leads to decision made for the case with Re = 50. At randomly chosen points, Table1 provides the solutions of both U- and V-velocity components compared against the corresponding exact values. It can be seen that the solutions are in good agreement for both choices of shape parameter; ε
=
0.001 and ε
=
0.1. A bigger picture of results obtained in this work is shown in Table 2 and Table 3, with cases of Reynolds number starting from 1 up to 400.Table 1
Numerical solutions for Re = 50,Δ
t
=
0.001, t
=
0.5, and N
=
12
×
12.
x
y
ε
=
0.001
ε
=
0.1
U-velocity
V-velocity
U-velocity
V-velocity
Hermite
Exact
Hermite
Exact
Hermite
Exact
Hermite
Exact
0
0.3
0.6873
0.6872
0.8127
0.8128
0.6873
0.6872
0.8127
0.8128
0.1
0.2
0.6653
0.6628
0.8347
0.8372
0.6553
0.6631
0.8247
0.8367
0.3
0.8
0.7381
0.7395
0.7719
0.7605
0.7281
0.7398
0.7719
0.7601
0.4
0.6
0.6938
0.6943
0.8162
0.8057
0.6838
0.6946
0.8162
0.8052
0.6
0.7
0.6553
0.6628
0.8347
0.8372
0.6653
0.6631
0.8247
0.8367
0.7
0.2
0.5149
0.5105
0.9951
0.9895
0.5049
0.5109
0.9951
0.9890
0.7
0.6
0.5892
0.5872
0.9008
0.9128
0.5892
0.5875
0.9108
0.9124
0.8
0.1
0.5014
0.5031
0.9986
0.9969
0.5014
0.5034
0.9986
0.9964
1
0.4
0.5027
0.5027
0.9973
0.9973
0.5027
0.5027
0.9973
0.9973
1
0.7
0.5164
0.5164
0.9836
0.9836
0.5164
0.5164
0.9836
0.9836Table 2
AveragedL
r
m
s-error obtained with N
=
100 and ε
=
0.001 for low to moderate Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
1
100
5E-04
1.11E-06
1.56E-05
3.55E-05
3.12E-03
10
100
0.0050
1.08E-03
1.58E-04
5.65E+26
2.91E-02
200
900
0.0010
8.97E-03
1.19E-03
3.66E-01
1.32E-01
400
900
0.01
1.55E-02
1.77E-03
4.48E-01
1.39E-01Table 3
AveragedL
∞-error obtained with N
=
100, ε
=
0.001 for low to moderate Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
1
100
5E-04
2.10E-06
1.95E-05
7.51E-05
3.91E-03
10
100
0.0050
2.20E-03
1.95E-04
1.10E+27
3.87E-02
200
900
0.0010
2.38E-02
3.90E-03
8.36E-01
2.49E-01
400
900
0.01
4.11E-02
7.80E-03
1.03E+00
2.50E-01From these two Tables, it can be clearly seen that, with different values of time stepΔ
t and the number of collocation nodes N, both methods produce approximately the same quality of results at both points of measured times: t
=
0.01 and t
=
2.0.For a validation purpose, a comparison of computed solution obtained from using Hermite to other numerical works available in literature is done and presented in Tables4 and 5. In this case, the number of nodes is 441 and the measurement is taken at Re = 100, Δ
t
=
0.0001, and t
=
2.0. The comparison clearly shows that good agreement in terms of solution quality is reached using the method under investigation.Table 4
Computed value of U-velocity at for Re = 100,Δ
t
=
0.0001, t
=
2.0, and N
=
21
×
21
.
x
y
Exact
Hermite
[20]
[8]
[21]
[5]
0.1
0.1
0.500482
0.502321
0.500470
0.50035
0.50012
0.49983
0.3
0.3
0.500482
0.504311
0.500441
0.50042
0.50042
0.49977
0.5
0.5
0.500482
0.502553
0.500414
0.50046
0.50041
0.49973
0.3
0.7
0.555675
0.552258
0.554805
0.55609
0.55587
0.55429
0.1
0.9
0.744256
0.741302
0.744197
0.74409
0.74416
0.74340
0.5
0.9
0.555675
0.556811
0.554489
0.55604
0.55637
0.55413Table 5
Computed value of V-velocity at Re = 100 withΔ
t
=
0.0001, t
=
2.0, and N
=
21
×
21
.
x
y
Exact
Hermite
[20]
[8]
[21]
[5]
0.1
0.1
0.999518
0.99695
0.999530
0.99936
0.99946
0.99826
0.3
0.3
0.999518
0.99771
0.999559
0.99951
0.99938
0.99861
0.5
0.5
0.999518
0.99882
0.999586
0.99958
0.99941
0.99821
0.3
0.7
0.944325
0.94502
0.945195
0.94387
0.94387
0.94409
0.1
0.9
0.755744
0.75499
0.755803
0.75592
0.75558
0.75500
0.5
0.9
0.944325
0.94102
0.945511
0.94392
0.94345
0.94345At Reynolds number of 500, solutions produced witht
=
0.5
,
Δ
t
=
0.001
,
N
=
30
×
30 are depicted in Figures 5 and 6. Once again, reasonably good quality of computed solutions is obtained.Figure 5
U-velocity profile produced by (left) the Hermite collocation scheme and (right) the exact, at Re = 500 witht
=
0.5
,
Δ
t
=
0.001, and N
=
30
×
30.Figure 6
V-velocity profile Average produced by (left) the Hermite collocation scheme and (right) The exact, at Re = 500 witht
=
0.5
,
Δ
t
=
0.001, and N
=
30
×
30.For all the information presented so far, it strongly indicates that Hermite is now successfully applied to this type of PDEs. This gives some confidence to move to more challenging cases with higher Reynolds number and the results are to be presented in the following sections.


## 4.2. Moderate- to High-Reynolds Number Cases

It is well known that the problem becomes more difficult when the Reynolds number increases as the instability in the numerical process can completely ruin the final solution [47]. This problem is generally encountered by any study that concerns Reynolds number above 800. In this work, the numerical investigation starts from Re = 600 and moves to the highest of 1,300. With using only 900 collocation nodes and a fixed value of shape ε
=
0.001, Tables 6 and 7, respectively, show L
r
m
s and L
∞ of solutions computed by Hermite and Kansa at two different times t
=
0.01 and t
=
2.0. At the smaller time t
=
0.01, both tables reveal that both Kansa and Hermite produce approximately the same results quality with those obtained from Hermite being only slightly better by an order of magnitude. This is no longer the case when the computation reaches the time of t
=
2.0 where it can clearly been seen from both tables that Kansa has a strong growth in errors with L
r
m
s
≈ 1.0E-02 growing to L
r
m
s
≈ 1.0E+04 and L
∞
≈ 1.0E-02 to L
∞
≈ 1.0E+04. However, the Hermite loses only slightly its accuracy with L
r
m
s
≈ 1.0E-03 to L
r
m
s
≈ 1.0E-01 and L
∞
≈ 1.0E-02 to L
∞
≈ 1.0E-01.Table 6
AveragedL
r
m
s- error obtained with ε
=
0.001 for moderate to high Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
600
900
0.001
1.52E-02
2.27E-03
4.23E-01
1.41E-01
800
900
0.001
4.94E-02
2.87E-03
1.04E+04
1.42E-01
1000
900
0.001
4.88E-02
3.54E-03
1.11E+04
1.42E-01
1200
900
0.001
4.87E-02
4.21E-03
1.20E+04
1.42E-01Table 7
AveragedL
∞- error obtained with ε
=
0.001 for moderate to high Reynolds number.
Re
N
Δ
t
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
600
900
0.001
4.37E-02
1.17E-02
9.55E-01
2.50E-01
800
900
0.001
1.50E-01
1.55E-02
2.87E+04
2.50E-01
1000
900
0.001
1.52E-01
1.94E-02
3.09E+04
2.50E-01
1200
900
0.001
1.55E-01
2.32E-02
3.34E+04
2.50E-01The influence of the shape parameter to the solution accuracy is further studied also for these cases of high Reynolds number and the general observation can be made based on information given in Tables8 and 9. At a very small number of collocation nodes, N
=
100, and a fixed time step Δ
t
=
0.001, both tables reveal the same conclusion as previously shown. Kansa is capable of reproducing good results quality only for the case with small time, i.e., t
=
0.01. At t
=
2.0, on the other hand, while Kansa is seen to fail to provide good solutions, Hermite clearly manages to keep the quality of numerical solutions around 1.00E-01 for both error norms and for all values of shapes under investigation, ε
∈
0.01,0.1,1.0. At a particular Re = 1,000, a wider range of shape parameter, 0
<
ε
≤
3, is further investigated and the measurement of L
r
m
s is shown in Figures 7 and 8 where a clear observation on the effect of shape parameter is illustrated. The figures indicate clearly that the Hermite scheme is much less sensitive to shape parameter when the Reynolds number increases when compared with the conventional Kansa method. In addition, it is also observed that the Hermite scheme studied and applied in this work is comparatively less sensitive to the progression of time even when the Reynolds number is taken up to 1,300. This observation is illustrated in Figures 9 and 10 with surface plots of solutions obtained from Kansa and Hermite being depicted in Figure 11 together with the exact ones. It must be noted here that results of U-velocity component are not shown here as they are in similar trend to those for V-velocity component.Table 8
AveragedL
r
m
s- error obtained with N
=
100, Δ
t
=
0.001 for high Reynolds number.
ε
Re
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
0.01
1000
6.70E-03
5.49E-03
4.22E+01
1.22E-01
0.1
1100
7.55E-03
6.02E-03
2.90E+02
1.22E-01
1
1200
3.91E-02
6.58E-03
2.27E+08
1.85E-01Table 9
AveragedL
∞- error obtained with N
=
100
,
Δ
t
=
0.001 for high Reynolds number.
ε
Re
t
=
0.01
t
=
2.0
Kansa
Hermite
Kansa
Hermite
0.01
1000
2.49E-02
1.94E-02
1.31E+02
2.50E-01
0.1
1100
2.91E-02
2.13E-02
1.03E+03
2.50E-01
1
1200
1.35E-01
2.38E-02
8.82E+08
5.86E-01Figure 7
AverageL
r
m
s-error obtained from using Kansa collocation scheme at Re = 1,000 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 8
AverageL
r
m
s- error obtained from using Hermite collocation scheme at Re = 1,000 with t
=
0.5
,
Δ
t
=
0.005, and N
=
10
×
10.Figure 9
L
∞-error progression in time 0.00
≤
t
≤
2.00of the U-maximum error comparison computed at Re = 1,300 with ε
=
0.005
,
Δ
t
=
0.0005, and N
=
10
×
10.Figure 10
L
∞-error progression in time 0.00
≤
t
≤
2.00of the V-maximum error comparison computed at Re = 1,300 with ε
=
0.005
,
Δ
t
=
0.0005, and N
=
10
×
10.Figure 11
V-velocity profile at Re = 1,300 withε
=
0.005
,
t
=
0.5
,
Δ
t
=
0.0005
,
N
=
30
×
30; (a) Conventional Kansa, (b) Hermite scheme, and (c) exact solution.
(a)
(b)
(c)


## 5. Conclusion

In this investigation, a further developed version of Kansa collocation method called “Hermite” scheme was applied to one of the most complex types of PDEs. To evaluate the general effectiveness of the mesh-free method, the Burgers’ equations are chosen. The Gaussian type of radial basis function was chosen for this task. With applying the differential operator twice to the governing equation before being approximated via a linear summation of the basis function, as in the Hermite scheme, the solutions obtained have revealed some interesting aspects and some findings can be listed as follows:(i)
Judging from the low Reynolds number cases, the Hermite is successfully applied to the transient, nonlinear, and coupled-Burgers’ equations with a good agreement in results’ accuracy, compared with the exact solution and other numerical works.
(ii)
It is found that, with using the Hermite method, the computing process is much less sensitive to the change of the shape parameterε making it more practical when in use.
(iii)
By using the Hermite scheme, the growth of error accumulated in time is found to be significantly smaller than the conventional Kansa method particularly at high Reynolds number.
(iv)
With a suitable value of shape parameterε, at moderate to very high Reynolds number, the numerical method is found to produce reasonable results.
(v)
Across the whole range of Reynolds number under investigation,1
≤
R
e
≤
1,300, the Hermite scheme is found to outperform the conventional version of collocation or Kansa method.Even though some desirable aspects have been discovered in this study regarding the use of collocation mesh-free methods, there are still some interesting questions to be further addressed and investigated, for instance, the process of choosing the optimal shape parameter where several factors shall be taken into consideration: the effect of viscous force, the number of centers, the time-dependent physics, and so on. This is all subject to further studies of the authors.



---

*Source: 3416860-2018-09-27.xml*