Computers and Mathematics with Applications 56 (2008) 2175–2188

Contents lists available at ScienceDirect

Computers and Mathematics with Applications

journal homepage: www.elsevier.com/locate/camwa

The method of lines for solution of the one-dimensional wave equation

subject to an integral conservation condition
Fatemeh Shakeri, Mehdi Dehghan ∗
Department of Applied Mathematics, Faculty of Mathematics and Computer Sciences, Amirkabir University of Technology, No. 424, Hafez Avenue,
Tehran 15914, Iran

article info a b s t r a c t

Article history: Hyperbolic partial differential equations with an integral condition serve as models in
Received 14 November 2007 many branches of physics and technology. Recently, much attention has been expended
Received in revised form 22 March 2008 in studying these equations and there has been a considerable mathematical interest in
Accepted 25 March 2008
them. In this work, the solution of the one-dimensional nonlocal hyperbolic equation is
presented by the method of lines. The method of lines (MOL) is a general way of viewing
Keywords:
a partial differential equation as a system of ordinary differential equations. The partial
Hyperbolic partial differential equation
Non-classic boundary condition
derivatives with respect to the space variables are discretized to obtain a system of ODEs
Method of lines in the time variable and then a proper initial value software can be used to solve this ODE
System of ordinary differential equations system. We propose two forms of MOL for solving the described problem. Several numerical
Integral conservation condition examples and also some comparisons with finite difference methods will be investigated
to confirm the efficiency of this procedure.
© 2008 Elsevier Ltd. All rights reserved.

1. Introduction

Over the last few years, it has become increasingly apparent that many physical phenomena can be described in terms of
hyperbolic partial differential equations with an integral condition replacing the classic boundary condition [1]. This type of
equations arises, for example in the study of thermoelasticity [2,3], plasma physics [4], chemical heterogeneity [5,6] and etc.
Growing attention is being paid to the development, analysis and implementation of numerical methods for the solution of
these problems. Hyperbolic initial boundary value problems in one dimension that involve nonlocal boundary conditions
have been studied by several authors [7,8,1,9–12]. For parabolic equations subject to nonlocal boundary conditions the
interested reader can see [13] and the references therein. Also nonlocal problems include the problems with an integral
term in initial condition [14]. In the current work we will not discuss on this group.
In this research, we consider the following problem of this family of equations

∂ 2v ∂ 2v
− 2 = q(x, t ), 0 ≤ x ≤ l, 0 < t ≤ T , (1.1)
∂t 2 ∂x
with initial conditions
v(x, 0) = f1 (x), 0 ≤ x ≤ l, (1.2)
and
vt (x, 0) = f2 (x), 0 ≤ x ≤ l, (1.3)

∗ Corresponding author.
E-mail address: [email protected] (M. Dehghan).

0898-1221/$ – see front matter © 2008 Elsevier Ltd. All rights reserved.
doi:10.1016/j.camwa.2008.03.055
2176 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

and Dirichlet boundary condition

v(0, t ) = g1 (t ), 0 < t ≤ T, (1.4)
and the nonlocal condition
Z l
v(x, t )dx = g2 (t ), 0 < t ≤ T, (1.5)
0

where q, f1 , f2 , g1 and g2 are known functions. We assume that the functions q, f1 , f2 , g1 , and g2 satisfy the conditions in
order that the solution of this equation exists and is unique. The existence and uniqueness of the solution of this problem
are discussed in [15].
Author of [1] presented several finite difference schemes for the numerical solution of problem (1.1)–(1.5). These three-
level techniques are based on two second-order (one explicit and one weighted) schemes and a fourth-order technique (a
weighted explicit) [1]. Also in [16] the shifted Legendre Tau technique is developed for the solution of the studied model.
The approach in their work consists of reducing the problem to a set of algebraic equations by expanding the approximate
solution as a shifted Legendre function with unknown coefficients. The integral and derivative operational matrices are
given. These matrices together with the tau method are then utilized to evaluate the unknown coefficients of shifted
Legendre functions. Author of [17] developed a numerical technique based on an integro-differential equation and local
interpolating functions for solving the one-dimensional wave equation subject to a nonlocal conservation condition and
suitably prescribed initial boundary conditions. Authors of [18] combined finite difference and spectral methods to solve
the one-dimensional wave equation with an integral condition. The time variable is approximated using a finite difference
scheme. But the spectral method is employed for discretizing the space variable. The main idea behind their approach is
that high-order results can be obtained. They have tested the new method for two examples from the literature [18]. A
numerical scheme for solving the second-order wave equation with given initial conditions and a boundary condition and
an integral condition in place of the classical boundary condition is presented in [19]. The cubic B-spline scaling function
together with boundary scaling function on interval [0, 1] are employed for solving the model. The obtained results show
that this approach can solve the problem effectively. It is worth pointing out that the matrices are not sparse such as those
that can be obtained when one uses the finite difference methods, but here [19] we need smaller matrices to get some results
to be compared with finite difference methods.
In this work a different approach is used, the solution of the above equation is computed by the method of lines. Method
of lines is an alternative computational approach which involves making an approximation to the space derivatives and
by reducing the problem to that of solving a system of initial value ordinary differential equations and then using a time
integrator for solving the ODE system. One can increase the accuracy of the method by the use of highly efficient reliable
initial value ODE solvers which means that comparable orders of accuracy can also be achieved in the time integration
without using extremely small time steps.
This paper is organized in the following way. In Section 2, we introduce the method of lines briefly and apply it to
Eqs. (1.1)–(1.5) in two different ways. Some numerical results and comparisons with the finite difference method presented
in [1] are given in Section 3 and finally a conclusion is drawn in Section 4.

