MLPG method based on moving kriging interpolation for solving convection–diffusion equations with integral condition

S. Khankham; A. Luadsong; N. Aschariyaphotha

doi:10.1016/j.jksus.2015.03.001

View/Download PDF

Buy Reprints

PDF

Translate this page into:

Original article

27 (

); 292-301

doi:

10.1016/j.jksus.2015.03.001

MLPG method based on moving kriging interpolation for solving convection–diffusion equations with integral condition

S. Khankham^{^a,^b}, A. Luadsong^{^a,⁎}, N. Aschariyaphotha^{^c}

a Department of Mathematics, Faculty of Science, King Mongkut’s University of Technology Thonburi (KMUTT), 126 Pracha-utid Road, Bangmod, Toongkru, Bangkok 10140, Thailand

b Science Achievement Scholarship of Thailand in Mathematics, Faculty of Science at Mahidol University, 272 Rama VI Road, Ratchathewi, Bangkok 10400, Thailand

c Ratchaburi Learning Park, King Mongkut’s University of Technology Thonburi (KMUTT), Rang Bua, Chom Bueng, Ratchaburi 70150, Thailand

⁎Corresponding author. anirut.lua@kmutt.ac.th (A. Luadsong)

Received: 2014-6-16, Accepted: 2015-3-2, Published: 2015-10-01

Disclaimer:
This article was originally published by Elsevier and was migrated to Scientific Scholar after the change of Publisher.

Peer review under responsibility of King Saud University.

Abstract

A formulation of the meshless local Petrov–Galerkin (MLPG) method based on the moving kriging interpolation (MK) is presented in this paper. The method is used for solving time-dependent convection–diffusion equations in two-dimensional spaces with the Dirichlet, Neumann, and non-local boundary conditions on a square domain. The method is developed based on the moving kriging interpolation method for constructing shape functions which have the Kronecker delta property. In the method, the test function in each sub-domain is chosen as the indicator function. The Crank–Nicolson method is chosen for temporal discretization. Two test problems are presented which demonstrate the easiness and accuracy of this method as shown by the relative error.

Keywords

Convection–diffusion equations

Crank–Nicolson

Moving kriging interpolation

MLPG

Integral condition

Show Related Articles from PubMed

1 Introduction

In recent decades, meshless method in mathematical modeling has attracted much attention, due to their flexibility in solving practical science and engineering problems. The meshless method has been developed and achieved remarkable progress in mathematical modeling and related fields. The main reason for interest in meshless method is it can save computational time because no mesh is needed. They can easily produce high accuracy, since in areas where more refinement is required; nodes can quite easily be added. In addition, they can easily construct high-order shape functions. They can easily solve large deformation and strong nonlinear problem. The connectivity among the nodes is generated as a part of the computation and can change with time.

Based on the form of system equation, the meshfree methods are classified into three categories. The first category covers the meshfree methods based on strong forms of system equations, in which discretization is performed directly from the governing differential equations, such as the general finite difference method (Liszka and Orkisz, 1980), the smooth particle hydrodynamic (SPH) method (Liu and Liu, 2003; Lucy, 1977), and other meshfree collocation methods. The second category includes meshfree methods based on weak forms of system equations, such as the element-free Galerkin (EFG) method (Belytschko et al., 1994), the meshless local Petrov–Galerkin (MLPG) method (Atluri and Zhu, 1998; Atluri and Lin, 2000; Abbasbandy and Shirzadi, 2011; Mirzaei, 2010; Mahmoodabadi et al., 2011 ), the point interpolation method (PIM) (Liu, 2002; Liu and Gu, 2001), etc. The third category concerns meshfree methods based on the combination of weak and strong forms, such as the meshfree weak–strong-form (MWS) method. The major difference in these meshless methods is in the choice of interpolation techniques.

