# A New Numerical Method of Particular Solutions for Inhomogeneous Burgers’ Equation

Authors: Huantian Xie; Dingfang Li; Feng Li
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Based on the finite difference scheme in time, the method of particular solutions using radial basis functions is proposed to solve one-dimensional time-dependent inhomogeneous Burgers’ equations. Two numerical examples with good accuracy are given to validate the proposed method.


---

## Body

## 1. Introduction

In this paper we consider the one-dimensional nonlinear evolutionary partial differential equation(1)∂u(x,t)∂t+u(x,t)∂u(x,t)∂x-1R∂2u(x,t)∂x2=f(x,t).
The corresponding homogeneous equation
(2)∂u(x,t)∂t+u(x,t)∂u(x,t)∂x-1R∂2u(x,t)∂x2=0
was first introduced by Bateman [1] who considered its steady state solutions. Later, Burgers [2, 3] treated it as a mathematical model for free turbulence, and subsequently this equation is widely referred to as Burgers’ equation. Burgers’ equation can model several physical phenomena such as traffic, shock waves, and continuous stochastic processes. It can also be used to test various numerical algorithms. Due to its wide range of applicability, several researchers have been interested in the properties of its solution.Burgers’ equation has been solved analytically for a restricted set of arbitrary initial and boundary conditions [4, 5]. Benton and Platzman [6] surveyed about 35 distinct exact solutions of the one-dimensional Burgers-like equations and their classifications. It is well known that the exact solution of Burgers’ equation can only be computed for restricted values of R. Therefore, various numerical methods were employed to obtain the solution of Burgers’ equation. It is not our purpose to exhaust all these numerical schemes. Nevertheless, the solution methodologies commonly fall into the following classes: finite difference method (FDM), finite element method (FEM), and spectral methods. A survey of these techniques is given in [7, 8]. The previous numerical methods all depend on the mesh of the studied domain [9–12].To alleviate the difficulty of mesh generation, various meshless techniques have been introduced during the past two decades. In a meshless (meshfree) method, a set of scattered nodes are selected in the computational domain. Meshless schemes include the method of fundamental solutions (MFS) [13], the method of particular solutions (MPS) [14–17], the element-free Galerkin method [18], local point interpolation [12], and boundary knot method [19]. It is known that the MFS is a boundary-type meshless method which is highly accurate for solving homogeneous equations if the fundamental solution of the given differential operator is known [20]. However, the fundamental solution of a given differential equation is not always available and often very difficult to derive. The ill-conditioning of the matrix resulting from the formulation of using the MFS and the location of source points are still outstanding research problems. To extend the MFS to inhomogeneous equations or time-dependent problems, the MPS has been introduced to evaluate the particular solution of the given differential equation. Since the particular solution is not unique, there is a rich variety of numerical techniques developed for this purpose.Radial basis functions (RBFs), polynomial functions, trigonometric functions, and so forth [20–24], have been employed as the basis functions to approximate the particular solutions for the given differential equation. Once a particular solution has been evaluated, the given inhomogeneous equations can be reduced to the homogeneous equation. The original differential equation can be recovered by adding the homogeneous solution and the particular solution. This is a two-stage numerical scheme and is a well-known procedure for solving linear partial differential equations.In general, the fundamental solution can be viewed as a special type of particular solution. When the inhomogeneous term is replaced by the delta function, a particular solution becomes a fundamental solution. The main idea of the MFS is that a fundamental solution satisfies the homogeneous equation inside the domain, and one only needs to enforce the fundamental solution on the boundary conditions to obtain the solution of the given homogeneous problem. Motivated by a similar idea, the particular solution can be used to solve inhomogeneous problems; that is, since the particular solution satisfies the given inhomogeneous equation through the domain without satisfying the boundary conditions, one only needs to impose the boundary conditions to obtain the solution of the given inhomogeneous problem.It is the purpose of this paper to extend the MPS to a one-stage numerical scheme for solving one-dimensional time-dependent Burgers’ equations through the use of RBFs and then solve the time-dependent problems without the need of a two-stage numerical scheme [25] for obtaining homogeneous solution.


## 2. The Method of Particular Solutions