2. Method of lines

Method of lines is a semi-discrete method [20–25] which involves reducing an initial boundary value problem to a system
of ordinary differential equations (ODE) in time through the use of a discretization in space. The resulting ODE system
can be solved using the standard initial value software, which may use a variable time-step/variable order approach with
time local error control. The most important advantage of the MOL approach is that it has not only the simplicity of the
explicit methods [26] but also the superiority (stability advantage) of the implicit ones unless a poor numerical method for
solution of ODEs is employed. It is possible to achieve higher-order approximations in the discretization of spatial derivatives
without significant increases in the computational complexity. This technique has the broad applicability to physical and
chemical systems modeled by PDEs. The models that include the solution of mixed systems of algebraic equations, ODEs
and PDEs, the resolution of steep moving fronts, parameter estimation and optimal control, other problems such as delay
differential equations [27], two-dimensional sine-Gordon equation [28], the Nwogu one-dimensional extended Boussinesq
equation [29], partial differential equation problems describing nonlinear wave phenomena, e.g., a fully nonlinear third-
order Korteweg–de Vries (KdV) equation, the fourth-order Boussinesq equation, the fifth-order Kaup–Kupershmidt equation
and an extended KdV5 equation [30], nonlinear inverse heat conduction problem [31], interface problem [32], multi-
component atmospheric pollutant propagation model with pollutants phase transformation consideration [33], a Bingham
problem in cylindrical pipes [34], a mathematical model for capillary formation [35], three-dimensional transient radiative
transfer equation [36], elliptic partial differential equations which describe steady-state mass and energy transport in
solids [37] and many other physical problems.
Author of [38] used MOL to transform the initial boundary value problem associated with the nonlinear hyperbolic
Boussinesq equation, into a first-order nonlinear initial value problem. Numerical methods are developed by replacing
the matrix-exponential term into a recurrence relation by rational approximants. MOL is employed in [25] to solve the
 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188 2177

Korteweg–de Vries equation. Authors of [39] investigated the performance of various terms of upwinding to provide some
guidance in the selection of upwind methods in the MOL solution of strongly convective partial differential equations. MOL
is used in [40] to obtain numerical solution of a mathematical model for capillary formation in tumor angiogenesis. Authors
of [41] applied moving and adaptive mesh methods for the numerical solution of an extended fifth-order Korteweg–de
Vries model motivated by water waves in the presence of surface tension. Also they investigated the dynamics and
interaction of embedded solitons. Authors of [42] applied MOL for solving parabolic partial differential equations. They
showed (experimentally) that increasing the points in the formula in the spatial approximation cannot bring the MOL and
exact solutions into close agreement(for parabolic problems but not necessarily for hyperbolic or elliptic equations). The
method of lines and finite difference method were tested in [43] from the viewpoints of solution accuracy and central
processing unit time by applying them to the solution of time-dependent two-dimensional Navier–Stokes equations for
transient laminar flow without/with sudden expansion and comparing their results with steady-state numerical predictions
and measurements previously reported in the literature. MOL was found to be superior to finite difference method with
respect to CPU and set-up times and its flexibility for incorporation of other conservative equations [43]. The method of
lines concept has been combined with the boundary element based elimination of the spatial derivatives to obtain a solution
method for partial differential equations which are parabolic in time [44]. The boundary element method alleviates the
need for spatial discretization and casts the problem in an integral format. Hence errors associated with the numerical
approximation of the spatial derivatives are totally eliminated. Applications of the method of lines to the hyperbolic partial
differential equations with stiff nonlinear source terms is investigated in [45]. Some of the options available for time
integration when using a moving grid method of lines code is surveyed in [46]. A new technique is proposed in [47] for
the numerical integration of the system of ordinary differential equations that arises in the method of lines solution of time-
dependent partial differential equations. This system is usually stiff, so it is desirable for the numerical method to solve it to
have good properties concerning stability.
In order to use this approach for solving (1.1)–(1.5), we discretize the coordinate x with M (M even) uniformly spaced
grid points xi = xi−1 + h, x0 = 0, xM = l, i = 1, 2, . . . , M. Note that h = l/M and we can also write xi = ih. At first, in part (I),
we use a second-order difference approximation for the second derivative in x in grid points xi , i = 1, 2, . . . , M − 1 and in
part (II), we apply a fourth-order difference approximation to the second derivative in x in grid points xi , i = 2, 3, . . . , M − 2.