The MLPG method was first discovered by Atluri and Zhu (1998). Meshless local Petrov–Galerkin (MLPG) method is a truly meshless method. This method is based on a weak form computed over a local sub-domain and the moving least squares (MLS) approximation. The MLPG method is one of the most viable methods in which the MLS approach is used to construct the shape functions. Although the MLPG method has been applied to many problems, there exists an inconvenience or disadvantage when using the MLPG because of the difficulty in implementing essential boundary conditions. This is because the MLS shape functions lack the Kronecker delta property. Therefore the MK interpolation method (Yimnak and Luadsong, 2014; Kaewumpai and Luadsong, 2014) has been proposed to overcome this problem. It uses the nodal values in the local support domain to construct shape functions with the Kronecker delta property. The MK interpolation method works well for practical problems.

Finding numerical solutions with non-local boundary conditions is important in applications such as chemical diffusion, heat conduction processes, population dynamics thermoelasticity, medical science, electrochemistry and control theory. Because of the complex problems, some analytic solutions are difficult to find in the sense of boundary conditions and geometrical shapes. So, recently much attention has been paid in the literature to the analysis and implementation of accurate methods for the numerical solution of time-dependent partial differential equations with non-local boundary conditions. Abbasbandy and Shirzadi (2010) researched on the MPLG method for the two-dimensional diffusion equation with the Neumann boundary condition and non-classical boundary condition, and a meshless method for the two-dimensional diffusion equation with an integral condition. Mohyud-Din and Yildirim (2010) employed the homotopy analysis method (HAM) for the solutions of two-dimensional diffusion equations subject to non-standard boundary specifications. The HAM method does not need transformation techniques, linearization and discretization. The approximate solutions by the HAM method were compared with exact solutions. It was shown that the method is reliable, efficient and requires less computation. Sajavicius (2013) researched optimization, conditioning and accuracy of radial basis function method for partial differential equations with nonlocal boundary conditions in case of two-dimensional Poisson equation. Techapirom and Luadsong (2013) have used an MLPG method to study the two-dimensional heat equation with Dirichlet, Neumann and non-local boundary conditions in a square domain. Their study demonstrated the good accuracy of the proposed method and indicates that it can be easily extrapolated to other problems. Sataprahm and Luadsong (2014) have used the Meshless Local Petrov–Galerkin (MLPG) method for the incompressible Navier–Stokes equations. The numerical examples presented the local symmetric weak form (LSWF) and the local unsymmetric weak form (LUSWF) with a classical Gaussian weight and an improved Gaussian weight on both regular and irregular nodes is demonstrated. It is found that LSWF1 with a classical Gaussian weight order 2 gives the most accurate result. Goodrich (2015) considers the coupled systems of boundary value problems with nonlocal boundary conditions. The result shows that this system can possess at least one positive solution even if no growth conditions. Islam et al. (2015) present two numerical methods that are analyzed for the solution of two-dimensional Poisson equation with two different types of nonlocal boundary conditions. The first numerical method is a collocation method based on Haar wavelet whereas the second numerical method is a meshless method based on different types of radial basis functions (RBFs). A two-point boundary condition and an integral boundary condition are the two types of nonlocal boundary conditions considered in the present work. Two numerical results have shown that the accuracy and efficiency wise performance is confirmed through application of the algorithms on the benchmark tests.

The purpose of this article is to present a very efficient MLPG based on MK interpolation for solving the two-dimensional time-dependent convection–diffusion equation with integral conditions. To deal with time derivatives, using the Crank–Nicolson implicit method. The solution is approximated by the MK interpolation method which holds the Kronecker delta property, thereby enhancing the nodal shape construction accuracy. It was shown that the method is reliable and efficient. Two test problems are presented illustrating the easiness and accuracy of this method as shown by the relative error.

2 Moving kriging interpolation method

The kriging interpolation is well-known geostatic technique for spatial interpolation in geology and mining (Lucy, 1977). The formulation of the construction of meshless shape function by MK interpolation is introduced briefly in the following. Similar to the MLS approximation, consider the function $T (x)$ defined in the domain $Ω$ discretized by a set of properly scattered nodes $x_{i}, i = 1, 2, \dots, N$ , where $N$ is the total number of nodes in the whole domain. It is assumed that only $N$ nodes surrounding point $x$ have the effect on $T (x)$ . The sub-domain $Ω_{x}$ that encompasses these surrounding nodes is called the interpolation domain of point $x$ . The MK interpolation $T^{h} (x)$ at point $x$ is defined as presented in Liu and Gu (2001). Therefore, the formulation of the meshless shape function using MK interpolation is given by:

