Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 2, pp. 571-582, Warsaw 2017
DOI: 10.15632/jtam-pl.55.2.571

MESHLESS LOCAL RADIAL POINT INTERPOLATION (MLRPI) FOR
GENERALIZED TELEGRAPH AND HEAT DIFFUSION EQUATION WITH

NON-LOCAL BOUNDARY CONDITIONS

Elyas Shivanian, Arman Khodayari

Imam Khomeini International University, Department of Mathematics, Qazvin, Iran

e-mail address: shivanian@sci.ikiu.ac.ir

In this paper, themeshless local radial point interpolation (MLRPI)method is formulated to
the generalized one-dimensional linear telegraph and heat diffusion equation with non-local
boundary conditions. TheMLRPImethod is categorized under meshless methods in which
any background integration cells are not required, so that all integrations are carried out
locally over small quadrature domains of regular shapes, such as lines in one dimensions,
circles or squares in two dimensions and spheres or cubes in three dimensions. A technique
based on the radial point interpolation is adopted to construct shape functions, also called
basis functions, using the radial basis functions. These shape functions have delta function
property in the framework of interpolation, therefore they convince us to impose boundary
conditions directly. The time derivatives are approximated by the finite difference time-
-steppingmethod.We also apply Simpson’s integration rule to treat the non-local boundary
conditions. Convergency and stability of theMLRPImethod are clarified by surveying some
numerical experiments.

Keyword: non-local boundary condition, meshless local radial point interpolation (MLRPI)
method, local weak formulation, radial basis function, telegraph equation

1. Introduction

The telegraph equation is one of the important equations of mathematical physics with ap-
plications to many different fields such as transmission and propagation of electrical signals
(Gonzalez-Velasco, 1995; Jordan and Puri, 1999), vibrational systems (Boyce and DiPrima,
1977), randomwalk theory (Banasiak andMika, 1998) andmechanical systems (Tikhonov and
Samarskii, 1990), etc. The heat diffusion andwave propagation equations are particular cases of
the telegraph equation. Recently, increasing attention has been paid to the development, ana-
lysis and implementation of stable methods for numerical solutions of second-order hyperbolic
equations. There have beenmany numericalmethods for hyperbolic equations, such as the finite
difference, the finite element, and the collocationmethods, etc. (seeAlmenar et al., 1997; Ciment
and Leventhal, 1978) and literatures therein.

On the other hand, many of natural phenomena in science and engineering have been mo-
delled by non-local boundary value problems. In these non-local problems, some integral terms
often appear in the boundary conditions. These types of problems constitute a special class of
boundary value problems which widely appear in mathematical modelling of various processes
of physics, heat transfer, ecology, thermoelasticity, chemistry, biology and industry.

According to the numerical results obtained, the present methods can be considered as
practical andeffective numerical techniques to solve telegraph equationswithnon-local boundary
conditions.



572 E. Shivanian, A. Khodayari

Let Ω = [0,1]. Consider the 1D linear telegraph equation

∂2u

∂t2
+ c

∂u

∂t
+ bu−p

∂2u

∂x2
= f(x,t) (x,t)∈ Ω × [0,T ] (1.1)

with the initial and non-local boundary conditions

u(x,0)= u0(x)
∂u

∂t
(x,0)= ψ(x)

u(0, t) = γ1

1∫

0

u(x,t) dx+µ1(t)

u(1, t) = γ2

1∫

0

u(x,t) dx+µ2(t)

(1.2)