2.1. MOL I

Let vi (t ) approximate v(xi , t ). In Fig. 1, the lines along which the approximations vi (t ) are defined, are shown. Using the
second-order central difference approximation for the second derivative in x results in
d2 v i vi−1 − 2vi + vi+1
− = q(xi , t ), i = 1, . . . , M − 1, 0 < t ≤ T . (2.1)
dt 2 h2
Also the conditions (1.2)–(1.4) are reduced to
vi (0) = f1 (xi ), 0 ≤ x ≤ l, (2.2)
(vi )t (0) = f2 (xi ), 0 ≤ x ≤ l, (2.3)
and
v0 (t ) = g1 (t ), 0 < t ≤ T. (2.4)
Applying Simpson’s numerical integration rule, we approximate the nonlocal condition (1.5) as follows
Z l M
X
v(x, t )dx = h ci vi (t ) = g2 (t ), 0 < t ≤ T, (2.5)
0 i=0

where c0 = cM = 31 , c2k−1 = 43 for k = 1, . . . , M

2
and c2k = 2
3
, for k = 1, . . . , M
2
− 1. Now, using (2.4) and (2.5), we can
obtain the following formula for vM (t )
" # " #
M −1 M −1
1 X 1 X
vM (t ) = − h ci vi (t ) − g2 (t ) =− hc0 g1 (t ) + h ci vi (t ) − g2 (t ) . (2.6)
hcM i=0
hcM i=1

Before we solve the system of ODE (2.1) with initial conditions (2.2) and (2.3), we apply some transformations to convert
this system to a system of the first-order equations
dv1
u1 (t ) = ,
dt
dv2
u2 (t ) = ,
dt
..
.
dvM −1
uM −1 (t ) = ,
dt
2178 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

Fig. 1. The colored strip is the domain on which the PDE is defined. The approximations vi (t ) are defined along the dashed lines.

Fig. 2. Plot of the exact solution in Example 1.

which results in the following system of 2M − 2 equations

dui vi−1 − 2vi + vi+1
− = q(xi , t ), i = 1, . . . , M − 2, 0 < t ≤ T ,
dt h2
M −3
cM −2 cM −1
vM −2 (1 − ) − vM −1 (2 + )− 1
ci vi (t )
P
cM cM cM
duM −1 i=1 1
− = q(xM −1 , t ) − [hc0 g1 (t ) − g2 (t )],
dt h2 h3 cM
dvi
= ui (t ), i = 1, . . . , M − 1, 0 < t ≤ T ,
dt
with initial conditions
vi (0) = f1 (xi ), ui (0) = f2 (xi ).

2.2. MOL II

In this part, for the second derivative in x in grid points xi , i = 2, 3, . . . , M − 2, we use the following fourth-order
difference approximation
vi+2 − 16vi+1 + 30vi − 16vi−1 + vi−2
(vi )xx = − , i = 2, . . . , M − 3,
12h2
cM −1 cM −2
1
hcM
[hc0 g1 (t ) − g2 (t )] + vM −1 (16 + cM
) − vM −2 (30 + cM
)
(vM −2 )xx =
12h2
M −5
cM −3 cM −4
vM −3 (16 + ) − vM −4 (1 + )+ 1
ci vi (t )
P
cM cM cM
i=1
+ .
12h2
Using transformations introduced in previous part, an ODE system is obtained as follows
du1 v0 − 2v1 + v2
− = q(x1 , t ), 0 < t ≤ T ,
dt h2
dui vi+2 − 16vi+1 + 30vi − 16vi−1 + vi−2
+ = q(xi , t ), i = 2, . . . , M − 3, 0 < t ≤ T ,
dt 12h2
 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188 2179

Fig. 3. Plot of the vi (t ) for i = 5k, k = 1, . . . , 10 in Example 1, the graphs of vi (t ) for i = 5k, k = 11, . . . , 19, are symmetries of vi (t ), i = 5k, k = 9, . . . , 1,
regarding the horizontal axis respectively.

Fig. 4. Plot of the v(xi , t ) − vi (t ) for (a1) xi = 1

4
, (a2) xi = 1
2
, (a3) xi = 3
4
, (a4) xi = 1 in Example 1.

M −5
cM −3 cM −4
vM −3 (16 − ) − vM −4 (1 + )+ 1
ci vi (t )
P
cM −1 cM −2
duM −2 vM −1 (16 − cM
) − vM −2 (30 + cM
) cM cM cM
i=1
− −
dt 12h2 12h2
1
= q(xM −2 , t ) + [hc0 g1 (t ) − g2 (t )], 0 < t ≤ T,
12h3 cM
2180 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

Table 1
Computational results for Example 1
xi Exact value Absolute error by [1] Absolute error by MOL I Absolute error by MOL II