(1)

T^{h} (x) = \sum_{i = 1}^{N} ϕ_{i} (x) T_{i} = Φ (x) T, x \in Ω_{x}

where

T = {[T (x_{1}) T (x_{2}) \dots T (x_{N})]}^{T}

is a vector value of the function in the domain

Ω

Φ (x)

is a

1 \times N

vector of shape functions, expressed as:

(2)

Φ (x) = p^{T} (x) A + r^{T} (x) B,

where matrices A and B are defined as:

(3)

A = {(P^{T} R^{- 1} P)}^{- 1} P^{T} R^{- 1}

(4)

B = R^{- 1} (I - PA)

In which I is a unit matrix of size N × N, and vector $p (x)$ is:

(5)

p^{T} (x) = [p_{1} (x) p_{2} (x) \dots p_{N} (x)] .

For matrix P with the size N × M, values of the polynomial basis functions at the given set of nodes are collected:

(6)

P = [\begin{matrix} p_{1} (x_{1}) & \dots & p_{M} (x_{1}) \\ ⋮ & ⋱ & ⋮ \\ p_{1} (x_{N}) & \dots & p_{M} (x_{N}) \end{matrix}] .

Matrices R and vector $r (x)$ are defined by the following equations:

(7)

R = [\begin{matrix} r (x_{1}, x_{1}) & \dots & r (x_{1}, x_{N}) \\ ⋮ & ⋱ & ⋮ \\ r (x_{N}, x_{1}) & \dots & r (x_{N}, x_{N}) \end{matrix}]

(8)

r^{T} (x) = [r (x, x_{1}) r (x, x_{2}) \dots r (x, x_{N})],

where

r (x_{i}, x_{j})

is the correlation function between any pair of nodes located at X_i and X_j representing the covariance of the field value

T (x)

, i.e.

(9)

r (x_{i}, x_{j}) = E [T (x_{i}) T (x_{j})] .

Similarly, the covariance $E [T (x_{i}) T (x_{j})]$ can be replaced by $r (x_{i}, x_{j})$ . It can be seen from the foregoing formulations that the values of matrices R and $r (x)$ play important roles in the computation. A simple and frequently- used correlation function is a Gaussian function:

(10)

r (x_{i}, x_{j}) = e^{- α_{c} {(\frac{r_{ij}}{d_{c}})}^{2}}

where

r_{ij} = ‖ x_{i} - x_{j} ‖, d_{c}

and

α_{c} > 0

are the correlation parameters used to fit the model and are assumed to be given.

The first-order partial derivatives of the shape function $Φ (x)$ against the coordinates $x_{i}, i = 1, 2$ can be easily obtained from Eq. (2)

(11)

Φ_{, i} (x) = p_{, i}^{T} (x) A + r_{, i}^{T} (x) B,

where

{(\cdot)}_{i}

denotes

\partial (\cdot) / \partial x^{i}

3 Formulation of the moving kriging interpolation for convection–diffusion equation

3.1

3.1 Governing equation and boundary conditions

The two-dimensional time-dependent convection–diffusion equation subject to non-local boundary conditions can be written as:

(12)

\frac{\partial T}{\partial t} + u \cdot \nabla T - \nabla \cdot (K \nabla T) = f, x \in Ω, t > 0

where

T

is the unknown function,

K

is the diffusivity coefficient,

u

is the flow velocity,

f

is a source term. The boundary conditions are assumed to be:

–

The non-local boundary condition:
(13)
$T = h_{0} (x) μ (t), x \in Γ_{n} .$
–

The essential boundary condition:
(14)
$T = h_{1} (x, t), x \in Γ_{d} .$
–

The natural boundary condition:
(15)
$n \cdot \nabla T = g (x, t), x \in Γ_{q} .$
–

The initial condition assumed to be
(16)
$T (x, 0) = T_{0} (x), x \in Ω,$ and the integral condition
(17)
$\int_{Ω} T (x, t) d Ω = m (t),$

where

μ (t)

is the unknown function,

g, h_{0}, h_{1}

and

m