where c, b and p are positive constants, γ1 and γ2 are constants and the functions f, ψ, µ1(t)
and µ2(t) are assumed to be sufficiently smooth. Many partial differential equations are too
complex to be solved by analyticalmethods. This causedmathematicians and engineers to come
up with numerical methods such as the finite difference method (FDM) and the finite element
method (FEM) to solve the equations. Although, these methods have been successfully applied
to computational fluid dynamics problems, their accuracy depends critically on mesh quality
and they have many difficulties in dealing with some complex problems. These difficulties can
be overcome bymeshless methods which have attracted considerable interest over the past few
years (Kochmann and Venturini, 2014; Liu and Gu, 2005; Pan and Yuan, 2009; Sladek et al.,
2006).Thesemeshlessmethodsdonot requiremesh fordiscretisation of theproblemdomain, and
they construct approximate functions only via a set of nodes, so-called field nodes. In general,
the meshless methods can be grouped into two categories. The first category is based on weak
forms such as the element free Galerkin (EFG) method (Belytschko et al., 1995; Singh et al.,
2007), the second category is based on strong forms such as meshless methods based on the
radial basis functions (RBFs) (Dehghan and Shokri, 2008; Kansa, 1990). In addition, ameshless
methodbased on combination of the strongandweak formhas also beendeveloped and is known
as the meshless weak strong (MWS) form method. Due to the ill-conditioning of the resultant
linear systems in the RBF-collocation method, various approaches are proposed to circumvent
this problem (Libre et al., 2008; Ling and Schaback, 2008), being among them. Theweak forms
are used to derive a set of algebraic equations through a numerical integration process using
a set of quadrature domain that may be constructed globally or locally in the domain of the
problem. In the global formulation, background cells are required for the integration of theweak
form. Strictly speaking, thesemeshlessmethods are not trulymeshlessmethods.But inmethods
based on the local weak form formulation, numerical integrations are carried out over a local
quadrature domains, therefore, no cells are required. As a result, they are referred to as truly
meshlessmethods such as themeshless local Petrov-Galerkin (MLPG)method (Atluri and Zhu,
1998; Dehghan and Mirzaei, 2008; Shirzadi, 2014; Shivanian, 2015b). In the literature, several
meshless weak form methods have been proposed such as the diffuse element method (DEM)
(Nayroles et al., 1992), smooth particle hydrodynamic (SPH) (Bratsos, 2008; Dashtimanesh and
Ghadimi, 2013), reproducing kernel particle method (RKPM) (Liu et al., 1995), boundary node
method (BNM) (Mukherjee and Mukherjee, 1997), partition of unity finite element method
(PUFEM) (Melenh and Babuska, 1996), finite sphere method (FSM) (De and Bathe, 2000),
boundary point interpolation method (BPIM) (Gu and Liu, 2002) and boundary radial point
interpolation method (BRPIM) (Gu and Liu, 2003). Liu applied the concept of MLPG and
developed themeshless local radial point interpolation (MLRPI) method (Hosseini et al., 2015;
Liu and Gu, 2001; Shivanian, 2013, 2015a; Shivanian and khodabandehlo, 2014). In this paper,



Meshless local radial point interpolation (MLRPI) for generalized telegraph... 573

we concentrate on the numerical solution of Eqs. (1.1) and (1.2) using themeshless local radial
point interpolation (MLRPI) method. Besides, we use Simpson’s integration rule to impose the
non-local boundary condition.

2. Approximation of field variables using the radial point interpolation method

Consider a continuous function u(x) defined in a domain Ω, which is represented by a set of
field nodes. The u(x) at the point of interest x is approximated as follows

u(x)=
n∑

i=1

Ri(x)ai+
m∑

j=1

pj(x)bj =R
T(x)a+PT(x)b (2.1)

whereRi(x) is thea radial basis function (RBF), n is the number ofRBFs, pj(x) is themonomial
in the 1-D space x and m is the number of themonomials. In the present work, we have applied
thin plate spline (TPS) multiquadrics (MQ) as the radial basis functions in Eq. (2.1). In order
to determine ai and bj in Eq. (2.1), a support domain is needed for the point of interest at x so
that n field nodes are included in the support domain. Then, coefficients ai and bj in Eq. (2.1)
can be determined by the following system of n linear equations

Us =Rna+Pmb (2.2)

in which the vector Us is

Us = {u1,u2,u3, . . . ,un}
T (2.3)

moreover, Rn and Pm are the RBFs and polynomial moment matrices, respectively. On the
other hand, Eq. (2.1) can be rewritten as

u(x) =RT(x)a+PT(x)b=
{
RT(x),PT(x)

}[a
b

]
(2.4)

and then, by using that, we obtain

u(x) =
{
RT(x),PT(x)

}[Rn Pm
PTm 0

]
−1

Ũs

=
{
RT(x),PT(x)

}
G−1Ũs = Φ̃

T(x)Ũs

(2.5)

where Φ̃T(x) can be be introduced by

Φ̃T(x)=
{
RT(x),PT(x)

}
G−1 = {φ1(x),φ2(x), . . . ,φn(x),φn+1(x), . . . ,φn+m(x)} (2.6)

The first n functions of the above vector function are called the RPIM shape functions corre-
sponding to the nodal displacements.We show them by the vector Φ̃T(x), so that it is

Φ̃T(x)= {φ1(x),φ2(x), . . . ,φn(x)} (2.7)

