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

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.

1. Introduction

In this paper we consider the one-dimensional nonlinear evolutionary partial differential equation 𝜕𝑢(𝑥,𝑡)()()𝜕𝑡+𝑢𝑥,𝑡𝜕𝑢𝑥,𝑡−1𝜕𝑥𝑅𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=𝑓(𝑥,𝑡). (1) The corresponding homogeneous equation 𝜕𝑢(𝑥,𝑡)()()𝜕𝑡+𝑢𝑥,𝑡𝜕𝑢𝑥,𝑡−1𝜕𝑥𝑅𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=0 (2) 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 𝑅. 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)(4) where Δ and 𝐵 are the Laplace operator and boundary differential operators, respectively, Ω⊆ℝ^𝑑 is the solution domain, Γ is its boundary, and 𝑓(𝐱) and 𝑏(𝐱) are given functions.

Approximate 𝑓(𝐱) by a finite series of RBFs {𝜙_𝑗}^𝑁₁ through interpolation, and interpolants 𝑓_𝑛(𝐱) to 𝑓(𝐱) can be constructed as 𝑓(𝐱)≃𝑓_𝑁(𝐱)=^𝑁_𝑗=1𝜆_𝑗𝜙𝑟_𝑗, (5) in which 𝑟_𝑗=‖𝐱−𝐱_𝑗‖ is Euclidean distance, {𝐱_𝑗}^𝑁₁ is a set of interpolation points, and the real coefficients {𝜆_𝑗}^𝑁₁ are to be determined by solving ^𝑁_𝑗=1𝜆_𝑗𝜙𝑟_𝑗=𝑓𝐱_𝑘,1≤𝑘≤𝑁 (6) if the 𝑛×𝑛 real coefficient matrix 𝐴_𝜙=(𝜙(𝑟_𝑘𝑗)) is positive definite.

Therefore, from (6) an approximate particular solution 𝑢(𝐱) to (3) is given by 𝑢(𝐱)≃^𝑁_𝑗=1𝜆_𝑗Φ𝑟_𝑗, (7) where Φ(𝑟_𝑗) is obtained analytically by solving ΔΦ𝑟_𝑗=𝜙𝑟_𝑗. (8) If we impose 𝑢(𝐱) in (7) to satisfy the governing equation in (3) and boundary conditions in (4), then 𝑢(𝐱) becomes an approximate solution of the original partial differential equations (3)-(4). To be more specific, we have ^𝑁_𝑗=1𝜆_𝑗𝜙𝑟_𝑗=𝑓𝐱_𝑘,𝐱_𝑘∈Ω,_𝑁𝑗=1𝜆_𝑗𝐵Φ𝑟_𝑗=𝑏𝐱_𝑘,𝐱_𝑘∈Γ. (9) For the numerical implementation, we let {𝐱_𝑘}^𝑛𝐼₁ be the interior points, {𝐱_𝑘}^𝑁_𝑛𝐼+1 the boundary points, and 𝑁=𝑛_𝐼+𝑛_𝐵. From (9) we have ^𝑁_𝑗=1𝜆_𝑗𝜙𝑟_𝑗=𝑓𝐱_𝑘,1≤𝑘≤𝑛_𝐼,_𝑁𝑗=1𝜆_𝑗𝐵Φ𝑟_𝑗=𝑏𝐱_𝑘,𝑛_𝐼+1≤𝑘≤𝑁. (10) The above system of equations can be easily solved using a standard matrix solver. Once the {𝜆_𝑗}^𝑁₁ are determined, the approximate particular solution becomes the approximate solution 𝑢(𝐱) of (3)-(4); that is, 𝑢(𝐱)=^𝑁_𝑗=1𝜆_𝑗Φ𝑟_𝑗. (11) Note that an accurate approximation of the particular solution 𝑢(𝐱) 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 𝑓(𝐱)≃𝑓_𝑛(𝐱)=^𝑁_𝑗=1𝜆_𝑗𝜙𝑟_𝑗+^𝐾_𝑘=1𝜇_𝑘𝑝_𝑘(𝐱), (12) along with the constraints ^𝑁_𝑗=1𝜆_𝑗𝑝_𝑘𝐱_𝑗=0,1≤𝑘≤𝐾, (13) where {𝑝_𝑘}^𝐾₁ is a basis of 𝒫_𝑚−1, the space of 𝑑-variate polynomials of order not exceeding 𝑚−1, and 𝑑𝐾=𝑚−1+𝑑 (14) is the dimension of 𝒫_𝑚−1.