0.10 0.95105651629515 3.3 × 10−5 1.5970339260263 × 10−5 9.1695313453322 × 10−8

0.20 0.80901699437495 3.0 × 10−5 2.2782454249470 × 10−5 8.0261219004285 × 10−8
0.30 0.58778525229247 3.2 × 10−5 2.0901699376741 × 10−5 5.9441112476577 × 10−8
0.40 0.30901699437495 3.1 × 10−5 1.2281581424911 × 10−5 3.3507949526168 × 10−8
0.50 0.00000000000000 3.3 × 10−5 1.612449505348880 × 10−17 1.087626749329525 × 10−16
0.60 −0.30901699437495 3.4 × 10−5 1.2281581424800 × 10−5 3.3507949692702 × 10−8
0.70 −0.58778525229247 3.1 × 10−5 2.0901699376519 × 10−5 5.9441111255332 × 10−8
0.80 −0.80901699437495 3.2 × 10−5 2.2782454249470 × 10−5 8.0261216561794 × 10−8
0.90 −0.95105651629515 3.4 × 10−5 1.5970339258931 × 10−5 9.1695316228879 × 10−8
1.00 −1.00000000000000 3.2 × 10−5 8.437694987151190 × 10−15 2.664535259100376 × 10−15

M −3
cM −2 cM −1
vM −2 (1 − ) − vM −1 (2 + )− 1
ci vi (t )
P
cM cM cM
duM −1 i=1 1
− = q(xM −1 , t ) − [hc0 g1 (t ) − g2 (t )],
dt h2 h3 cM
0 < t ≤ T,
dvi
= ui (t ), i = 1, . . . , M − 1, 0 < t ≤ T .
dt
Now we can solve the resulting ODE system using an ODE solver with time local error control. As is mentioned in [23]
the basic idea of the method of lines is to replace the spatial derivatives in the partial differential equation with algebraic
approximations. Once this is carried out, the spatial derivatives are no longer stated explicitly in terms of the spatial
independent variables. Thus, in effect only the initial value variable, typically time in a physical problem, remains. In other
words, with only one remaining independent variable, we have a system of ordinary differential equations that approximate
the original PDE. The challenge, then, is to formulate the approximating system of ODEs. Once this is carried out, we can apply
any integration algorithm for initial value ODEs to compute an approximate numerical solution to the PDE. Thus, one of the
salient features of the MOL is the use of existing, and generally well-established, numerical methods for ODEs [22]. In the
current work, to solve the obtained ODE system, we use ode45 solver in MATLAB which is based on the explicit Runge–Kutta
(4, 5) formula, the Dormand–Prince pair.

3. Numerical examples

In this section, some examples will be investigated to show the reliability and efficiency of the proposed schemes in this
paper.

3.1. Test 1

As the first example, consider the Eqs. (1.1)–(1.5) with l = 1, T = 4,

q(x, t ) = 0,
f1 (x) = 0, f2 (x) = π cos(π x),
g1 (t ) = sin(π t ), g2 (t ) = 0.
The exact solution of this equation is
v(x, t ) = cos(π x) sin(π t ).
We compare the results obtained by the procedure in previous section with finite difference method introduced in [1] in
Table 1. For comparison purpose, we set h = 0.01 and calculate the absolute errors for t = 21 (the final time T in [1] is 12 ).
Also the maximum norms of errors which are considered in the following are reported in Tables 2 and 3.
kvexact − vMOL kxi ,∞ = max |v(xi , t ) − vi (t )|.
0<t ≤T

kvexact − vMOL k∞ = max max |v(xi , t ) − vi (t )|.

1≤i≤M 0<t ≤T

Some other numerical results are shown in Figs. 2–4. It is worth noting that since discretization in MOL is only applied to
the spatial variable, increasing the final time T does not decrease the accuracy of the solution and we can obtain an improved
solution by using a more accurate integrator for solving the ODE system, but in finite difference methods, there is usually
limitations when final time T increases.
As it may be seen from Table 1, the results obtained using MOL I, have generally the same accuracy as those obtained
from finite difference method in [1], but MOL II is more accurate than these two schemes.
 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188 2181

Table 2
The norm kvexact − vMOL kxi ,∞ in Example 1

xi MOL I MOL II

0.10 1.9404299881232 × 10 −5
1.0078024292870 × 10−7
0.20 2.9082139621384 × 10−5 1.0395382088468 × 10−7
0.30 2.7979999464800 × 10−5 9.519132848634 × 10−8
0.40 1.6949536840560 × 10−5 6.254244772075 × 10−8
0.50 1.109874786267721 × 10−14 2.177963786387436 × 10−14
0.60 1.6949536832678 × 10−5 6.254244200310 × 10−8
0.70 2.7979999469574 × 10−5 9.519131310975 × 10−8
0.80 2.9082139622272 × 10−5 1.0395382143980 × 10−7
0.90 1.9404299884340 × 10−5 1.0078024492710 × 10−7
1.00 5.933031843596837 × 10−13 2.680078381445128 × 10−13