are given functions. The non-local boundary condition is variable-separable, with spatial dependence given by

h_{0} (x)

and time dependence given by

μ (t)

n

is the outward unit normal vector to

\partial Ω

, where

\partial Ω

is the boundary of the domain.

Γ_{n}, Γ_{d}

and

Γ_{q}

are the part of the boundary

\partial Ω

, satisfying

Γ_{n} \cap Γ_{d} \cap Γ_{q} = ϕ

and

Γ_{n} \cup Γ_{d} \cup Γ_{q} = \partial Ω

3.2

3.2 Local weak-forms

In implementation of the present MK interpolation method, a weak form is first constructed over a local sub-domain $Ω_{s}$ bounded by $\partial Ω_{s}$ , where $Ω_{s}$ is located entirely inside the global domain $Ω$ . Using the local weighted residual technique, the local weak form of Eq. (12) over a local sub-domain $Ω_{s}$ can be written as:

(18)

\int_{Ω_{s}^{i}} [\frac{\partial T}{\partial t} + u \cdot \nabla T - \nabla \cdot (K \nabla T)] v_{i} d Ω = \int_{Ω_{s}^{i}} {fv}_{i} d Ω

where

T

is the trial function and

v

is the test function.

Using $\nabla \cdot (K \nabla T) v_{i} = \nabla \cdot (K \nabla (T) v_{i}) - K \nabla T \cdot \nabla v_{i}$ and the divergence theorem, we obtain the following local weak formulation (LWF):

(19)

\int_{Ω_{s}^{i}} \frac{\partial T}{\partial t} v_{i} d Ω + \int_{Ω_{s}^{i}} u \cdot \nabla {Tv}_{i} d Ω - \int_{\partial Ω_{s}^{i}} K v_{i} \frac{\partial T}{\partial n} d Γ + \int_{Ω_{s}^{i}} K \nabla T \cdot \nabla v_{i} d Ω = \int_{Ω_{s}^{i}} {fv}_{i} d Ω .

In general, the boundary $\partial Ω_{s}^{i}$ of the local sub-domain $Ω_{s}^{i}$ may intersect with the boundary of the global domain $Ω$ . Therefore, $\partial Ω_{s}^{i} = Γ_{sn}^{i} \cup Γ_{sd}^{i} \cup Γ_{sq}^{i} \cup L_{s}^{i}$ where, $L_{s}^{i}$ is a part of $\partial Ω_{s}^{i}$ . By imposing the natural boundary condition, we obtain:

(20)

\int_{Ω_{s}^{i}} \frac{\partial T}{\partial t} v_{i} d Ω + \int_{Ω_{s}^{i}} u \cdot \nabla {Tv}_{i} d Ω - \int_{Γ_{sn}^{i}} K v_{i} \frac{\partial T}{\partial n} d Γ - \int_{Γ_{sd}^{i}} K v_{i} \frac{\partial T}{\partial n} d Γ - \int_{L_{s}^{i}} K v_{i} \frac{\partial T}{\partial n} d Γ + \int_{Ω_{s}^{i}} K \nabla T \cdot \nabla v_{i} d Ω = \int_{Ω_{s}^{i}} {fv}_{i} d Ω + \int_{Γ_{sq}^{i}} {gv}_{i} d Γ .

The test function in each sub-domain is chosen as the indicator function.

(21)