Consider the following boundary value problem:(3)Δu(x)=f(x),x∈Ω,(4)Bu(x)=b(x),x∈Γ,
where Δ and B are the Laplace operator and boundary differential operators, respectively, Ω⊆ℝd is the solution domain, Γ is its boundary, and f(x) and b(x) are given functions.Approximatef(x) by a finite series of RBFs {ϕj}1N through interpolation, and interpolants f^n(x) to f(x) can be constructed as
(5)f(x)≃f^N(x)=∑j=1N‍λjϕ(rj),
in which rj=∥x-xj∥ is Euclidean distance, {xj}1N is a set of interpolation points, and the real coefficients {λj}1N are to be determined by solving
(6)∑j=1N‍λjϕ(rj)=f(xk),1≤k≤N
if the n×n real coefficient matrix Aϕ=(ϕ(rkj)) is positive definite.Therefore, from (6) an approximate particular solution u(x) to (3) is given by
(7)u(x)≃∑j=1N‍λjΦ(rj),
where Φ(rj) is obtained analytically by solving
(8)ΔΦ(rj)=ϕ(rj).
If we impose u(x) in (7) to satisfy the governing equation in (3) and boundary conditions in (4), then u(x) becomes an approximate solution of the original partial differential equations (3)-(4). To be more specific, we have
(9)∑j=1N‍λjϕ(rj)=f(xk),xk∈Ω,∑j=1N‍λjBΦ(rj)=b(xk),xk∈Γ.
For the numerical implementation, we let {xk}1nI be the interior points, {xk}nI+1N the boundary points, and N=nI+nB. From (9) we have
(10)∑j=1N‍λjϕ(rj)=f(xk),1≤k≤nI,∑j=1N‍λjBΦ(rj)=b(xk),nI+1≤k≤N.
The above system of equations can be easily solved using a standard matrix solver. Once the {λj}1N are determined, the approximate particular solution becomes the approximate solution u(x) of (3)-(4); that is,
(11)u(x)=∑j=1N‍λjΦ(rj).
Note that an accurate approximation of the particular solution u(x) depends on the appropriate choice of radial basis function ϕ. In the RBF literature [20, 26], some of the globally defined RBFs are only conditionally positive definite [26]. The unique solvability of the interpolation problem can be obtained by adding a polynomial term to the interpolation (5), giving
(12)f(x)≃f^n(x)=∑j=1N‍λjϕ(rj)+∑k=1K‍μkpk(x),
along with the constraints
(13)∑j=1N‍λjpk(xj)=0,1≤k≤K,
where {pk}1K is a basis of 𝒫m-1, the space of d-variate polynomials of order not exceeding m-1, and
(14)K=(m-1+dd)
is the dimension of 𝒫m-1.There are many types of globally defined RBFs [20], and the most popular RBFs areInverse multiquadric (IMQ),(15)ϕ(r)=(r2+c2)q/2,q<0,c>0,Multiquadric (MQ),(16)ϕ(r)=(r2+c2)q/2,q∈ℝ+∖2ℤ,c>0,m>[q2].Gaussian (G),(17)φ(r)=e-cr2,c>0.Polyharmonic (PH),(18)ϕ(r)=rq,q∈ℝ+∖2ℤ,c>0,m>[q2].Polyharmonic (PH),(19)ϕ(r)=rqlog(r),q∈2ℤ,m>q2.


## 3. Methodology