Equation (2.5) is then transformed into

u(x) = Φ̃T(x)Us =
n∑

i=1

φi(x)ui (2.8)



574 E. Shivanian, A. Khodayari

3. Finite differences approximation

The followingfinitedifferenceapproximations of theorderO(∆t)2 areused for timediscretization

∂2u(x, t)

∂t2
∼=
1

∆t2
(
u(k+1)(x)−2u(k)(x)+u(k−1)(x)

)

∂u(x, t)

∂t
∼=
1

2∆t

(
u(k+1)(x)−u(k−1)(x)

)
(3.1)

Also, we employ the following approximation using the Crank-Nicolson technique

u(x, t)∼=
1

3

(
u(k+1)(x)+u(k)(x)+u(k−1)(x)

)

∂2u(x, t)

∂x2
∼=
1

3

(∂2u(k+1)(x, t)
∂x2

+
∂2u(k)(x, t)

∂x2
+
∂2u(k−1)(x, t)

∂x2

) (3.2)

where uk(x)= u(x,k∆t).

Using the above approximations, Eq. (1.1) can be written as

1

∆t2
(
u(k+1)(x)−2u(k)(x)+u(k−1)(x)

)
+

c

2∆t

(
u(k+1)(x)−u(k−1)(x)

)

+
b

3

(
u(k+1)(x)+u(k)(x)+u(k−1)(x)

)

−
p

3

(∂2u(k+1)(x)
∂x2

+
∂2u(k)(x)

∂x2
+
∂2u(k−1)(x)

∂x2

)
=
1

3

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)

(3.3)

Supposing the notations λ =1/∆t2 and µ = c/(2∆t), we obtain

(
λ+µ+

b

3

)
u(k+1)−

p

3

∂2u(k+1)(x)

∂x2
=
(
2λ−

b

3

)
u(k)+

p

3

∂2u(k)(x)

∂x2

+
(
−λ+µ−

b

3

)
u(k−1)+

p

3

∂2u(k−1)(x)

∂x2
+
1

3

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)
.

(3.4)

4. The meshless local weak form formulation

Instead of setting the global weak form, theMLRPImethod sets up theweak formover the local
quadrature cell such as Ωq, which is a small region taken for each node in the global domain Ω.
The local quadrature cells overlap with each other and cover the whole global domain Ω. The
local quadrature cells could be of any geometric shape and size. In one dimensional problems,
they are lines (intervals). The local weak form of Eq. (3.4) for xi ∈ Ω

i
q = (xi− rq,xi+ rq) can

be constructed as

∫

Ωiq

[
(
(
λ+µ+

b

3

)
u(k+1)−

p

3

∂2u(k+1)(x)

∂x2

]
ν(x) dx =

∫

Ωiq

[(
2λ−

b

3

)
u(k)+

p

3

∂2u(k)(x)

∂x2

]
ν(x) dx

+

∫

Ωiq

[
(
(
−λ+µ−

b

3

)
u(k−1)+

p

3

∂2u(k−1)(x)

∂x2

]
ν(x) dx

+

∫

Ωiq

[1
3

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)]
ν(x) dx

(4.1)



Meshless local radial point interpolation (MLRPI) for generalized telegraph... 575

where Ωiq is the local quadrature domain corresponding to the point i, and ν(x) is theHeaviside
step function defined by (Hu et al., 2006; Liu et al., 2006)

ν(x)=

{
1 x ∈ Ωq

0 x /∈ Ωq
(4.2)

as the test function in each local quadrature domain. Hence, we obtain

(
λ+µ+

b

3

)∫

Ωiq

u(k+1)ν(x) dx−
p

3

∫

Ωiq

∂2u(k+1)(x)

∂x2
ν(x) dx

=
(
2λ−

b

3

)∫

Ωiq

u(k)ν(x) dx+
p

3

∫

Ωiq

∂2u(k)(x)

∂x2
ν(x) dx

+
(
−λ+µ−

b

3

)∫

Ωiq

u(k−1)ν(x) dx+
p

3

∫

Ωiq

∂2u(k−1)(x)

∂x2
ν(x) dx

+
1

3

∫

Ωiq

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)
ν(x) dx

(4.3)

Using integration by parts, one obtains

∫

Ωiq

∂2u(k)(x)

∂x2
ν(x) dx = ν(x)

∂u(k)(x)

∂x

∣∣∣
x=xi+rq

x=xi−rq
−

∫

Ωiq

∂u(k)(x)

∂x

∂ν(x)

∂x
dx (4.4)

Then, by applying the test function, the following local weak equation is obtained

(
λ+µ+

b

3

)∫

Ωiq

u(k+1) dx−
p

3

(
∂u(k+1)(x)

∂x

∣∣∣
x=xi+rq

x=xi−rq

)

=
(
2λ−

b

3

)∫

Ωiq

u(k) dx+
p

3

(
∂u(k)(x)

∂x

∣∣∣
x=xi+rq

x=xi−rq

)

+
(
−λ+µ−

b

3

)∫

Ωiq

u(k−1) dx+
p

3

(
∂u(k−1)(x)

∂x

∣∣∣
x=xi+rq

x=xi−rq

)

+
1

3

∫

Ωiq

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)
dx

