Extended B-Spline Collocation Method For KdV-Burgers Equation

0.56 MB·

by  Ozlem Ersoy Hepson

Extended B-Spline Collocation Method For
KdV-Burgers Equation
Ozlem Ersoya, Alper Korkmazb,∗, Idiris Dagc
a
DepartmentofMathematics&Computer,EskisehirOsmangaziUniversity,26480,Eskisehir,Turkey.
b
7 DepartmentofMathematics,C¸ankırıKaratekinUniversity,18200,C¸ankırı,Turkey.
1 c
0 DepartmentofComputerEngineering,EskisehirOsmangaziUniversity,26480,Eskisehir,Turkey.
2
January 12, 2017
n
a
J

1
1

Abstract
]
A
The extended form of the classical polynomial cubic B-spline function is
N
usedtosetupacollocationmethodforsomeinitialboundaryvalueproblems
.
h derivedfortheKorteweg-deVries-Burgersequation. Havingnonexistenceof
t
a third order derivatives of the cubic B-splines forces us to reduce the order of
m
thetermu togiveacoupledsystemofequations. Thespacediscretization
xxx
[ of this system is accomplished by the collocation method following the time

1 discretization with Crank-Nicolson method. Two initial boundary value
v problems, one having analytical solution and the other is set up with a non
3
analytical initial condition, have been simulated by the proposed method.
9
8
2 Keywords: KdV-Burgers Equation; Extended cubic B-spline; collocation; motion
0
of waves.
.
1
0
7
1 1 Introduction
:
v
i Consider the Korteweg-de Vries - Burgers (KdVB) equation of the form
X
r
a u +εuu −ϑu +µu = 0 a ≤ x ≤ b (1)
t x xx xxx
where u = u(x,t), subscripts denote the derivatives and ε,ν and µ are real co-
efficients with the property ενµ (cid:54)= 0 [1–3]. The KdVB equation for dispersive
∗[email protected]
1
and dissociative media is derived in various cases such as existence and absence
of a complete system of eigenvector or for waves in plasma with Hall dispersion
and Joule dissipation by Ruderman [4]. Including both viscosity and inertia terms
simultaneously causes long gravity waves to be governed by the KdVB equation.
This result was obtained by deriving the balance equations for an incompressible
viscous fluid starting from the column model approximation in the study of Bampi
and Morro [5].
The Weierstrass P-function solution, which can be expressed in terms of Jacobi
cosine functions, for the KdVB 1 was found by Kalinowski and Grundland [3]. In
the study, the Riemann invariants for non-homogenous systems of the first order
partial differential equations method was implemented to obtain that solution.
Brugarino and Pantano [6] obtaned the Jacobi elliptic type solution of the nonho-
mogenous KdVB equation with variable coefficients. Some travelling solitary wave
solutions to compound KdVB equation was developed by using automated method
based on an ansatz of a polynomial in terms of a tanh function [7]. A particular
analytical solution for the KdVB equation was also obtained by using variable
transformations and proofs of theorems [8]. The solutions in the stationary waves
form were determined for the generalized KdVB equation containing nonconser-
vative terms of linear pump, linear HF dissipation, and nonlinear dissipation [9].
Malkov [10] constructed the asymptotic travelling wave solution, representing a
shock-train, of the KdVB equation driven by the long scale periodic driver. The
stability of travelling wave solutions of the KdVB equation was studied for various
values of the viscosity coefficient α [11].
The exact solutions consisting of powers of some hyperbolic functions of the KdVB
equation were derived by homogenous balance technique in [12]. These solutions
arethecombinationofabelly-shapeandkink-typesolitarywaves. Zhao[13]imple-
mented the hyperbolic function method and the Wu elimination technique for the
newtype solitarywave andperiodic solutionscontainingsomeblow-upsolutions of
the KdVB. After reducing the KdVB equation to an ordinary differential equation
by the compatible wave transformation, he determined the coefficients and param-
eters in the predicted solution by substituting it into this equation. Yuanxi and
Jiashi [14] developed many solitary wave solutions for the KdVB equation by the
superpositionmethodbasedontheanalysisonthefeaturesoftheBurgers,theKdV
and the KdV-Burgers equations. A generalized tanh function method based on
Riccati equation was derived to obtain multiple soliton solutions containing some
trigonometric, hyperbolic and complex functions for the KdVB equation [15]. Soli-
man [16] obtained the tanh-type, coth-type exact solutions of the equation by the
modified extended tanh method. tanh-type solution for the KdVB equation was
obtained by using G(cid:48)/G expansion method [17]. Some compact and non-compact
solutions for variants of the KdVB equation were constructed by Wazwaz [18].
2
Kudryashov showed that many of the solutions of the KdVB equation listed above
can be converted to each other by using [19]. He also showed that the traveling
wave solutions of the Fisher equation and the KdVB are in the same form and
derived an exponential-type solution with the assistance of Weierstrass function
for the KdVB equation.
Besides proposed various exact solutions in the literature summarized above, some
numerical techniques were also developed to solve some problems constructed with
KdVB equation. One of the frontier numerical studies on KdVB equation is based
on Bubnov-Galerkin’s finite element method using cubic B-splines [20]. The tem-
poral evaluation of a Maxwellian and the time evaluation of the solutions were
studied in details in that study. A collocation method based on quintic B-spline
functions were also derived for the numerical solutions of the KdVB equation [21].
A linear Galerkin-Fourier spectral technique was implemented for the numerical
solutionsoftheKdVBequtionwithperiodicinitialcondition[22]. Aspectralcollo-
cation method based on differentiated Chebyshev polynomials combined with the
fourth order Runge-Kutta method was proposed for the numerical solution of the
initial boundary value problem modeling tanh-type solitary wave solutions [23].
Unconditionally stable septic B-spline collocation method was proposed to solve
the KdVB equation numerically by El-Danaf [24].
In this study, we consider some initial boundary value problems defined as
u +εuu −ϑu +µu = 0 a < x < b (2)
t x xx xxx
subject to the initial condition
u(x,0) = f(x), a ≤ x ≤ b
in the finite problem interval [a,b]. We choose the boundary conditions from the
set
u(a,t) = g (t), u(b,t) = g (t),
1 2
u (a,t) = g (t), u (b,t) = g (t),
x 3 x 4
u (a,t) = g (t), u (b,t) = g (t),
xx 5 xx 6
u (a,t) = g (t), u (b,t) = g (t)
xxx 7 xxx 8
where g (t),i = 1,2,...,8 denote x independent functions.
i
2 Method of Solution
2.1 Extended Cubic B-splines
Let π be a uniform grid distribution of the finite interval [a,b] such as π : a = x <
0
x < ... < x = b with equal sub interval length h = (b−a)/N. The extended
1 N
3
form of a cubic B-spline E is defined as [25,26]
i

4h(1−λ)(x−x )3 +3λ(x−x )4, [x ,x ],
 i−2 i−2 i−2 i−1

1  (−41−2hλ()xh−4 +xi1−21h)33(−x3−λ(xxi−−1)x+i−61h)42(2+λ)(x−xi−1)2 [xi−1,xi],
E (x) = (4−λ)h4 −12h3(x−x )+6h2(2+λ)(x−x )2
i 24h4 i+1 i+1 [x ,x ],
 +12h(x−x )3 −3λ(x−x )4 i i+1
 i+1 i+1

 4h(λ−1)(x−x )3 +3λ(x−x )4, [x ,x ],
 i+2 i+2 i+1 i+2

 0 otherwise.
(3)
where the real λ is the free parameter. In fact, classical cubic polynomial B-splines
are the particular form of extended cubic B-splines when λ = 0. The set of the
extended cubic B-splines {E (x)}N+1 is a basis for the functions defined in the
i i=−1
interval [x ,x ] [25,26]. When the free parameter is chosen different from zero,
0 N
the shape of the function changes slightly, Fig 1. The nonzero functional and
derivative values of each extended cubic B-spline E (x) at the grids of πin the
i
[a,b] can be summarized as in Table 1.
Figure 1: Extended B-splines for various values of the free parameter λ
4
Table 1: E (x) and its derivatives at the grid points
i
x x x x x x
i−2 i−1 i i+1 i+2
24E (x) 0 4−λ 16+2λ 4−λ 0
i
2hE(cid:48)(x) 0 −1 0 1 0
i
2h2E(cid:48)(cid:48)(x) 0 2+λ −4−2λ 2+λ 0
i
Having only first two derivatives of each extended B-spline E (x) forces us to
i
reduce the derivative order of the KdVB equation (1). Thus, defining a new
variable v(x,t) = u (x,t) reduces the KdVB equation (1) to a coupled system of
x
equations
u +εuv −ϑv +µv = 0
t x xx
(4)
u −v = 0
x
2.2 Collocation Method
Let U(x,t) and V(x,t) be the approximate solutions to u(x,t) and v(x,t), respec-
tively that
N+1
(cid:88)
U(x,t) = δ E (x),
i i
i=−1
(5)
N+1
(cid:88)
V(x,t) = φ E (x)
i i
i=−1
where δ and φ are time dependent parameters that will be calculated via the
i i
extended cubic B-splines and the complementary conditions. The approximate
solutions given in (5) and their first and second derivatives at a grid x can be
i
determined by using the Table (1). Thus, the approximate solutions and their
derivatives take the form
4−λ 8+λ 4−λ
U = U(x ,t) = δ + δ + δ ,
i i i−1 i i+1
24 12 24
−1
U(cid:48) = U(cid:48)(x ,t) = (δ −δ )
i i 2h i−1 i+1
2+λ
U(cid:48)(cid:48) = U(cid:48)(cid:48)(x ,t) = (δ −2δ +δ )
i i 2h2 i−1 i i+1
(6)
4−λ 8+λ 4−λ
V = V(x ,t) = φ + φ + φ
i i i−1 i i+1
24 12 24
−1
V(cid:48) = V(cid:48)(x ,t) = (φ −φ )
i i 2h i−1 i+1
2+λ
V(cid:48)(cid:48) = V(cid:48)(cid:48)(x ,t) = (φ −2φ +φ )
i i 2h2 i−1 i i+1
5
To integrate the system (4) in time variable we first use forward finite difference
and the Crank-Nicolson methods to give
Un+1 −Un (UV)n+1 +(UV)n (V )n+1 +(V )n (V )n+1 +(V )n
x x xx xx
+ε −ϑ +µ = 0
∆t 2 2 2
Un+1 +Un Vn+1 +Vn
x x − = 0
2 2
(7)
where Un+1 = U(x,(n+1)∆t) represent the solution at the (n+1).th time level.
Here tn+1 = tn +t, and ∆t is the time step, superscripts denote n.th time level ,
tn = n∆t. Linearizing the terms (UV)n+1and (UV)n in (7) as
(UV)n+1 = Un+1Vn +UnVn+1 −UnVn
(8)
(UV)n = UnVn
gives
Un+1 −Un Un+1Vn +UnVn+1 Vn+1 +Vn Vn+1 +Vn
+ε( )−ϑ( x x )+µ( xx xx) = 0
∆t 2 2 2
Un+1 +Un Vn+1 +Vn
x x − = 0
2 2
(9)
Substitution of the approximate solutions (5) and their derivatives (6) into (9)
yields
ν δn+1 +ν φn+1 +ν δn+1 +ν φn+1 +ν δn+1 +ν φn+1 (10)
m1 m−1 m2 m−1 m3 m m4 m m1 m+1 m5 m+1
= ν δn +ν φn +ν δn +ν φn +ν δn +ν φn
m6 m−1 m7 m−1 m8 m m9 m m6 m+1 m10 m+1
ν δn+1 +ν φn+1 +0δn+1 +ν φn+1 +ν δn+1 −ν φn+1 (11)
m11 m−1 m12 m−1 m m13 m m11 m+1 m12 m+1
= −ν δn −ν φn +0δn −ν φn −ν δn +ν φn
m11 m−1 m12 m−1 m m13 m m11 m+1 m12 m+1
The coefficients in the system (10) and (11) are
(cid:18) (cid:19)
2 2
ν = +εL α ν = α
m1 1 m8 2
∆t ∆t
ν = εKα −ϑβ +µγ ν = −µγ
m2 1 1 1 m9 2
(cid:18) (cid:19)
2
ν = +εL α ν = −ϑβ −µγ
m3 2 m10 1 1
∆t
ν = εKα +µγ ν = β
m4 2 2 m11 1
ν = εKα +ϑβ +µγ ν = −α
m5 1 1 1 m12 1
2
ν = α ν = −α
m6 1 m13 2
∆t
ν = ϑβ −µγ
m7 1 1
6
where
K = α δ +α δ +α δ
1 i−1 2 i 1 i+1
L = α φ +α φ +α φ
1 i−1 2 i 1 i+1
and
4−λ 8+λ
α = , α =
1 2
24 12
1 2+λ 4+2λ
β = − , γ = , γ = −
1 2h 1 2h2 2 2h2
The system (10-11) can be written in the following matrices system
Axn+1 = Bxn (12)
where
 
ν ν ν ν ν ν
m1 m2 m3 m4 m1 m5
ν ν 0 ν ν −ν 
m11 m12 m13 m11 m12
 
 ν ν ν ν ν ν 
m1 m2 m3 m4 m1 m5
 
A = νm11 νm12 0 νm13 νm11 −νm12 
 
 ... ... ... ... ... ... 
 
 
ν ν ν ν ν ν
 m1 m2 m3 m4 m1 m5 
ν ν 0 ν ν −ν
m11 m12 m13 m11 m12
and
 
ν ν ν ν ν ν
m6 m7 m8 m9 m6 m10
−ν −ν 0 −ν −ν ν 
m11 m12 m13 m11 m12
 
 ν ν ν ν ν ν 
m6 m7 m8 m9 m6 m10
 
B = −νm11 −νm12 0 −νm13 −νm11 νm12 
 
 ... ... ... ... ... ... 
 
 
ν ν ν ν ν ν
 m6 m7 m8 m9 m6 m10
−ν −ν 0 −ν −ν ν
m11 m12 m13 m11 m12
There are 2N +2 equations with 2N +6 unknown parameters
xn+1 = (δn+1,φn+1,δn+1,φn+1...,δn+1,φn+1)
−1 −1 0 0 N+1 N+1
in this system. A unique solution of the system can be obtained by imposing the
boundary conditions U (a,t) = 0,U (b,t) = 0,V (a,t) = 0,V (b,t) = 0 to have the
x x x x
following the equations:
δ = δ , φ = φ , δ = δ , φ = φ
−1 1 −1 1 N−1 N+1 N−1 N+1
7
Elimination of the parameters δ ,φ ,δ ,φ , using the equations (10) from
−1 −1 N+1 N+1
the the system gives a solvable system of 2N+2 linear equations including 2N+2
unknown parameters.
Since the right handside of this equation consists of only known values at the n.th
time level, the solution of the system gives the solution at the (n + 1).th time
level. In order to initialize this iteration system, we need the initial vector x0.The
initial parameter vectors d = (δ ,δ ,..δ ,δ ), d = (φ ,φ ,..φ ,φ ) can
1 −1 0 N N+1 2 −1 0 N N+1
be determined by using
U (a,0) = 0 = γ δ0 +γ δ0 +γ δ0,
xx 1 −1 2 0 1 1
U (x ,0) = γ δ0 +γ δ0 +γ δ0 = U (x ,0),i = 1,...,N −1
xx i 1 i−1 2 i 1 i+1 xx i
U (b,0) = 0 = γ δ0 +γ δ0 +γ δ0 ,
xx 1 N−1 2 N 1 N+1
V (a,0) = 0 = φ0 −φ0
x −1 1
V (x ,0) = φ0 −φ0 = V (x ,0),i = 1,...,N −1
x i i−1 i+1 x i
V (b,0) = φ0 −φ0
x N−1 N+1
3 Numerical tests
We chose some initial boundary value problems constructed on the the KdV-
Burgers equation to check the accuracy and validity of the proposed method. The
discrete maximum error norm
(cid:12) (cid:12)
L = |u−U| = max(cid:12)u −Un(cid:12)
∞ ∞ j j
j
between the numerical and exact solution is computed for different values of the
viscosity parameter ϑ, different times with various time and space step sizes ∆t,
h.
3.1 Example 1: Initial boundary value problem with ana-
lytical solution
Consider the initial boundary value problem constructed with the KdV-Burgers’
equation (1) with the initial condition
6ϑ2 (cid:20) (cid:18) ϑx (cid:19) 1 (cid:18)ϑx (cid:19)(cid:21)
u(x,0) = − 1+tanh + sech2 0µ (13)
25µ 10µ 2 1
and the boundary conditions
6ϑ2 (cid:20) (cid:18) ϑ 6ϑ2 (cid:19) 1 (cid:18) ϑ 6ϑ2 (cid:19)(cid:21)
u(x,a) = − 1+tanh (a+ t) + sech2 (a+ t)
25µ 10µ 25µ 2 10µ 25µ
6ϑ2 (cid:20) (cid:18) ϑ 6ϑ2 (cid:19) 1 (cid:18) ϑ 6ϑ2 (cid:19)(cid:21)
u(x,b) = − 1+tanh (b+ t) + sech2 (b+ t)
25µ 10µ 25µ 2 10µ 25µ
8
These complementary conditions have been derived from the Fan and Zhang’s [27]
analytical solution of the form
6ϑ2 (cid:20) (cid:18) ϑ (cid:18) 6ϑ2 (cid:19)(cid:19) 1 (cid:18) ϑ (cid:18) 6ϑ2 (cid:19)(cid:19)(cid:21)
u(x,t) = − 1+tanh x+ t + sech2 x+ t
25µ 10µ 25µ 2 10µ 25µ
(14)
for the particular choice of ε = 1. This solution has three components, a constant,
the tanh component and the sech component and it represents a traveling wave
moving along the x−axis as time goes.
Forthesakeofcomparisonwithsomeearlierworks,wefixtheparametersµ = 0.01,
h = 0.5, ∆t = 0.001 and run the algorithm in the finite problem interval [−20,20]
uptotheterminatingtimet = 1. Whenthecoefficientischosenasϑ = 0.001,0.005
and 0.01, the terminating time profiles of the waves are given in Fig 2(a)-2(c). The
increase of the dispersion parameter value dominates the Burgers-type solution
and the wave gets steeper. The maximum error distrubitions for each case are
also plotted in Fig 3(a)-3(c) for the same parameter values. The error distribution
plots in all cases show that the error is greater near the descent part of the wave.
A comparison of discrete maximum norms with the results obtained by the radial
basiscollocationmethodswithmultiquadric(MQ),inversequadric(IQ)andGaus-
sian (GA) formsin the study [28] at thetime t = 1 andt = 10 istabulated in Table
2. For compatibility we chose the parameters as ε = 1, ϑ = 0.004, ∆t = 0.001 and
h = 0.5. This comparison is an indicator of the efficiency of the proposed method
based on the extended B-spline functions. Both the particular selections of the
extention parameter of the extended B-splines as λ = 0 and λ = −1.969 gives
at least two decimal digits better results than the results of radial basis colloca-
tion methods. Even the proposed method for both values of extension parameter
generate the results in the same decimal digit accuracy at the time t = 1, the
nonzero value of the extension parameter keeps the accuracy digits at t = 10 but
the classical cubic B-spline based method lose one decimal digit in accuracy. The
method derived from radial basis functions can not catch the same accuracy even
though the same parameters are used during the implementation of the proposed.
Table 2: Comparison of the maximum error with radial basis collocation methods
t Present (λ = 0) Present (λ = −1.969) MQ [28] GA [28] IQ [28]
1 9.441×10−11 7.271×10−11 6.822×10−9 7.913×10−9 4.077×10−7
10 1.269×10−10 9.509×10−11 2.479×10−8 3.294×10−6 1.270×10−6
9
3.2 Example 2: Non analytical initial data
In the second problem, we consider the particular form of the KdVB when ϑ =
0. Eliminating the dispersion term reduces the KdVB to the Korteweg-de Vries
equation (KdVE). In order to check the validity of the method, we chose a non
analytical solution describing the split of the initial pulse into a train of waves.
This initial pulse is described by the data
(cid:20) (cid:21)
1 |x|−25
u(x,0) = 1−tanh (15)
2 5
and the routine is run up to the terminating time t = 800 in the finite interval
[−50,150] for the compatibility to the results obtained in [29–31]. During this long
process time we will observe the lowest four conserved quantities [30,32]
∞
(cid:90)
C = udx
1
−∞
∞
(cid:90)
C = u2dx
2
−∞
∞
(cid:90) (cid:18) (cid:19)
3µ
C = u3 − (u(cid:48))2 dx
3
ε
−∞
∞
(cid:90) (cid:18) µ µ2 (cid:19)
C = u4 −12 uu2 +7.2 u2 dx
4 ε x ε2 xx
−∞
tobecalculatedforthisinitialboundaryvalueproblem. Theapproximatevaluesof
the conserved quantities are computed by modified trapezoid rule that the mean of
thefunctionalvaluesattheconsecutivepointsintheproblemintervalinsteadofthe
functionalvaluesatthegridpoints. Theparametersarechosenasε = 0.2, µ = 0.1,
∆t = 0.05, h = 0.4 in accordance with the earlier study of Korkmaz [30]. The
conserved quantities for various choices of the extension parameter λ are compared
with the ones obtained by the differential quadrature methods based on cosine
expansion (CDQ) and the Lagrange polynomials (LPDQ) combined with fourth
orderRunge-Kuttatechniqueatvariousdiscretetimesduringthesimulation,Table
3.
It can be concluded that the lowest two quantities C and C that do not con-
1 2
tain any derivatives of the solution are conserved successfully for all values of the
extension parameter λ even its zero. They both are in good agreement with the
ones reported in Korkmaz’ study [30]. Although the third lowest quantity C is
3
10