There are many types of globally defined RBFs [20], and the most popular RBFs are Inverse multiquadric (IMQ), 𝜙(𝑟)=𝑟²+𝑐^2𝑞/2,𝑞<0,𝑐>0, (15) Multiquadric (MQ), 𝜙(𝑟)=𝑟²+𝑐^2𝑞/2,𝑞∈ℝ⁺𝑞⧵2ℤ,𝑐>0,𝑚>2. (16) Gaussian (G), 𝜑(𝑟)=𝑒^−𝑐𝑟2,𝑐>0. (17) Polyharmonic (PH), 𝜙(𝑟)=𝑟^𝑞,𝑞∈ℝ⁺𝑞⧵2ℤ,𝑐>0,𝑚>2. (18) Polyharmonic (PH), 𝜙(𝑟)=𝑟^𝑞(𝑟)𝑞log,𝑞∈2ℤ,𝑚>2. (19)

3. Methodology

In order to solve such whole-space problems by numerical methods, we limit our consideration to a finite subdomain [𝑎,𝑏]. In other words, function 𝑢(𝑥,𝑡) satisfies the following general nonlinear one-dimensional time-dependent Burgers’ equation: 𝜕𝑢(𝑥,𝑡)=1𝜕𝑡𝑅𝜕²𝑢(𝑥,𝑡)𝜕𝑥²()()−𝑢𝑥,𝑡𝜕𝑢𝑥,𝑡()()∈[]×[],𝜕𝑥+𝑓𝑥,𝑡𝑥,𝑡𝑎,𝑏0,𝑇 (20) with the initial condition 𝑢(𝑥,0)=𝑔(𝑥),𝑥∈[𝑎,𝑏], (21) and the Dirichlet boundary condition 𝑢(𝑎,𝑡)=𝑝(𝑡),𝑢(𝑏,𝑡)=𝑞(𝑡),𝑡∈[0,𝑇], (22) where 𝑅>0 is interpreted as the Reynolds number, and 1/𝑅 is the kinematic viscosity, and 𝑓(𝑥,𝑡),  𝑔(𝑥),  𝑝(𝑡), and 𝑞(𝑡) are known functions.

In the following section, a generalized trapezoidal method (𝜃-method) is used to approximate the time derivative in (20). Let 𝜏=𝑡^𝑛+1−𝑡^𝑛 be the time step and 𝑡^𝑛=𝑛𝜏. For any 𝑡^𝑛≤𝑡≤𝑡^𝑛+1 and 0≤𝜃≤1,  𝑢(𝑥,𝑡) can be approximated as follows: 𝑢(𝑥,𝑡)≃𝜃𝑢𝑥,𝑡^𝑛+1+(1−𝜃)𝑢𝑥,𝑡^𝑛,𝑓(𝑥,𝑡)≃𝜃𝑓𝑥,𝑡^𝑛+1+(1−𝜃)𝑓𝑥,𝑡^𝑛. (23) Then, Δ𝑢(𝑥,𝑡)≃𝜃Δ𝑢𝑥,𝑡^𝑛+1+(1−𝜃)Δ𝑢𝑥,𝑡^𝑛,𝜕𝑢(𝑥,𝑡)≃𝑢𝜕𝑡𝑥,𝑡^𝑛+1−𝑢𝑥,𝑡^𝑛𝜏, (24) where Δ=𝜕²/𝜕𝑥². For simplicity, we denote 𝑢^𝑛≡𝑢(𝑥,𝑡^𝑛),  𝑢^𝑛_𝑥≡𝜕𝑢(𝑥,𝑡^𝑛)/𝜕𝑥, and 𝑓^𝑛≡𝑓(𝑥,𝑡^𝑛). Substituting (23)–(24) into (20)–(22), we obtain the following equation: 𝑢^𝑛+1−𝑢^𝑛𝜏1=𝜃𝑅Δ𝑢^𝑛+1−𝑢^𝑛+1𝑢_𝑥^𝑛+1+𝑓^𝑛+1+(11−𝜃)𝑅Δ𝑢^𝑛−𝑢^𝑛𝑢^𝑛_𝑥+𝑓^𝑛. (25) Rewrite the previous equation as follows: 𝑢^𝑛+11−𝜃𝜏𝑅Δ𝑢^𝑛+1−𝑢^𝑛+1𝑢_𝑥^𝑛+1+𝑓^𝑛+1=𝑢^𝑛+()𝜏11−𝜃𝑅Δ𝑢^𝑛−𝑢^𝑛𝑢^𝑛_𝑥+𝑓^𝑛. (26) Then, Δ𝑢^𝑛+1=𝑅𝑢𝜃𝜏^𝑛+1+𝑅𝑢^𝑛+1𝑢_𝑥^𝑛+1−𝑅𝑓^𝑛+1−𝑅𝑢𝜃𝜏^𝑛−𝑅(1−𝜃)𝜃1𝑅Δ𝑢^𝑛−𝑢^𝑛𝑢^𝑛_𝑥+𝑓^𝑛. (27)