Fig. 5. Plot of the exact solution in Example 2.

Fig. 6. Plot of the vi (t ) for i = 10k, k = 1, . . . , 5 in Example 2, the graphs of vi (t ) for i = 10k, k = 6, . . . , 9, are symmetries of vi (t ), i = 10k, k = 4, . . . , 1,
regarding the horizontal axis respectively.

3.2. Test 2

Consider the problem (1.1)–(1.5) with l = 1, T = 5,

q(x, t ) = 0,
f1 (x) = cos(π x), f2 (x) = 0,
g1 (t ) = cos(π t ), g2 (t ) = 0.

v(x, t ) = (cos(π (x + t )) + cos(π (x − t ))) is the exact solution of this system. The computational results obtained by MOL
1
2
I and MOL II with h = 0.01 are compared with those investigated in [1] in Table 4. Because of the comparison, the absolute
errors are tabulated at t = 41 . From this table, we see that the solution based on the MOL II is more accurate as compared
to MOL I or the finite difference method developed in [1]. The norm of errors resulted by MOL I and MOL II, graphs of the
2182 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

Fig. 7. Plot of the v(xi , t ) − vi (t ) for (a1) xi = 1

4
, (a2) xi = 1
2
, (a3) xi = 3
4
, (a4) xi = 1 in Example 2.

Table 3
The norm kvexact − vMOL k∞ in Example 1
MOL I MOL II

2.9082139622272 × 10−5 1.0395382143980 × 10−7

Table 4
Computational results for Example 2
xi Exact value Absolute error by [1] Absolute error by MOL I Absolute error by MOL II

0.10 0.67249851196396 5.2 × 10−5 1.2923596914405 × 10−5 8.6007681532330 × 10−8

0.20 0.57206140281768 5.1 × 10−5 1.7467564661144 × 10−5 8.9389747826019 × 10−8
0.30 0.41562693777745 5.1 × 10−5 1.3423564604209 × 10−5 1.731384369208 × 10−9
0.40 0.21850801222441 5.3 × 10−5 7.057215508699 × 10−6 9.28636417763 × 10−10
0.50 0.00000000000000 5.0 × 10−5 1.312806175591734 × 10−17 6.396931098296924 × 10−17
0.60 −0.21850801222441 5.2 × 10−5 7.057215509060 × 10−6 9.28636584296 × 10−10
0.70 −0.41562693777745 5.4 × 10−5 1.3423564604320 × 10−5 1.731384868808 × 10−9
0.80 −0.57206140281768 5.3 × 10−5 1.7467564658924 × 10−5 8.9389752599978 × 10−8
0.90 −0.67249851196396 5.5 × 10−5 1.2923596913073 × 10−5 8.6007680755174 × 10−8
1.00 −0.70710678118655 5.4 × 10−5 3.952393967665557 × 10−14 2.198241588757810 × 10−14

exact solution, approximate solution vi (t ) for several values of i and errors obtained by MOL II are given in Tables 5 and 6
and Figs. 5–7 respectively.
 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188 2183

Table 5
The norm kvexact − vMOL kxi ,∞ in Example 2

xi MOL I MOL II

0.10 3.204176544324699 × 10 −5
1.984972012314401 × 10−7
0.20 4.558331840376351 × 10−5 1.679581814739706 × 10−7
0.30 4.182611406278181 × 10−5 1.219577727695764 × 10−7
0.40 2.458643023162122 × 10−5 1.137402102918683 × 10−7
0.50 1.242518910410138 × 10−14 2.811945523278168 × 10−14
0.60 2.458643024139118 × 10−5 1.137402084530614 × 10−7
0.70 4.182611405878500 × 10−5 1.219577708821973 × 10−7
0.80 4.558331840875951 × 10−5 1.679581880242864 × 10−7
0.90 3.204176544702175 × 10−5 1.984971881308084 × 10−7
1.00 3.066435994014682 × 10−13 3.248512570053208 × 10−13

Table 6
The norm kvexact − vMOL k∞ in Example 2
MOL I MOL II

4.558331840875951 × 10 −5
1.984972012314401 × 10−7

Fig. 8. Plot of the exact solution in Example 3.

Fig. 9. Plot of the vi (t ) for i = 10k, k = 1, . . . , 10 in Example 3.

3.3. Test 3

In this example, we apply the proposed methods in this paper to find the solution of (1.1)–(1.5) with l = 1, T = 4,

q(x, t ) = −2(x − t ) exp(−x − t ),

f1 (x) = 0, f2 (x) = x exp(−x),
g1 (t ) = 0, g2 (t ) = −2t exp(−t − 1) + t exp(−t ).
2184 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

Fig. 10. Plot of the v(xi , t ) − vi (t ) for (a1) xi = 1

4
, (a2) xi = 1
2
, (a3) xi = 3
4
, (a4) xi = 1 in Example 3.