(4.5)

5. Discretization in the MLRPI method

In this Section, we consider Eq. (4.5) to see how to obtain discrete equations. Consider N
regularly located points on the boundary and domain of the problem, i.e. interval [0,1], so
that the distance between two consecutive nodes in each direction is constant and equal to h.
Assuming that u(xi,k∆t), i = 1,2, . . . ,N are known, our aim is to compute u(xi,(k +1)∆t),
i =1,2, . . . ,N. So,we have N unknowns and to compute these unknowns,we need N equations.
Toobtain thediscrete equations from locallyweak forms(4.5) for thenodes located in the interior
of the domain, i.e., for xi ∈ interior Ω, we substitute approximation formulas (2.8) into local
integral equations (4.5) to have



576 E. Shivanian, A. Khodayari

[(
λ+µ+

b

3

) N∑

j=1

(∫

Ωiq

φj(x) dx

)
−
p

3

N∑

j=1

(
∂φj(x)

∂x

∣∣∣
x=xi+rq

−
∂φj(x)

∂x

∣∣∣
x=xi−rq

)]
u
(k+1)
j

=

[(
2λ−

b

3

) N∑

j=1

(∫

Ωiq

φj(x) dx

)
+
p

3

N∑

j=1

(
∂φj(x)

∂x

∣∣∣
x=xi+rq

−
∂φj(x)

∂x

∣∣∣
x=xi−rq

)]
u
(k)
j

+

[(
−λ+µ−

b

3

) N∑

j=1

(∫

Ωiq

φj(x)dx

)
+
p

3

N∑

j=1

(
∂φj(x)

∂x

∣∣∣
x=xi+rq

−
∂φj(x)

∂x

∣∣∣
x=xi−rq

)]
u
(k−1)
j

+
1

3

∫

Ωiq

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)
dx

(5.1)

6. Numerical implementation of the MLRPI method

By using Simpson’s integration rule for nodes which are located on the boundary, we have for
all k

u(k)(x1)= γ1
h

3

[
u(k)(x1)+4u

(k)(x2)+2u
(k)(x3)+ . . .+4u

(k)(xN−1)+u
(k)(xN)

]
+µ1(k∆t)

u(k)(xN)= γ2
h

3

[
u(k)(x1)+4u

(k)(x2)+2u
(k)(x3)+ . . .+4u

(k)(xN−1)+u
(k)(xN)

]
+µ2(k∆t)

(6.1)

where x1 =0 and xN =1.
The matrix forms of Eqs. (5.1) and (6.1) for all N nodal points in the domain and the

boundary of the problem are given below
[(
λ+µ+

b

3

) N∑

j=1

Ai,j −
p

3

N∑

j=1

Bi,j

]
u
(k+1)
j =

[(
2λ−

b

3

) N∑

j=1

Ai,j +
p

3

N∑

j=1

Bi,j

]
u
(k)
j

+

[(
−λ+µ−

b

3

) N∑

j=1

Ai,j +
p

3

N∑

j=1

Bi,j

]
u
(k−1)
j +Ei(k−1,k,k +1)

(6.2)

where

Ai,j =

∫

Ωiq

φj(x) dx

Bi,j =

(
∂φj(x)

∂x

∣∣
x=xi+rq

−
∂φj(x)

∂x

∣∣∣
x=xi−rq

)

Ei(k−1,k,k +1)=
1

3

∫

Ωiq

(
f(k+1)(x)+f(k)(x)+f(k−1)(x)

)
dx

(6.3)

Assuming

Ai,j =
(
λ+µ+

b

3

)
Ai,j −

p