Assume that 𝑢^𝑛+1(𝑥) is a sought solution to the elliptic PDE. We can represent the right hand side of (27) as a function 𝐹(𝑥). This means that (27) is a standard Poisson-type differential equation Δ𝑢^𝑛+1=𝐹(𝑥). (28) Therefore, if the fictitious function 𝐹(𝑥) is known, (27) is equivalent to the Poisson-type equation (28) under the same boundary conditions.

Uniformly choose 𝑁−2 collocation points {𝑥_𝑘}₂^𝑁−1 in the interior of domain [𝑎,𝑏] and two boundary points 𝑥=𝑎 and 𝑥=𝑏. For implementation, let 𝑥₁=𝑎,  𝑥_𝑁=𝑏. Approximating the function 𝐹(𝑥) by RBFs {𝜙_𝑗}^𝑁₁, we have 𝐹(𝑥)≃^𝑁_𝑗=1𝜆_𝑗^𝑛+1𝜙_𝑗𝑟_𝑘𝑗,𝑘,𝑗=1,2,…,𝑁, (29) where 𝑟_𝑘𝑗=‖(𝑥_𝑘−𝑥_𝑗)²‖. Then we can approximate 𝑢(𝑥) at time step 𝑛+1 as follows: 𝑢^𝑛+1(𝑥)≃^𝑁_𝑗=1𝜆_𝑗^𝑛+1Φ_𝑗𝑟_𝑗, (30) where Φ_𝑗(𝑥) is obtained by solving ΔΦ(𝑥)=𝜙(𝑥). (31) We use two RBFs, namely, IMQ and PH, where 𝜙_𝑗(𝑥)=1𝑟²_𝑗+𝑐²³,𝜙_𝑘(𝑥)=𝑟³_𝑘. (32) It is easy to obtain the following Φ(𝑥): Φ_𝑗(𝑥)=1𝑐²𝑟²_𝑗+𝑐²,Φ_𝑘(𝑥)=1𝑟20^5𝑘. (33)

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 Δ𝑢^𝑛 in (27). In this case, we can reformulate (27) as follows: 𝑢^𝑛+11−𝜏𝑅Δ𝑢^𝑛+1−𝑢^𝑛+1𝑢_𝑥^𝑛+1=𝑢^𝑛+𝜏𝑓^𝑛+1. (34) The nonlinear term 𝑢^𝑛+1𝑢_𝑥^𝑛+1 is linearized as follows: 𝑢^𝑛+1𝑢_𝑥^𝑛+1=𝑢^𝑛+1𝑢^𝑛_𝑥+𝑢^𝑛𝑢_𝑥^𝑛+1−𝑢^𝑛𝑢^𝑛_𝑥. (35) Substituting (35) into (34) and rearranging, we obtain 𝑢^𝑛+11−𝜏𝑅Δ𝑢^𝑛+1−𝑢^𝑛+1𝑢^𝑛_𝑥−𝑢^𝑛𝑢_𝑥^𝑛+1=𝑢^𝑛+𝜏𝑢^𝑛𝑢^𝑛_𝑥+𝜏𝑓^𝑛+1. (36)

Write (30) together with boundary condition (22) in matrix form [𝑢]^𝑛=𝐴[𝜆]^𝑛, (37) where [𝑢]^𝑛=[𝑢^𝑛₁,𝑢^𝑛₂,…,𝑢^𝑛_𝑁],  [𝜆]^𝑛=[𝜆^𝑛₁,𝜆^𝑛₂,…,𝜆^𝑛_𝑁], and 𝐴=[𝑎_𝑘𝑗],  1≤𝑘,  𝑗≤𝑁. There are 𝑁−2 internal points and 2 boundary points. The 𝑁×𝑁 matrix 𝐴 can be split into 𝐴=𝐴_𝑑+𝐴_𝑏, where 𝐴_𝑑=𝑎_𝑘𝑗𝑎if2≤𝑘≤𝑁−1,1≤𝑗≤𝑁_𝑘𝑗𝐴=0otherwise,_𝑏=𝑎_𝑘𝑗𝑎if𝑘=1,𝑁,1≤𝑗≤𝑁_𝑘𝑗=0otherwise. (38) Applying this to the domain points and boundary points, (37) and (22) can be reformulated in the following matrix form: 𝐀_𝑑−𝜏𝜀∇²𝐀_𝑑+𝜏diag𝐀_𝑑[𝜆]^𝑛∇𝐀_𝑑+diag∇𝐀_𝑑[𝜆]^𝑛𝐀_𝑑+𝐀_𝑏[𝜆]^𝑛+1=𝐀_𝑑[𝜆]^𝑛+𝜏𝐀_𝑑[𝜆]^𝑛⋅∗∇𝐀_𝑑[𝜆]^𝑛+[𝐹]^𝑛+1, (39) where ∇ is the gradient differential operator, diag(𝐀_𝑑[𝜆]^𝑛) is a diagonal matrix with 𝐀_𝑑[𝜆]^𝑛 as its main diagonal, [𝐹]^𝑛+1=[𝑝^𝑛+1,𝜏𝑓₂^𝑛+1,…,𝜏𝑓^𝑛+1_𝑁−1,𝑞^𝑛+1]^𝑇, 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 𝐿_∞=max_𝑗𝑢_𝑗−𝑈_𝑗,𝐿₂=^𝑁_𝑗=1𝑢_𝑗−𝑈_𝑗²,RMS=1𝑁^𝑁_𝑗=1𝑢_𝑗−𝑈_𝑗², (40) where 𝑢_𝑗=𝑢(𝑥_𝑗,𝑇) is the exact analytical solution, and 𝑈_𝑗 is the numerical solution of 𝑢_𝑗.