Table 7
The norm kvexact − vMOL kxi ,∞ in Example 3

xi MOL I MOL II

0.10 8.790697957105 × 10 −7
1.3567209786181 × 10−9
0.20 1.2140191431882 × 10−6 1.2268780713587 × 10−9
0.30 1.1133469118643 × 10−6 9.826361113685 × 10−10
0.40 7.043162150097 × 10−7 6.474974688364 × 10−10
0.50 4.225513185602 × 10−7 1.327798704320 × 10−10
0.60 7.709645756804 × 10−7 7.196475637627 × 10−10
0.70 1.0072074030434 × 10−6 1.0519506238316 × 10−9
0.80 1.2092548875610 × 10−6 1.2563831364165 × 10−9
0.90 1.1058150297394 × 10−6 1.3559769695970 × 10−9
1.00 1.3828775877189 × 10−6 2.4905782158857 × 10−9

v(x, t ) = xt exp(−(x + t )) is the exact solution of this problem which can be readily verified. Such as the previous examples,
the results are shown in Tables 7 and 8 and Figs. 8–10.

3.4. Test 4

In the fourth example, we solve the problem (1.1)–(1.5) with l = 1, T = 5,

q(x, t ) = 2x5 + 2x3 − 2x2 − (20x3 + 6x − 2)(t 2 − t ),
f1 (x) = 0, f2 (x) = −x5 − x3 + x2 ,
t (t − 1)
g1 (t ) = 0, g2 (t ) = .
12
 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188 2185

Table 8
The norm kvexact − vMOL k∞ in Example 3
MOL I MOL II

1.3828775877189 × 10−6 2.4905782158857 × 10−9

Fig. 11. Plot of the exact solution in Example 4.

Fig. 12. Plot of the vi (t ) for i = 5k, k = 1, . . . , 10 in Example 4.

Fig. 13. Plot of the vi (t ) for i = 5k, k = 11, . . . , 20 in Example 4.

The exact solution of this equation is

v(x, t ) = (x5 + x3 − x2 )(t 2 − t ).
Some numerical results are summarized in Tables 9 and 10 and Figs. 11–14.
2186 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

Fig. 14. Plot of the v(xi , t ) − vi (t ) for (a1) xi = 1

4
, (a2) xi = 1
2
, (a3) xi = 3
4
, (a4) xi = 1 in Example 4.

Table 9
The norm kvexact − vMOL kxi ,∞ in Example 4

xi MOL I MOL II

0.10 1.4589452381 × 10−4 2.27205158021 × 10−8

0.20 2.7172977818 × 10−4 1.59301172808 × 10−8
0.30 3.7488135133 × 10−4 2.98653940467 × 10−8
0.40 4.3322433196 × 10−4 4.39071972114 × 10−8
0.50 4.1659346493 × 10−4 5.65441311551 × 10−8
0.60 2.9994415140 × 10−4 7.65543484160 × 10−8
0.70 5.996250296 × 10−5 9.44377208656 × 10−8
0.80 3.3337581084 × 10−4 1.134464433505 × 10−7
0.90 8.9980978039 × 10−4 1.326242955457 × 10−7
1.00 1.66653416461 × 10−3 2.4346786631213 × 10−6

Table 10
The norm kvexact − vMOL k∞ in Example 4
MOL I MOL II

1.66653416461 × 10−3 2.4346786631213 × 10−6

The MOL approach usually enables us to solve quite general and complicated partial differential equations relatively
easily and with acceptable efficiency. It is applicable to a wide range of problems in many areas [47]. The reader can see the
Appendix of [22] for some problems in physics, fluid dynamics, reactor models, automatic control, and more.
 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188 2187

4. Conclusion

The method of lines is generally recognized as a comprehensive and powerful approach to the numerical solution of
time-dependent partial differential equations. This method proceeds in two separate steps: Spatial derivatives are first
replaced with finite difference, finite volume, finite element or other algebraic approximations and then the resulting
system of ordinary differential equations which is usually stiff, is integrated in time. The success of this method is
explained by the availability of high-quality numerical algorithms for the solution of stiff systems of ODEs. This paper
investigated MOL approach for solving the one-dimensional hyperbolic equation with an integral condition [1,48]. Also
the new algorithms were tested on several problems. The computational results confirmed the efficiency, reliability and
accuracy of this procedure. It is worth pointing out that the solutions for MOL 1 and MOL 2 (second- and fourth-order
finite differences, respectively) clearly demonstrate the improved accuracy of the higher-order (MOL 2) method. Also, this
superior performance is achieved with very little increased computational effort (since the fourth-order finite difference
approximation of MOL 2 requires only five terms in a weighted sum compared with three terms in MOL 1). Thus, clearly the
use of higher-order finite differences in the MOL solution of the wave equation, and generally for other partial differential
equations, is recommended. Finally we would like to mention that one issue of future work is to develop similar method
to solve the parabolic inverse problems investigated in [49]. Other task is to apply the technique proposed in the current
research to find the solutions of the problems studied in [50–53].