v_{i} (x) = \{\begin{matrix} 1, x \in Ω_{s}^{i} \\ 0, x \notin Ω_{s}^{i} \end{matrix}

so,

\nabla v (x) = 0, \forall x \in Ω_{s}^{i}

and the local weak form Eq. (22) is transformed into the following simple local integral equation:

(22)

\int_{Ω_{s}^{i}} \frac{\partial T}{\partial t} d Ω + \int_{Ω_{s}^{i}} u \cdot \nabla Td Ω - \int_{Γ_{sn}^{i}} K \frac{\partial T}{\partial n} d Γ - \int_{Γ_{sd}^{i}} K \frac{\partial T}{\partial n} d Γ - \int_{L_{s}^{i}} K \frac{\partial T}{\partial n} d Γ = \int_{Ω_{s}^{i}} fd Ω . + \int_{Γ_{sq}^{i}} g d Γ .

3.3

3.3 Discretization of the weak form

As a known test function is used in the LWF, the use of the LWF for one point (and hence for one local domain) will yield only one linear equation involving $\hat{T} (t)$ . Note that the trial function $T$ within the sub-domain $Ω_{s}$ , in the MK approximation, is determined by the fictitious nodal values ${\hat{T}}_{i} (t),$ within the domain of definition for all points $x$ falling within $Ω_{s}$ . The LWF in Eq. (22) gives one algebraic equation relating all these ${\hat{T}}_{i} (t)$ . Thus, we obtain as many equations as the number of nodes. Therefore, we need as many local domains $Ω_{s}$ as the number of nodes in the global domain to obtain as many equations as the number of unknowns. In the present implementation, the local domain is chosen as a circle, centered at a node $x_{i}$ .

To obtain the discrete equations from the LWF in Eq. (22), the MK approximation in Eq. (1), is alternatively used to approximate the trial function $T$ . Substituting of Eq. (1) into the LWF Eq. (22) for the entire node, leads to the following discretized system of linear equations:

(23)

\sum_{j = 1}^{N} (\int_{Ω_{s}^{i}} ϕ_{j} (x) d Ω)) {\hat{T}}_{j}^{'} (t) - \sum_{j = 1}^{N} [(\int_{Ω_{s}^{i}} u \cdot \nabla ϕ_{j} (x) d Ω - \int_{L_{s}^{i}} K ϕ_{j, n} (x) d Γ - \int_{Γ_{sn}^{i}} K ϕ_{j, n} (x) d Γ - \int_{Γ_{sd}^{i}} K ϕ_{j, n} (x) d Γ) {\hat{T}}_{j} (t)] = \int_{Ω_{s}^{i}} fd Ω + \int_{Γ_{sq}^{i}} g d Γ .

where

i = 1, 2, 3, \dots, N

Eq. (23) can be written in the matrix form as follows:

(24)

C \frac{\partial T}{\partial t} + DT = F (t),

where

C, D

and

F

are matrices described as follows:

(25)

C = [C_{ij}], C_{ij} = \int_{Ω_{s}^{i}} ϕ_{j} (x) d Ω

(26)

D = [D_{ij}], D_{ij} = \int_{Ω_{s}^{i}} u \cdot \nabla ϕ_{j} (x) d Ω - \int_{L_{s}^{i}} K ϕ_{j, n} (x) d Γ - \int_{Γ_{sn}^{i}} K ϕ_{j, n} (x) d Γ - \int_{Γ_{sd}^{i}} K ϕ_{j, n} (x) d Γ

(27)

F = [F_{i}], F_{i} = \int_{Ω_{s}^{i}} f d Ω + \int_{Γ_{sq}^{i}} g d Γ .

Setting a time-stepping scheme to overcome the time derivative and applying the Crank–Nicolson technique of approximation to Eq. (24) yields:

(28)

(2 C + Δ t D) T^{k + 1} = (2 C + Δ t D) T^{k} + 2 Δ t F^{\frac{1}{2}}

where

T

is a matrix described as follows:

(29)

T = (\begin{matrix} \hat{T} \\ \hat{μ} \end{matrix})

Assuming that ${\hat{T}}_{i}^{k}$ , for $i = 1, 2, . . ., N$ and ${\hat{μ}}^{k},$ are known, our aim is to compute ${\hat{T}}_{i}^{k + 1},$ for $i = 1, 2, \dots, N$ and ${\hat{μ}}^{k + 1}$ . Now, we have $N + 1$ unknowns, so to compute these unknowns we need one equation, which can be obtained from the non-local boundary condition from Eq. (17). Substituting trial function Eq. (1), we obtain:

(30)

\int_{Ω} {(T^{h})}^{k + 1} (x) d Ω = \int_{Ω} \sum_{j = 1}^{N} ϕ_{j} (x) {\hat{T}}_{j}^{k + 1} d Ω = \sum_{j = 1}^{N} [\int_{Ω} ϕ_{j} (x) d Ω] {\hat{T}}_{j}^{k + 1} = m^{k + 1}

which can be written in a matrix form as:

(31)

S {\hat{T}}^{k + 1} = m^{k + 1},

where

S

is a matrix described as follows:

(32)

S = [S_{j}], S_{j} = \int_{Ω} ϕ_{j} (x) d Ω .

For nodes on the natural, essential and non-local boundaries, we present an algorithm to impose boundary conditions as follows:

–

For node $x_{i}$ is on the natural boundary condition, $x_{i} \in Γ_{q}$
(33)
$\sum_{j = 1}^{N} [\int_{Γ_{sq}^{i}} ϕ_{j, n} (x) d Γ] {\hat{T}}_{j}^{k + 1} = \int_{Γ_{sq}^{i}} g (x^{l}, (k + 1) Δ t) d Γ$
–

For node $x_{i}$ is on the essential boundary condition, $x_{i} \in Γ_{d}$
(34)
$\sum_{j = 1}^{N} [\int_{Γ_{sd}^{i}} ϕ_{j} (x) d Γ] {\hat{T}}_{j}^{k + 1} = \int_{Γ_{sd}^{i}} h_{1} (x^{l}, (k + 1) Δ t) d Γ$
–

For node $x_{i}$ is on the non-local boundary condition, $x_{i} \in Γ_{n}$
(35)
$\sum_{j = 1}^{N} [\int_{Γ_{sn}^{i}} ϕ_{j} (x) d Γ] {\hat{T}}_{j}^{k + 1} - \int_{Γ_{sn}^{i}} h_{0} (x^{l}) d Γ {\hat{μ}}^{k + 1} = 0 .$

4 Numerical examples

In the computation, problems are considered in a unit square domain given by $0 ⩽ x ⩽ 1 .$ . The flow is unidirectional and constant with velocity components $u_{1} = 1$ and $u_{2} = 1, u = {(u_{1}, u_{2})}^{T}$ . The diffusivity coefficient is 1. The radius of the local sub-domain for each node is chosen as $r_{0} = 0.9 h$ , with $h$ is the distance between two consecutive nodes in each direction. The results are obtained for $T$ and $μ$ for several time levels with $Δ t = 0.001$ . The numbers of complete monomial basis (M) are 3, 6 and 10 and the numbers of nodal points (N) are 25, 36 and 81. The exact solutions are known for these problems and are used to test the accuracy of these numerical schemes. The numerical results in these tables show that the method converges by increasing the number of nodal points, and increasing the number of monomial basis for a fixed number of nodal points enhance accuracy.

4.1

4.1 Problem 1

We consider the two-dimensional time-dependent convection–diffusion equation with Dirichlet and non-local boundary condition: $\frac{\partial T}{\partial t} + u \cdot \nabla T - \nabla \cdot (K \nabla T) = f, (x, t) \in [0, 1] \times [0, 1] \times [0, T_{\max}], x= (x, y)$ subject to initial, Dirichlet boundary conditions and non-local boundary conditions: $T (x, 0) = (1 + y) \exp (x), m (t) = \frac{3}{2} \exp (t) (\exp (1) - 1), 0 \leq x \leq 1,$ $h_{0} (x) = \exp (x), h_{2} (x, t) = 2 \exp (x + t), t \geq 0 .$ $h_{3} (x, t) = (1 + y) \exp (t), h_{4} (x, t) = (1 + y) \exp (1 + t), t \geq 0, h_{1} (x, t) = h_{2} (x, t) \cup h_{3} (x, t) \cup h_{4} (x, t)$

The exact solution of this problem is: $T (x, t) = (1 + y) \exp (x + t), 0 \leq x \leq 1, t \geq 0, μ (t) = \exp (t) .$

The source term of this problem is: $f (x) = (1 + y) \exp (x + t) + u_{1} (1 + y) \exp (x + t) + u_{2} \exp (x + t) - (1 + y) \exp (x + t), 0 \leq x \leq 1 .$ The numerical results are shown in Tables 1 and 2, Figs. 1–12 .

Table 1 The relative error for

T (x, t)

at M = 3, 6 and 10, N = 25, 36 and 81.

N/m	3	6	10
25	2.0163e−2	3.1541e-3	2.5064e-3
36	1.7916e-2	2.6715e-3	1.9803e-3
81	6.2939e-3	2.5497e-3	1.5738e-3

Table 2 The relative error for

μ (t)

at M = 3, 6 and 10, N = 25, 36 and 81.

N/m	3	6	10
25	8.8847e-3	1.0669e-3	1.0410e-3
36	7.5965e-3	7.2043e-4	6.6892e-4
81	2.5943e-3	6.5504e-4	5.9344e-4

The relative error for T ( x , t ) at M = 3, 6 and 10 with N = 25.

The relative error for T ( x , t ) at N = 25, 36 and 81 with M = 3.

The relative error for T ( x , t ) at M = 3, 6 and 10 with N = 36.

The relative error for T ( x , t ) at N = 25, 36 and 81 with M = 6.

The relative error for T ( x , t ) at M = 3, 6 and 10 with N = 81.

The relative error for T ( x , t ) at N = 25, 36 and 81 with M = 10.

The relative error for μ ( t ) at M = 3, 6 and 10 with N = 25.

The relative error for μ ( t ) at N = 25, 36 and 81 with M = 3.

The relative error for μ ( t ) at M = 3, 6 and 10 with N = 36.

The relative error for μ ( t ) at N = 25, 36 and 81 with M = 6.

The relative error for μ ( t ) at M = 3, 6 and 10 with N = 81.

The relative error for μ ( t ) at N = 25, 36 and 81 with M = 10.

4.2

4.2 Problem 2

We consider the two-dimensional time-dependent convection–diffusion equation with Dirichlet, Neumann and non-local boundary condition: $\frac{\partial T}{\partial t} + u \cdot \nabla T - \nabla \cdot (K \nabla T) = f, (x, t) \in [0, 1] \times [0, 1] \times [0, T_{\max}], x = (x, y)$ subject to initial, Dirichlet boundary condition, Neumann boundary condition and non-local boundary condition: $T (x, 0) = \exp (x + y + 2 t), m (t) = \exp (2 t) (\exp (2) - 2 \exp (1) + 1), 0 \leq x \leq 1,$ $h_{0} (x) = \exp (x), h_{1} (x, t) = \exp (1 + x + 2 t), t \geq 0 .$ $g_{0} (x, t) = \exp (y + 2 t), g_{1} (x, t) = \exp (1 + y + 2 t), t \geq 0, g (x, t) = g_{0} (x, t) \cup g_{1} (x, t)$

The exact solution of this problem is: $T (x, t) = \exp (x + y + 2 t), 0 \leq x \leq 1, t \geq 0, μ (t) = \exp (2 t) .$

The source term of this problem is: $f (x) = u_{1} \exp (x + y + 2 t) + u_{2} \exp (x + y + 2 t), 0 \leq x ⩽ 1 .$ The numerical results are shown in Tables 3 and 4, Figs. 13–24 .

Table 3 The relative error for

T (x, t)

at M = 3, 6 and 10, N = 25, 36 and 81.

N/m	3	6	10
25	1.9372e-2	5.3279e-3	4.4636e-3
36	1.6748e-2	4.8862e-3	2.7373e-3
81	5.8101e-3	2.8265e-3	1.5231e-3

Table 4 The relative error for

μ (t)

at M = 3, 6 and 10, N = 25, 36 and 81.

N/m	3	6	10
25	6.1612e-3	1.0348e-3	1.0079e-3
36	2.5005e-3	5.6160e-4	5.4160e-4
81	1.7075e-3	3.1049e-4	3.0039e-4

5 Conclusions

A formulation of meshless local Petrov–Galerkin (MLPG) has been presented in the present work. We have proposed an idea of MK interpolation demonstrating how it can be used to construct the meshless shape functions for solving the two-dimensional time-dependent convection–diffusion equation with integral conditions. To deal with time derivatives, the Crank–Nicolson implicit method was used. The key feature of the MK interpolation is that the shape function possesses the Kronecker delta property which makes it easy to directly impose Dirichlet boundary conditions. The numerical results show a decrease in relative error with increasing number of nodal points. Similarly, the relative error decreases when the size of the monomial basis increases. Moreover, it can be concluded that the numerical results of two test problems have been used to verify the efficiency, easiness and accuracy of the method.

Acknowledgements

This research is partially supported by the Science Achievement Scholarship of Thailand in Mathematics.

References

Abbasbandy S., Shirzadi A., . A meshless method for two-dimensional diffusion equation with an integral condition. Eng. Anal. Bound. Elem.. 2010;34:1031-1037.
[Google Scholar]
Abbasbandy S., Shirzadi A., . MLPG method for two-dimensional diffusion equation with Neumann’s and non-classical boundary conditions. Appl. Numer. Math.. 2011;61:170-180.
[Google Scholar]
Atluri S.N., Lin H., . Meshless local Petrov–Galerkin (MLPG) method for convection–diffusion problems. CMES. 2000;1(2):45-60.
[Google Scholar]
Atluri S.N., Zhu T.L., . A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput. Mech.. 1998;22:117-127.
[Google Scholar]
Belytschko T., Lu Y.Y., Gu L., . Element-free Galerkin methods. Int. J. Numer. Methods Eng.. 1994;37(1):229-256.
[Google Scholar]
Goodrich C.S., . Coupled systems of boundary value problems with nonlocal boundary conditions. Appl. Math. Lett.. 2015;41:17-22.
[Google Scholar]
Islam S.U., Aziz I., Ahmad M., . Numerical solution of two-dimensional elliptic PDEs with nonlocal boundary conditions. Comput. Math. Appl.. 2015;69:180-205.
[Google Scholar]
Kaewumpai S., Luadsong A., . Two-field-variable meshless method based on moving kriging interpolation for solving simply supported thin plates under various loads. J. King Saud Univ. Sci. 2014:1018-3647.
[Google Scholar]
Liszka T., Orkisz J., . The finite difference method on arbitrary irregular grid and its application in applied mechanics. Comput. Struct.. 1980;11:83-95.
[Google Scholar]
Liu G.R., . A point assembly method for stress analysis for two-dimensional solids. Int. J. Solids Struct.. 2002;39:261-276.
[Google Scholar]
Liu G.R., Gu Y.T., . A point interpolation method for two-dimensional solid. Int. J. Numer. Methods Eng.. 2001;50:937-951.
[Google Scholar]
Liu G.R., Liu M.B., . Smoothed Particle Hydrodynamics – A Meshfree Practical Method. Singapore: World Scientific; 2003.
Lucy L., . A numerical approach to testing the fission hypothesis. Astrophys. J.. 1977;82:1013-1024.
[Google Scholar]
Mahmoodabadi M.J., Abedzadeh R.M., Bagheri A., Baradaran G.H., . Meshless local Petrov–Galerkin method for 3D steady-state heat conduction problems. Adv. Mech. Eng. 2011 (10 pages)
[CrossRef] [Google Scholar]
Mirzaei D., . Meshless local Petrov–Galerkin (MLPG) approximation to the two dimensional sine-Gordor equation. Appl. Math. Comput.. 2010;233:2737-2754.
[Google Scholar]
Mohyud-Din S., Yildirim A., . On two-dimensional diffusion with integral condition. J. King Saud Univ. Sci.. 2010;23(2):121-125.
[Google Scholar]
Sajavicius S., . Optimization, conditioning and accuracy of radial basis function method for partial differential equations with nonlocal boundary condition-A case of two-dimensional Poisson equation. Eng. Anal. Bound. Elem.. 2013;37:788-804.
[Google Scholar]
Sataprahm C., Luadsong A., . The meshless local Petrov–Galerkin method for simulating unsteady incompressible fluid flow. J. Egypt. Math. Soc.. 2014;23(3):501-510.
[Google Scholar]
Techapirom T., Luadsong A., . The MLPG with improved weight function for two-dimensional heat equation with non-local boundary condition. J. King Saud Univ. Sci.. 2013;25:341-348.
[Google Scholar]
Yimnak K., Luadsong A., . A local integral equation formulation based on moving kriging interpolation for solving coupled nonlinear reaction–diffusion equations. Adv. Math. Phys. 2014
[Google Scholar]