Example 1. Consider the following nonlinear one-dimensional time-dependent Burgers’ equation with a large Reynolds number 𝑅=10000, in the square domain [0,2𝜋], 𝜕𝑢(𝑥,𝑡)()()𝜕𝑡+𝑢𝑥,𝑡𝜕𝑢𝑥,𝑡−1𝜕𝑥𝑅𝜕²𝑢(𝑥,𝑡)𝜕𝑥²1=−2sin(2𝑥)𝑒^{−2𝑡/𝑅}, (41) with the initial condition 𝑢(𝑥,0)=cos(𝑥), (42) and the boundary condition 𝑢(𝑥,𝑡)=𝑢(𝑥+2𝜋,𝑡). (43) The analytical solution is given as 𝑢(𝑥,𝑡)=cos(𝑥)𝑒^−𝑡/𝑅. (44) We choose two RBFs, namely, IMQ and PH, as defined in (32). The 𝐿_∞,  𝐿₂, and root-mean-square (RMS) errors for our numerical solutions are shown in Table 1 for 𝑡=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 for 𝑡=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 

𝐿_∞, 𝐿₂, and RMS errors, with 𝜏=0.1, 𝑑_𝑥=𝜋/10, and 𝑥∈ [0, 2𝜋].

Figure 1 

Analytical and numerical solution at 𝑡=2, with IMQ as RBF, 𝜏=0.1,  𝑑𝑥=𝜋/10,  and 𝑥∈[0,2𝜋], for Example 1.

Figure 2 

Absolute error at 𝑡=2, with IMQ as RBF, 𝜏=0.1,  𝑑𝑥=𝜋/10,  and 𝑥∈[0,2𝜋], for Example 1.

Figure 3 

Space-time graph of the solution up to 𝑡=2, with IMQ as RBF, 𝜏=0.1,  𝑑𝑥=𝜋/10,  and 𝑥∈[0,2𝜋], for Example 1.

Example 2. We consider the second-order nonlinear Burgers’ equation 𝜕𝑢(𝑥,𝑡)()()𝜕𝑡+𝑢x,𝑡𝜕𝑢𝑥,𝑡−1𝜕𝑥𝑅𝜕²𝑢(𝑥,𝑡)𝜕𝑥²=𝑓(𝑥,𝑡), (45) where 𝑓(𝑥,𝑡) depends upon the exact solution of (45) as follows: 𝑢()=()+1𝑥,𝑡sin2𝜋𝑥2sin(𝜋𝑥)𝑒^𝑡. (46) 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 𝐿_∞,  𝐿₂ errors and RMS errors at 𝑡=0.5,  1,  1.5, and 2, with the Reynolds number 𝑅=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 𝑡=2. The space-time graph of the numerical solution is shown in Figure 6.

Table 2 

𝐿_∞, 𝐿₂, and RMS errors, with 𝜏=0.1, 𝑑_𝑥=0.05, and 𝑥∈[0, 1].

Figure 4 

Analytical and numerical solution at 𝑡=2, with PH as RBF, 𝜏=0.1,  𝑑𝑥=0.05,  and 𝑥∈[0,1], for Example 2.

Figure 5 

Absolute error at 𝑡=2, with PH as RBF, 𝜏=0.1,  𝑑𝑥=0.05,  and 𝑥∈[0,1], for Example 2.

Figure 6 

Space-time graph of the solution up to 𝑡=2, with PH as RBF, 𝜏=0.1,  𝑑𝑥=𝜋/10,  and 𝑥∈[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 Examples 1 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.

Acknowledgments

The work has been supported by the National Natural Science Foundation of China (Grants nos. 51190094, 61271337, 11201211, and 11201212) and AMEP of Linyi University. Special thanks are given to the reviewers for constructive advice and thoughtful comments.