Acknowledgments

The authors are very grateful to the reviewer A for carefully reading the paper and for his(her) constructive comments
and suggestions which have improved the paper.

References

[1] M. Dehghan, On the solution of an initial-boundary value problem that combines Neumann and integral condition for the wave equation, Numer.
Methods Partial Differential Equations 21 (2005) 24–40.
[2] D.E. Carlson, Linear thermoelasticity, in: Encyclopedia of Physics, vol. 2, Springer, Berlin, 1972.
[3] W.A. Day, A decreasing property of solutions of a parabolic equation with applications to thermoelasticity and other theories, Quart. Appl. Math. 4l
(1983) 468–475.
[4] A.A. Samarskii, Some problems in modern theory of differential equations, Differents. Uravn. l6 (1980) 122l–1228.
[5] J.H. Cushman, T.R. Ginn, Nonlocal dispersion irregular media with continuously evolving scales of heterogeneity, Transport Porous Media l3 (1993)
123–138.
[6] J.H. Cushman, H. Xu, F. Deng, Nonlocal reactive transport with physical and chemical heterogeneity: Localization error, Water Resource Res. 31 (1995)
2219–2237.
[7] A. Bouziani, Mixed problem for certain nonclassical equations containing a small parameter, Acad. Roy. Belg. Bull. CI, Sci. 6 (5) (1994) 389–400.
[8] A. Bouziani, Strong solution for a mixed problem with nonlocal condition for certain pluriparabolic equations, Hiroshima Math. J. 27 (3) (1997)
373–390.
[9] I.S. Gordeziani, G.A. Avalishvili, On the constructing of solutions of the-nonlocal initial boundary problems for one-dimensional medium oscillation
equations, Matem. Modelirovanie 12 (1) (2000) 94–103.
[10] N.I. Kavalloris, D.S. Tzanetis, Behaviour of critical solutions of a nonlocal hyperbolic problem in ohmic heating of foods, Appl. Math. E-Notes 2 (2002)
59–65.
[11] S. Mesloub, A. Bouziani, On a class of singular hyperbolic equation with a weighted integral condition, Internat. J. Math. Math. Sci. 22 (1999) 511–520.
[12] L.S. Pulkina, On solvability in L2 of nonlocal problem with integral conditions for hyperbolic equations, Differets. Uravn. VN 2 (2000) 1–6.
[13] M. Dehghan, A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications, Numer. Methods
Partial Differential Equations 22 (2006) 220–257.
[14] M. Dehghan, Implicit collocation technique for heat equation with non-classic initial condition, Int. J. Non-Linear Sci. Numer. Simul. 7 (2006) 447–450.
[15] S.A. Beilin, Existence of solutions for one-dimensional wave equations with nonlocal conditions, Electron. J. Differential Equations 76 (2001) 1–8.
[16] A. Saadatmandi, M. Dehghan, Numerical solution of the one-dimensional wave equation with an integral condition, Numer. Method Partial Differential
Equations 23 (2007) 282–292.
[17] W.T. Ang, A numerical method for the wave equation subject to a non-local conservation condition, Appl. Numer. Math. 56 (2006) 1054–1060.
[18] M. Ramezani, M. Dehghan, M. Razzaghi, Combined finite difference and spectral methods for the numerical solution of hyperbolic equation with an
integral condition, Numer. Methods Partial Differential Equations 24 (2008) 1–8.
[19] M. Dehghan, M. Lakestani, The use of cubic B-spline scaling functions for solving the one-dimensional hyperbolic equation with a nonlocal
conservation condition, Numer. Methods Partial Differential Equations 23 (2007) 1277–1289.
[20] G. Hall, J.M. Watt (Eds.), Modern Numerical Methods for Ordinary Differential Equations, Clarendon Press, Oxford, 1976.
[21] A.M. Loeb, W.E. Schiesser, Stiffness and accuracy in the method of lines integration of partial differential equations, in: Proc. of the 1973 Summer
Computer Simulation Conf., 2, 1973, pp. 25–39.
[22] W.E. Schiesser, The Numerical Method of Lines, Academic Press, New York, 1991.
[23] S. Hamdi, W.E. Schiesser, G.W. Griffiths, Method of lines, from scholarpedia. http://www.scholarpedia.org/article/Method_of_Lines.
[24] S. Hamdi, W.H. Enright, W.E. Schiesser, J.J. Gottlieb, Exact solutions and conservation laws for coupled generalized Korteweg deVries and quintic
regularized long wave equations, Nonlinear Anal. 63 (2005) 1425–1434.
[25] W.E. Schiesser, Method of lines solution of the Korteweg–de Vries equation, Comput. Math. Appl. 28 (1994) 147–154.
[26] M. Dehghan, Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices, Math. Comput.
Simulation 71 (2006) 16–30.
[27] T. Koto, Method of lines approximations of delay differential equations, Comput. Math. Appl. 48 (2004) 45–59.
[28] A.G. Bratsos, The solution of the two-dimensional sine-Gordon equation using the method of lines, J. Comput. Appl. Math. 206 (2007) 251–277.
[29] S. Hamdi, W.H. Enright, Y. Ouellet, W.E. Schiesser, Method of lines solutions of the extended Boussinesq equations, J. Comput. Appl. Math. 183 (2005)
327–342.
[30] P. Saucez, A.V. Wouwer, W.E. Schiesser, P. Zegeling, Method of lines study of nonlinear dispersive waves, J. Comput. Appl. Math. 168 (2004) 413–423.
[31] J. Taler, P. Duda, Solution of non-linear inverse heat conduction problems using the method of lines, Heat Mass Transfer 37 (2001) 147–155.
[32] H. Han, Z. Huang, The direct method of lines for the numerical solution of interface problem, Comput. Methods Appl. Mech. Engrg. 171 (1999) 61–75.
2188 F. Shakeri, M. Dehghan / Computers and Mathematics with Applications 56 (2008) 2175–2188