3
Bi,j Bi,j =

(
2λ−

b

3

)
Ai,j +

p

3
Bi,j

Ci,j =
(
−λ+µ−

b

3

)
Ai,j +

p

3
Bi,j U= {ui}N×1

Ek = [E1(k−1,k,k+1),E2(k−1,k,k +1), . . . ,EN(k−1,k,k +1)]
T

(6.4)

yields



Meshless local radial point interpolation (MLRPI) for generalized telegraph... 577

AU(k+1) =BU(k)+CU(k−1)+Ek (6.5)

Furthermore, to satisfy Eqs. (6.1), for both nodes belong to the boundary, i.e., {x1,xN}, we set

Eki =

{
µ1(k∆t) i =1

µ2(k∆t) i = N

∀j : Bi,j =Ci,j =0 i =1,N

A1 =
[
1−γ1

h

3
,−4γ1

h

3
,−2γ1

h

3
, . . . ,−4γ1

h

3
,−γ1

h

3

]

AN =
[
−γ2

h

3
,−4γ2

h

3
,−2γ2

h

3
, . . . ,−4γ2

h

3
,1−γ2

h

3

]

(6.6)

whereA1 andAN are the first and N-th rows of the matrixA, respectively.
At the first time level, when n = 0, according to the initial conditions that are introduced

in Eq. (1.2), we apply the following assumptions

u(0) = u0 u
(−1) ∼= u

(1)−2∆tψ(x)

where

u0 = [u0(x1),u0(x2), . . . ,u0(xN)]
T ψ = [ψ(x1),ψ(x2), . . . ,ψ(xN)]

T

7. Numerical experiments

In this Section, two numerical expriments for application of the meshless local radial point
interpolation method (MLRPI) in solving the one-dimensional linear telegraph equation with
non-local boundary conditions are presented. In both examples, the domain integrals are evalu-
ated with 3 points Gaussian quadrature rule. In these problems, the regular distributed nodal
points are used. Also, in order to implement the meshless local weak form in these cases, the
radius of the local quadrature domain rq =0.8h is selected, where h is the distance between the
nodes in the x direction (h = ∆x). The size of rq is such that the union of these sub-domains
must cover the whole global domain. The radius of the support domain to the local radial point
interpolation method is rs = 4rq. This size is significant enough to have a sufficient number of
nodes (n) to give appropriate shape functions. Also, in Eq. (2.1), we set m =5.

Example 1. We set c =20, b =25 and p =1. The exact solution of the first example is taken
as u(x,t) = t3(2x3 −x +4), (x,t) ∈ [0,1]× [0,1]. According to this exact solution, f(x,t) is
given by

f(x,t)= (6t+3ct2)(2x3−x+4)+ bt3(2x3−x+4)−12pxt3

the initial conditions are

u(x,0)= 0
∂u

∂t
(x,0)= 0

and the non-local boundary conditions take the form

u(0, t) =

1∫

0

u(x,t) dx t ­ 0

u(1, t) =
5

4

1∫

0

u(x,t) dx t ­ 0



578 E. Shivanian, A. Khodayari

Tables 1 and 2 as well as Fig. 1a show the results of the MLRPI method to solve Example 1
using TPS as the radial basis function. Also, Tables 3 and 4 as well as Fig. 1b illustrate the
results of the current method to solve Example 1 using MQ as the radial basis function. As
it is seen, the MLRPI method is of high accuracy. Furthermore, it is seen that the method is
convergent with respect to the spatial and time variable using both TPS andMQ.

Fig. 1. Numerical solutions and the exact solution at time t =1.0 for Example 1: (a) using TPS,
(b) using MQ. The solid line corresponds to the exact solution, the starred line corresponds to the

numerical solution of theMLRPI with ∆t =0.0001 and ∆x =0.0125

Table 1. The L1, L2 and L∞ errors calculated by MLRPI for Example 1 with different ∆x
and ∆t at time t =1.0 (using TPS)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞

0.001 0.25 4.070835E-03 2.136594E-03 1.508984E-03

0.001 0.125 3.470766E-04 1.792837E-04 1.587210E-04

0.001 0.1 1.718355E-04 7.995307E-05 6.922815E-05

0.001 0.05 6.156571E-05 1.401474E-05 5.949385E-06

Table 2.The L1, L2 and L∞ errors calculated byMLRPI for Example 1 with different ∆x and
∆t at time t =1.0 (using TPS)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞

0.001 0.0125 1.945099E-04 2.169317E-05 2.997031E-06

0.0005 0.0125 4.863614E-05 5.424111E-06 7.494035E-07

0.00025 0.0125 1.217382E-05 1.357562E-06 1.875924E-07

0.0001 0.0125 1.964876E-06 2.192364E-07 3.026972E-08

Table 3.The L1, L2 and L∞ errors calculated byMLRPI for Example 1 with different ∆x and
∆t at time t =1.0 (usingMQ)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞

0.001 0.25 4.070835E-03 2.136594E-03 1.508984E-03

0.001 0.125 3.459171E-04 1.780828E-04 1.575668E-04

0.001 0.1 1.715535E-04 7.972770E-05 6.898580E-05

0.001 0.05 6.154538E-05 1.401788E-05 5.974310E-06



Meshless local radial point interpolation (MLRPI) for generalized telegraph... 579

Table 4.The L1, L2 and L∞ errors calculated byMLRPI for Example 1 with different ∆x and
∆t at time t =1.0 (usingMQ)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞

0.001 0.0125 1.942840E-04 2.166659E-05 2.993550E-06

0.0005 0.0125 4.841030E-05 5.397581E-06 7.459259E-07

0.00025 0.0125 1.194834E-05 1.331264E-06 1.841192E-07

0.0001 0.0125 1.740619E-06 1.944979E-07 2.682201E-08

Example 2. We set c =20, b =25 and p =1. The exact solution of the this example is taken
as u(x,t) = t3 sin(x+1), (x,t)∈ [0,1]× [0,1]. According to this exact solution, f(x,t) is given
by

f(x,t)= (6t+3ct2)sin(x+1)+ bt3 sin(x+1)+pt3 sin(x+1)

the initial conditions are

u(x,0)= 0
∂u

∂t
(x,0)= 0

and the non-local boundary conditions take the form

u(0, t) = 0.8797864387

1∫

0

u(x,t) dx t ­ 0

u(1, t) = 0.950701283

1∫

0

u(x,t) dx t ­ 0

Tables 5 and 6 as well as Fig. 2a show the results of the MLRPI method to solve Example 2
using TPS as the radial basis function. Besides, Tables 7 and 8 as well as Fig. 2b demonstrate
the results of the presentmethod to solve Example 2 usingMQ as the radial basis function. As
it is seen, the MLRPI method is of high accuracy. Moreover, we see that the convergence with
respect to both the time step (∆t) and the number of nodal points (N) are hold, no matter
which kind of RBFwe use.

Fig. 2. Numerical solutions and the exact solution at time t =1.0 for Example 2: (a) using TPS,
(b) using MQ. The solid line corresponds to the exact solution, the starred line corresponds to the

numerical solution of theMLRPI with ∆t =0.0001 and ∆x =0.0125

On the top of that, theMLRPImethod can be used to solve complex engineering problems
with lower computational cost, higher accuracy, simpler construction of higher-order shape func-
tions and easier handling of large deformation problems.



580 E. Shivanian, A. Khodayari

Table 5.The L1, L2 and L∞ errors calculated byMLRPI for Example 2 with different ∆x and
∆t at time t =1.0 (using TPS)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞

0.001 0.25 4.467206E-04 2.263494E-04 1.434185E-04

0.001 0.125 2.548359E-05 1.304920E-05 1.148951E-05

0.001 0.1 1.502825E-05 6.430589E-06 5.320887E-06

0.001 0.05 1.175112E-05 2.629557E-06 8.956239E-07

Table 6.The L1, L2 and L∞ errors calculated byMLRPI for Example 2 with different ∆x and
∆t at time t =1.0 (using TPS)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞
0.001 0.0125 4.418773E-05 4.914641E-06 5.709157E-07

0.0005 0.0125 1.104494E-05 1.228475E-06 1.426643E-07

0.00025 0.0125 2.7603974E-06 3.070796E-07 3.561701E-08

0.0001 0.0125 4.405399E-07 4.920947E-08 6.985131E-09

Table 7.The L1, L2 and L∞ errors calculated byMLRPI for Example 2 with different ∆x and
∆t at time t =1.0 (usingMQ)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞
0.001 0.25 4.467206E-04 2.263494E-04 1.434185E-04

0.001 0.125 2.542164E-05 1.296910E-05 1.140565E-05

0.001 0.1 1.500617E-05 6.414381E-06 5.302278E-06

0.001 0.05 1.174746E-05 2.629195E-06 8.987330E-07

Table 8.The L1, L2 and L∞ errors calculated byMLRPI for Example 2 with different ∆x and
∆t at time t =1.0 (usingMQ)

∆t ∆x ‖E‖1 ‖E‖2 ‖E‖∞
0.001 0.0125 4.413908E-05 4.909225E-06 5.701902E-07

0.0005 0.0125 1.099625E-05 1.223055E-06 1.419426E-07

0.00025 0.0125 2.711815E-06 3.016704E-07 3.548134E-08

0.0001 0.0125 3.928077E-07 4.390304E-08 5.957403E-09

8. Conclusions

In the aforementioned discussion, we applied the meshless local radial point interpolation
(MLRPI) method to solve the linear telegraph equation with non-local boundary conditions.
The radial point interpolation method is adopted for approximating the field variable. Also the
weak form of the discretized equations has been constructed on local subdomains. So, this me-
thod requires neither domain element nor background cells in either the interpolation or the
intergration. Itmeans thismethod is a trulymeshlessmethod. Furthermore, time discretization
hasbeendoneusingfinitedifference techniques.Theprincipal benefit of themethod is to capture
the behavior of the solution for similar problemswith non-local boundary conditionswheremost
of schemes fail. Also, theMLRPImethod can easily handle the damage of the components, such
as fracture which is very useful to simulate material breakage. Finally, accuracy and usefulness
of the proposedmethod are illustrated by two examples.



Meshless local radial point interpolation (MLRPI) for generalized telegraph... 581

Acknowledgments

The authors are grateful to anonymous reviewers for carefully reading this paper and for their com-

ments and suggestions which have improved the paper.

References

1. Almenar P., Jodar L., Martin J.A., 1997,Mixed problems for the time-dependent telegraph
equation: Continuous numerical solutions with a priori error bounds,Mathematical and Computer
Modelling, 25, 11, 31-44

2. Atluri S., Zhu T., 1998, A new meshless local Petrov-Galerkin (MLPG) approach in computa-
tional mechanics,Computational Mechanics, 22, 117-127

3. Banasiak J.,Mika J.R., 1998, Singularly perturbed telegraph equationswith applications in the
randomwalk theory, Journal of Applied Mathematics and Stochastic Analysis, 11, 1, 9-28

4. Belytschko T., Lu Y.Y., Gu L., 1995, Element free Galerkin methods for static and dynamic
fracture, International Journal of Solids and Structures, 32, 2547-2570

5. BoyceW.E., DiPrimaR.C., 1977,Differential Equations Elementary and Boundary Value Pro-
blems, Wiley, NewYork

6. Bratsos A.G., 2008,An improved numerical scheme for the sine-Gordon equation in 2+1 dimen-
sions, International Journal for Numerical Methods in Engineering, 75, 787-799

7. Ciment M., Leventhal S.H., 1978, A note on the operator compact implicit method for the
wave equation,Mathematics of Computation, 32, 143-147

8. Dashtimanesh A., Ghadimi P., 2013,A three-dimensional SPHmodel for detailed study of free
surface deformation, just behind a rectangular planing hull, Journal of the Brazilian Society of
Mechanical Sciences and Engineering, 35, 4, 369-380

9. De S., Bathe K.J., 2000, Themethod of finite spheres,Computational Mechanics, 25, 329-345

10. Dehghan M., Mirzaei D., 2008, The meshless local Petrov-Galerkin (MLPG) method for the
generalized two-dimensional non-linear Schrödinger equation,Engineering Analysis with Boundary
Elements, 32, 747-756

11. Dehghan M., Shokri A., 2008, A numerical method for solution of the two dimensional sine-
Gordon equation using the radial basis functions,Mathematics and Computers in Simulation, 79,
700-715

12. Gonzalez-VelascoE.A., 1995,FourierAnalysis andBoundaryValueProblems,AcademicPress,
NewYork

13. Gu Y., Liu G., 2002, A boundary point interpolation method for stress analysis of solids, Com-
putational Mechanics, 28, 47-54

14. GuY.T., LiuG.R., 2003,A boundary radial point interpolationmethod (BRPIM) for 2-D struc-
tural analyses, Structural Engineering and Mechanics, 15, 535-550

15. Hosseini V.R., Shivanian E., Chen W., 2015, Local integration of 2-D fractional telegraph
equation via local radial point interpolant approximation, European Physical Journal Plus, 130,
33-54

16. Hu D., Long S., Liu K., Li G., 2006, A modified meshless local Petrov-Galerkin method to
elasticity problems in computer modeling and simulation, Engineering Analysis with Boundary
Elements, 30, 399-404

17. Jordan P.M., Puri A., 1999, Digital signal propagation in dispersivemedia, Journal of Applied
Physics, 85, 3, 1273-1282

18. Kansa E., 1990,Multiquadrics-a scattered data approximation scheme with applications to com-
putational fluid-dynamics. I. Surface approximations and partial derivative estimates, Computers
and Mathematics with Applications, 19, 8/9, 127-145



582 E. Shivanian, A. Khodayari

19. Kochmann D.M., Venturini G.N., 2014, A meshless quasicontinuum method based on local
maximum-entropy interpolation,Modelling and Simulation in Materials Science and Engineering,
22, 3, 034007

20. Libre N., Emdadi A., Kansa E., ShekarchiM., Rahimian M., 2008,A fast adaptivewavelet
scheme in RBF collocation for nearly singular potential PDEs,Computer Modeling in Engineering
and Sciences, 38, 3, 263-284

21. Ling L., Schaback R., 2008, Stable and convergent unsymmetric meshless collocationmethods,
SIAM Journal of Numerical Analysis, 46, 3, 1097-1115

22. Liu G.R., Gu Y.T., 2001, A local radial point interpolationmethod (LR-PIM) for free vibration
analyses of 2-D solids, Journal of Sound and Vibration, 246, 1, 29-46

23. Liu G., Gu Y., 2005,An Introduction to Meshfree Methods and their Programing, Springer

24. Liu K., Long S., Li G., 2006, A simple and less-costly meshless local Petrov-Galerkin (MLPG)
method for thedynamic fractureproblem,EngineeringAnalysis withBoundaryElements,30, 72-76

25. LiuW.K., JunS.,ZhangY.F., 1995,Reproducingkernelparticlemethods, International Journal
for Numerical Methods in Engineering, 20, 1081-1106

26. Melenk J.M., Babuska I., 1996, The partition of unity finite elementmethod: Basic theory and
applications,Computer Methods in Applied Mechanics and Engineerin, 139, 289-314

27. Mukherjee Y.X., Mukherjee S., 1997, Boundary node method for potential problems, Inter-
national Journal for Numerical Methods in Engineering, 40, 797-815

28. Nayroles B., Touzot G., Villon P., 1992, Generalizing the finite element method: diffuse
approximation and diffuse elements,Computational Mechanics, 10, 307-318

29. Pan X., Yuan H., 2009, Applications of meshless methods for damage computations with finite
strains,Modelling and Simulation in Materials Science and Engineering, 17, 4, 045005

30. Shirzadi A., 2014, Solving 2D reaction-diffusion equations with nonlocal boundary conditions by
the RBF-MLPGmethod, Computational Mathematics and Modeling, 25, 4, 521-529

31. Shivanian E., 2013, Analysis of meshless local radial point interpolation (MLRPI) on a nonline-
ar partial integro-differential equation arising in population dynamics, Engineering Analysis with
Boundary Elements, 37, 1693-1702

32. ShivanianE., 2015a,Anew spectralmeshless radial point interpolation (SMRPI)method:Awell-
behaved alternative to the meshless weak forms, Engineering Analysis with Boundary Elements,
54, 1-12

33. Shivanian E., 2015b,Meshless local Petrov-Galerkin (MLPG)method for three-dimensional non-
linearwave equations viamoving least squares approximation,EngineeringAnalysis withBoundary
Elements, 50, 249-257

34. Shivanian E., Khodabandehlo H.R., 2014,Meshless local radial point interpolation (MLRPI)
on the telegraph equation with purely integral conditions, European Physical Journal Plus, 129,
241-251

35. Singh I.V.,TanakaM.,EndoM., 2007,Meshlessmethod for nonlinearheat conduction analysis
of nano-composites,Heat and Mass Transfer, 43, 10, 1097-1106

36. Sladek J., Sladek V., Zhang C., Tan C.L., 2006, Evaluation of fracture parameters for crack
problems in FGM by a meshless method, Journal of Theoretical and Applied Mechanics, 44, 3,
603-636

37. Tikhonov A.N., Samarskii A.A., 1990,Equations of Mathematical Physics, Dover, NewYork

Manuscript received May 24, 2016; accepted for print December 5, 2016