In order to solve such whole-space problems by numerical methods, we limit our consideration to a finite subdomain[a,b]. In other words, function u(x,t) satisfies the following general nonlinear one-dimensional time-dependent Burgers’ equation:
(20)∂u(x,t)∂t=1R∂2u(x,t)∂x2-u(x,t)∂u(x,t)∂x+f(x,t)(x,t)∈[a,b]×[0,T],
with the initial condition
(21)u(x,0)=g(x),x∈[a,b],
and the Dirichlet boundary condition
(22)u(a,t)=p(t),u(b,t)=q(t),t∈[0,T],
where R>0 is interpreted as the Reynolds number, and 1/R is the kinematic viscosity, and f(x,t),  g(x),  p(t), and q(t) are known functions.In the following section, a generalized trapezoidal method (θ-method) is used to approximate the time derivative in (20). Let τ=tn+1-tn be the time step and tn=nτ. For any tn≤t≤tn+1 and 0≤θ≤1,  u(x,t) can be approximated as follows:
(23)u(x,t)≃θu(x,tn+1)+(1-θ)u(x,tn),f(x,t)≃θf(x,tn+1)+(1-θ)f(x,tn).
Then,
(24)Δu(x,t)≃θΔu(x,tn+1)+(1-θ)Δu(x,tn),∂u(x,t)∂t≃u(x,tn+1)-u(x,tn)τ,
where Δ=∂2/∂x2. For simplicity, we denote un≡u(x,tn),  uxn≡∂u(x,tn)/∂x, and fn≡f(x,tn). Substituting (23)–(24) into (20)–(22), we obtain the following equation:
(25)un+1-unτ=θ[1RΔun+1-un+1uxn+1+fn+1]+(1-θ)[1RΔun-unuxn+fn].
Rewrite the previous equation as follows:
(26)un+1-θτ[1RΔun+1-un+1uxn+1+fn+1]=un+(1-θ)τ[1RΔun-unuxn+fn].
Then,
(27)Δun+1=Rθτun+1+Run+1uxn+1-Rfn+1-Rθτun-R(1-θ)θ[1RΔun-unuxn+fn].Assume thatun+1(x) is a sought solution to the elliptic PDE. We can represent the right hand side of (27) as a function F(x). This means that (27) is a standard Poisson-type differential equation
(28)Δun+1=F(x).
Therefore, if the fictitious function F(x) is known, (27) is equivalent to the Poisson-type equation (28) under the same boundary conditions.Uniformly chooseN-2 collocation points {xk}2N-1 in the interior of domain [a,b] and two boundary points x=a and x=b. For implementation, let x1=a,  xN=b. Approximating the function F(x) by RBFs {ϕj}1N, we have
(29)F(x)≃∑j=1N‍λjn+1ϕj(rkj),k,j=1,2,…,N,
where rkj=∥(xk-xj)2∥. Then we can approximate u(x) at time step n+1 as follows:
(30)un+1(x)≃∑j=1N‍λjn+1Φj(rj),
where Φj(x) is obtained by solving
(31)ΔΦ(x)=ϕ(x).
We use two RBFs, namely, IMQ and PH, where
(32)ϕj(x)=1(rj2+c2)3,ϕk(x)=rk3.
It is easy to obtain the following Φ(x):
(33)Φj(x)=1c2rj2+c2,Φk(x)=120rk5.Note that (27) is a recursion formula, and we can solve each elliptic PDE step by step starting with initial condition (21). As is well known, it is difficult to obtain an accurate numerical derivative from scattered data. Therefore, we choose θ=1 in our method in order to avoid evaluating Δun in (27). In this case, we can reformulate (27) as follows:
(34)un+1-τ[1RΔun+1-un+1uxn+1]=un+τfn+1.
The nonlinear term un+1uxn+1 is linearized as follows:
(35)un+1uxn+1=un+1uxn+unuxn+1-unuxn.
Substituting (35) into (34) and rearranging, we obtain
(36)un+1-τ[1RΔun+1-un+1uxn-unuxn+1]=un+τunuxn+τfn+1.Write (30) together with boundary condition (22) in matrix form
(37)[u]n=A[λ]n,
where [u]n=[u1n,u2n,…,uNn],  [λ]n=[λ1n,λ2n,…,λNn], and A=[akj],  1≤k,  j≤N. There are N-2 internal points and 2 boundary points. The N×N matrix A can be split into A=Ad+Ab, where
(38)Ad={(akj)if2≤k≤N-1,1≤j≤Nakj=0otherwise,Ab={(akj)ifk=1,N,1≤j≤Nakj=0otherwise.
Applying this to the domain points and boundary points, (37) and (22) can be reformulated in the following matrix form:
(39){TTAd-τɛ∇2Ad+τ[diag(Ad[λ]n)∇Ad+diag(∇Ad[λ]n)Ad]+AbTT}[λ]n+1=Ad[λ]n+τ(Ad[λ]n)·*(∇Ad[λ]n)+[F]n+1,
where ∇ is the gradient differential operator, diag(Ad[λ]n) is a diagonal matrix with Ad[λ]n as its main diagonal, [F]n+1=[pn+1,τf2n+1,…,τfN-1n+1,qn+1]T, and the accent “·*” means component by component multiplication of two vectors.


## 4. Numerical Results

Two different problems are used to test the accuracy of our method. In order to evaluate the numerical errors, we adopt three kinds of errors defined by(40)L∞=maxj|uj-Uj|,L2=∑j=1N‍|uj-Uj|2,RMS=1N∑j=1N‍|uj-Uj|2,
where uj=u(xj,T) is the exact analytical solution, and Uj is the numerical solution of uj.Example 1.
Consider the following nonlinear one-dimensional time-dependent Burgers’ equation with a large Reynolds numberR=10000, in the square domain [0,2π],
(41)∂u(x,t)∂t+u(x,t)∂u(x,t)∂x-1R∂2u(x,t)∂x2=-12sin(2x)e-2t/R,
with the initial condition
(42)u(x,0)=cos(x),
and the boundary condition
(43)u(x,t)=u(x+2π,t).
The analytical solution is given as
(44)u(x,t)=cos(x)e-t/R.
We choose two RBFs, namely, IMQ and PH, as defined in (32). The L∞,  L2, and root-mean-square (RMS) errors for our numerical solutions are shown in Table 1 for t=0.5,  1,  1.5, and 2. It can be seen that the accuracy of the method using IMQ is much higher than that using PH.
The graph of the analytical and estimated solutions fort=2 is shown in Figure 1. The absolute error graph is shown in Figure 2. We also show the space-time graph of the estimated solution in Figure 3.Table 1
L
∞, L2, and RMS errors, with τ=0.1, dx=π/10, and x∈ [0, 2π].
t
L
∞
L
2
RMS
IMQ
PH
IMQ
PH
IMQ
PH
0.5
9.6592
×
10
-
6
7.7103
×
10
-
3
1.3706
×
10
-
5
8.2398
×
10
-
3
3.0648
×
10
-
6
1.8425
×
10
-
3
1.0
9.9640
×
10
-
6
1.2424
×
10
-
2
1.5692
×
10
-
5
1.7040
×
10
-
2
3.5088
×
10
-
6
3.8102
×
10
-
3
1.5
9.9985
×
10
-
6
1.5138
×
10
-
2
1.4249
×
10
-
5
2.3674
×
10
-
2
3.1863
×
10
-
6
5.2937
×
10
-
3
2.0
1.7471
×
10
-
6
2.7952
×
10
-
2
2.7217
×
10
-
5
4.0602
×
10
-
2
6.0859
×
10
-
6
9.0788
×
10
-
3Figure 1
Analytical and numerical solution att=2, with IMQ as RBF, τ=0.1,  dx=π/10,  and x∈[0,2π], for Example 1.Figure 2
Absolute error att=2, with IMQ as RBF, τ=0.1,  dx=π/10,  and x∈[0,2π], for Example 1.Figure 3
Space-time graph of the solution up tot=2, with IMQ as RBF, τ=0.1,  dx=π/10,  and x∈[0,2π], for Example 1.Example 2.
We consider the second-order nonlinear Burgers’ equation(45)∂u(x,t)∂t+u(x,t)∂u(x,t)∂x-1R∂2u(x,t)∂x2=f(x,t),
where f(x,t) depends upon the exact solution of (45) as follows:
(46)u(x,t)=[sin(2πx)+12sin(πx)]et.
We take the required initial and boundary functions from the exact solution in the domain [0,1]. Similar to Example 1, we also choose two RBFs, namely, IMQ and PH, as shown in (32). In Table 2, we compute the L∞,  L2 errors and RMS errors at t=0.5,  1,  1.5, and 2, with the Reynolds number R=0.001. It indicates that the case of IMQ also has higher accuracy than the case of PH. In Figure 4, we draw the graph of the analytical and estimated solutions, and in Figure 5 the absolute error graph is shown. Figures 4 and 5 are shown at t=2. The space-time graph of the numerical solution is shown in Figure 6.Table 2
L
∞, L2, and RMS errors, with τ=0.1, dx=0.05, and x∈[0, 1].
t
L
∞
L
2
RMS
IMQ
PH
IMQ
PH
IMQ
PH
0.5
1.0490
×
10
-
4
6.2473
×
10
-
3
3.8990
×
10
-
4
1.9479
×
10
-
2
5.5140
×
10
-
5
2.7547
×
10
-
3
1.0
2.6814
×
10
-
4
1.0352
×
10
-
2
1.0921
×
10
-
3
3.2281
×
10
-
2
1.5445
×
10
-
4
4.5652
×
10
-
2
1.5
7.1806
×
10
-
4
1.7210
×
10
-
2
2.8510
×
10
-
3
5.3690
×
10
-
2
4.0320
×
10
-
4
7.5929
×
10
-
3
2.0
1.9485
×
10
-
3
2.8760
×
10
-
2
7.8048
×
10
-
3
8.9872
×
10
-
2
1.1038
×
10
-
3
1.2710
×
10
-
2Figure 4
Analytical and numerical solution att=2, with PH as RBF, τ=0.1,  dx=0.05,  and x∈[0,1], for Example 2.Figure 5
Absolute error att=2, with PH as RBF, τ=0.1,  dx=0.05,  and x∈[0,1], for Example 2.Figure 6
Space-time graph of the solution up tot=2, with PH as RBF, τ=0.1,  dx=π/10,  and x∈[0,2π], for Example 2.Remark 3.
Note that, for the two examples shown previously, the RBFs are unconditionally positive definite to guarantee the solvability of the resulting systems. However, some RBFs are conditionally positive definite. These types of interpolation problems can be obtained by adding a polynomial term to the interpolation (12), and it is also easy to verify the efficiency of the proposed schemes for these cases.


## 5. Conclusions

In this paper we proposed and implemented the method of particular solutions to solve the one-dimensional time-dependent Burgers’ equation. The effectiveness of the computational scheme is well demonstrated. It must be emphasized that the choice of radial basis functions is a flexible feature of these methods. The radial basis functions can be globally supported, infinitely differentiable and contain free parameters, namely, shape parameters, which affect both accuracy of the solutions and conditioning of the collocation matrix. The optimal shape parameters in Examples1 and 2 using IMQ and PH, respectively, for all the calculations were found experimentally. The optimal choice of the shape parameters in RBFs is still an outstanding research problem [27–29]. A similar approach can be extended to solving 2D or 3D time-dependent partial differential equations. These research topics will be the focus of future investigation.



---

*Source: 974808-2013-07-17.xml*