[33] V.D. Egorov, The method of lines for condensation kinetics problem solving, J. Aerosol Sci. 30 (Suppl. I) (1999) S247–S248.
[34] G. Torres, C. Turner, Method of straight lines for a Bingham problem in cylindrical pipes, Appl. Numer. Math. 47 (2003) 543–558.
[35] A. Erdem, S. Pamuk, The method of lines for the numerical solution of a mathematical model for capillary formation: The role of tumor angiogenic
factor in the extra-cellular matrix, Appl. Math. Comput. 186 (2007) 891–897.
[36] I. Ayranci, N. Selçuk, MOL solution of DOM for transient radiative transfer in 3-D scattering media, J. Quantitative Spectroscopy Radiative Transfer 84
(2004) 409–422.
[37] V.R. Subramanian, R.E. White, Semi-analytical method of lines for solving elliptic partial differential equations, Chem. Eng. Sci. 59 (2004) 781–788.
[38] A.G. Bratsos, The solution of the Boussinesq equation using the method of lines, Comput. Methods Appl. Mech. Eng. 157 (1998) 33–44.
[39] P. Saucez, W.E. Schiesser, A. Vande Wouwer, Upwinding in the method of lines, Math. Comput. Simulation 56 (2001) 171–185.
[40] S. Pamuk, A. Erden, The method of lines for the numerical solution of a mathematical model for capillary formation: The role of endothelial cells in
the capillary, Appl. Math. Comput. 186 (2007) 831–835.
[41] P. Saucez, A. Vande Wouwer, P.A. Zegeling, Adaptive method of lines solutions for the extended fifth-order Korteweg deVries equation, J. Comput.
Appl. Math. 183 (2005) 343–357.
[42] Amr A. Sharaf, H.O. Bakodah, A good spatial discretisation in the method of lines, Appl. Math. Comput. 171 (2005) 1253–1263.
[43] N. Selcuk, T. Tarhan, S. Tanrikulu, Comparison of method of lines and finite difference solutions of 2-D Navier–Stokes equations for transient laminar
pipe flow, Internat. J. Numer. Methods Eng. 53 (2002) 1615–1628.
[44] P.A. Ramachandran, Method of lines with boundary elements for transient diffusion-reaction problems, Numer. Methods Partial Differential Equations
22 (2006) 831–846.
[45] I. Ahmad, M. Berzins, MOL solvers for hyperbolic PDEs with source terms, Math. Comput. Simulation 56 (2001) 115–125.
[46] J.R. Cash, Efficient time integrators in the numerical method of lines, J. Comput. Appl. Math. 183 (2005) 259–274.
[47] H. Ramos, J. Vigo-Aguiar, An almost L-stable BDF-type method for the numerical solution of stiff ODEs arising from the method of lines, Numer.
Methods Partial Differential Equations 23 (2007) 1110–1121.
[48] M. Dehghan, A. Saadatmandi, Variational iteration method for solving the wave equation subject to an integral conservation condition, Chaos Solitons
Fractals, 2008 (in press).
[49] M. Dehghan, Parameter determination in a partial differential equation from the overspecified data, Math. Comput. Modelling 41 (2005) 196–213.
[50] M. Dehghan, M. Tatari, Identifying an unknown function in a parabolic equation with overspecified data via He’s variational iteration method, Chaos
Solitons Fractals 36 (2008) 157–166.
[51] M. Dehghan, Efficient techniques for the second-order parabolic equation subject to nonlocal specifications, Appl. Numer. Math. 52 (2005) 39–62.
[52] M. Dehghan, Identification of a time-dependent coefficient in a partial differential equation subject to an extra measurement, Numer. Methods for
Partial Differential Equations 21 (2005) 611–622.
[53] M. Dehghan, The one-dimensional heat equation subject to a boundary integral specification, Chaos Solitons Fractals 32 (2007) 661–675.