Solution of the Coupled Burgers Equation Based on Operational Matrices of d-Dimensional Orthogonal Functions

(1)

Solution of the Coupled Burgers Equation Based on Operational Matrices of d-Dimensional Orthogonal Functions

Saeed Kazem^a, Malihe Shaban^b, and Jamal Amani Rad^c

a Department of Mathematics, Imam Khomeini International University, Ghazvin 34149-16818, Iran

b Department of Physics, Shahid Beheshti University, Evin, Tehran 19839, Iran

c Department of Computer Sciences, Shahid Beheshti University, Evin, Tehran 19839, Iran Reprint requests to S. K.; E-mail:saeedkazem@gmail.com

Z. Naturforsch.67a,267 – 274 (2012) / DOI: 10.5560/ZNA.2012-0026 Received October 7, 2011

This paper aims to construct a general formulation for thed-dimensional orthogonal functions and their derivative and product matrices. These matrices together with the Tau method are utilized to reduce the solution of partial differential equations (PDEs) to the solution of a system of algebraic equations. The proposed method is applied to solve homogeneous and inhomogeneous two-dimensional parabolic equations. Also, the mentioned method is employed to find the solution of the coupled Burgers equation. Illustrative examples are included to demonstrate the validity and applicability of the presented technique.

Key words:Coupled Burgers Equation; Operational Matrix; Chebyshev and Legendre Polynomials;

Tau Method.

Mathematics Subject Classification 2000:34A08

1. Introduction

The Burgers equation has been found to describe various kinds of phenomena such as a mathematical model of turbulence [1] and the approximate theory of flow through a shock wave travelling in viscous fluid [2]. The coupled Burgers system was derived by Esipov [3]. It is a simple model of sedimentation or evolution of scaled volume concentrations of two kinds of particles in fluid suspensions or colloids under the effect of gravity [4]. Using the Hopf–Cole transforma- tion, Fletcher [5] gave an analytical solution for the system of the two-dimensional Burgers equations. Sev- eral numerical methods for solving this equation have been given [6–13].

In this study, the system of coupled Burgers equations is investigated by applying the operational matrix. Accordingly, the orthogonal functions and their derivative and product matrices for d-dimensional time-depended partial differential equations are constructed. These matrices together with the Tau method are then utilized to reduce the solution of Burgers equations to the solution of a system of algebraic equations. The Tau approach is an approximation tech-

nique introduced by Lonczos [14] in 1938 to solve differential equations. The Tau method is based on expanding the required approximate solution as the elements of a complete set of orthogonal functions.

This method may be viewed as a special case of the so-called Petrov–Galerkin method. But, unlike the Galerkin approximation, the expansion functions are not required to satisfy the boundary constraint individ- ually [15–18].

2. Implementation of the Tau Method on the Heat Equation

In this section, we aim to convert the time-depended partial differential equations (PDEs) in the form of

u_t(x,t) =∇u(x,t) +f(x,t) inΩ×J, u(x,0) =g(x), x∈Ω,

Bu(x,t) =h(x,t) on∂ Ω×J

(1)

to a system of algebraic equations by applying the operational matrix of orthogonal functions in Ω×J.

Here,∇=∂²/∂x₁²+∂²/∂x²₂+· · ·+∂²/∂x²_d,Ω ⊆R^d is an open bounded domain with smooth boundary

c

2012 Verlag der Zeitschrift f¨ur Naturforschung, T¨ubingen·http://znaturforsch.com

(2)

∂ Ω,J= (0,T] withT >0, andg(x) andh(x,t)are known smooth functions.

In this way, we construct a family of orthogonal functions in Ω×J by using orthogonal polynomials such as Legendre and Chebyshev. At first, we assume that Ω =Ω1×Ω2× · · ·Ω_d, xi ∈Ωi ⊆R for i=1,2,· · ·,d, and the polynomials, denoted byφn^[i](x_i), are orthogonal with weight functionsw_i(x_i) overΩ_i, and alsoψn(t)is orthogonal with weight functionsw(t) overJ:

Z

Ω_i

φn^[i](x_i)φ_m^[i](x_i)w_i(x_i)dx_i=b^[i]nδ_nm, i=1,2,· · ·,d, Z T

0

ψ_n(t)ψ_m(t)w(t)dt=b_nδ_nm, (2) where δnm is the Kronecker function, b^[i]n andbn are known values for each polynomialφn^[i](x_i)andψn(t), respectively, that are obtained as

b^[i]_n = Z

Ωi

φn^[i](x_i)2

w_i(x_i)dx_i, i=1,2,· · ·,d, b_n=

Z _T 0

ψn(t)2

w(t)dt.

(3)

A function f_i(x_i)defined overΩican be expanded as fi(x_i) =

+∞

n=0

∑

a^[i]n φn^[i](x_i),

and the coefficientsa^[i]n are given by a^[i]n = 1

b^[i]_n Z

Ω_i

f_i(x_i)φ_n^[i](x_i)w_i(x_i)dx_i, n=0,1,2,· · ·. By defining

Φˆ^[i](x_i) =h

φ₀^[i](x_i),φ₁^[i](x_i),· · ·,φ_n^[i](x_i),· · ·iT

, Ψ(t) = [ψˆ ₀(t),ψ1(t),· · ·,ψ_n(t),· · ·]^T, the family of functions

Φˆ(x,t) =Φˆ^[1](x₁)⊗Φˆ^[2](x₂)⊗ · · · ⊗Φˆ^[d](x_d)⊗Ψˆ(t) are orthogonal with weight functions w(x,t) = w₁(x₁)w₂(x₂)· · ·w_d(x_d)w(t) over Ω×J; also ⊗ denotes the Kronecker product. It’s clear that the orthog- onality property appears as

Z

Ω×J

φˆn(x,t)φˆm(x,t)w(x,t)dxdt=c_nδnm. (4)

A functionu(x,t)can be expanded as u(x,t) =

+∞

n=0

∑

a_nφˆ_n(x,t), and the coefficientsa_nare given by

a_n= 1 cn

Z

Ω×J

u(x,t)φˆn(x,t)w(x,t)dxdt, n=0,1,2,· · ·.

(5) In practice, only the finite terms of ˆΦ(x,t)are con- sidered. Therefore,m_i-terms of ˆΦ^[i](x_i)ands-terms of ψ(t)ˆ are defined in the form

Φ^[i](x_i) =h

φ₀^[i](x_i),φ₁^[i](x_i),· · ·,φ_m^[i]

i−1(x_i)iT

,

Ψ(t) = [ψ₀(t),ψ1(t),· · ·,ψs−1(t)]^T, and consequently,N-terms of ˆφm(x)(N=s∏^d

i=1

m_i) can be obtained as

Φ(x,t) =Φ^[1](x₁)⊗Φ^[2](x₂)⊗ · · · ⊗Φ^[d](x_d)⊗Ψ(t).

These functions approximateu(x,t):

u_N(x,t) =

N−1

∑

i=0

a_iφ_i(x,t) =A^TΦ(x,t) with

A= [a₀,a₁,· · ·,aN−1]^T, φ_i(x,t) =φ_i^[1]

1 (x₁)φ_i^[2]

2 (x₂)· · ·φ_i^[d]

d (x_d)ψ_j(t),

i=i1m2m3· · ·m_ds+i2m3m4· · ·m_ds+· · · +id−1m_ds+i_ds+j,

(6)

i_k=

i m_k+1m_k+2· · ·m_ds

−i₁m₂m₃· · ·m_ks

−i₂m₃· · ·m_ks− · · · −ik−1m_ks, k=1,2,· · ·,d, (7) j=i−i₁m₂m₃· · ·m_ds+i₂m₃m₄· · ·m_ds+· · ·

+i_d−1m_ds+i_ds, (8)

where[·]denotes the integer part of the number.

2.1. Operational Matrix of Derivatives

At first, we describe a property of the Kronecker product as

(A₁⊗A2⊗ · · · ⊗An)(B₁⊗B2⊗ · · · ⊗Bn)

=A1B1⊗A2B2⊗ · · · ⊗AnBn. (9)

(3)

Lemma 1. The derivative operator of the vector Φ(x,t)can be expressed by

∂

∂xi

Φ(x,t)'D_iΦ(x,t), i=1,2,· · ·,d, where

D_i=I_m₁⊗I_m₂⊗ · · · ⊗I_m_i−1⊗D_i⊗I_m_i+1⊗ · · · ⊗I_m_d⊗I_s. D_iis the derivative matrix of vectorΦ^[i](x_i), and I_m_j is the identity matrix of order mj.

Proof.

∂

∂x_iΦ(x,t) =Φ^[1](x₁)⊗Φ^[2](x₂)⊗ · · ·

⊗Dx_iΦ^[i](x_i)⊗ · · · ⊗Φ^[d](x_d)⊗Ψ(t), where D_x_i is derivative operator of x_i. We assume D_x_iΦ^[i](x_i)'D_iΦ^[i](x_i)thatD_iis the derivative matrix of vectorΦ^[i](x_i).

Now by using (9), the Lemma can be proved.

Also, the time-derivative matrix is obtained as

∂

∂tΦ(x,t)'D_tΦ(x,t),

D_t=Im₁⊗I_m₂⊗ · · · ⊗Im_d⊗Dt,

whereD_tis the derivative matrix of vectorΨ(t).

Remark 1. The derivative matrix of vectorΦ^[i](x_i)for shifted Chebyshev and Legendre polynomials is obtained in [19–21]. In general, we can obtain the derivative matrix of vector Φ^[i](x_i)by using the orthogonal property (2)

D_i=<D_x_iΦ^[i](x_i),Φ^[i](x_i)^T>d⁻¹_i ,

where<Dx_iΦ^[i](x_i),Φ^[i](x_i)^T>andd_iare twomi×m_i matrices defined as

<D_x_iΦ^[i](x_i),Φ^[i](x_i)^T>

= Z

Ωi

D_x_iφn^[i](x_i)φ_m^[i](x_i)w_i(x_i)dx_i mi−1

n,m=0

, d_i=<Φ^[i](x_i),Φ^[i](x_i)^T>=diag{b^[i]_n}^m_n=0ⁱ⁻¹.

(10)

Denote thatφn^[i](x_i)are polynomials of degreen; one can write

φn^[i](x_i) =

n

∑

j=0

h^[i]_j,nx_i^j,

where

h^[i]_j,n= 1 j!

d^j

dx_i^jφn^[i](x_i) xi=0.

Subsequently,Dx_iφn^[i](x_i)is obtained as D_x_iφn^[i](x_i) =

n

∑

j=1

jh^[i]_j,nx_i^j−1.

Thus, by using above equation, the appeared integral in (10) can be obtained as

Z

Ωi

D^α_x

iφn^[i](x_i)φ_m^[i](x_i)w_i(x_i)dx_i

=

n

∑

j=1

jh^[i]_j,n Z

Ω_i

x_i^j−1φ_m^[i](x_i)w_i(x_i)dx_i.

In orthogonal polynomials one can say that the elements of the derivative matrixD_iare achieved as

D_i ^m_n,m=0ⁱ⁻¹ =











0, n−1<m,

n

∑

j=1 m

∑

k=0

j^h

[i]

j,nh^[i]_k,m b^[i]n

R

Ωix^k+j−1_i w_i(x_i)dx_i, n−1≥m.

For shifted Legendre polynomials defined in[0,L_i], we have

h^[i]_j,n=(−1)^j+n(j+n)!

(n−j)!(j!)²L_i^j , b^[i]_n = L_i

2n+1, w_i(x_i) =1, D_i ^mⁱ⁻¹

n,m=0=

( 0, n−1<morn+mis even,

4m+2

Li , otherwise.

Also, for shifted Chebyshev polynomials defined in [0,L_i], one has

h^[i]_j,n=n(−1)ⁿ⁻^j(n+j−1)!2^2j (n−j)!(2j)!L_i^j , b^[i]₀ =π, b^[i]_n =π

2, n≥1, w_i(x_i) = 1

pLix−x²,

D_i ^m_n,m=0ⁱ⁻¹ =







0, n−1<morn+mis even,

2n

L_i , m=0,

4n

L_i , m6=0.

Remark 2. The ∇-operational matrix can be constructed as

∇= ∂²

∂x²₁+ ∂²

∂x²₂+· · ·+ ∂²

∂x²_d,

∇=D²₁+D²₂· · ·+D²_d.

(4)

2.2. The Product Operational Matrix

The following property of the product of two orthogonal vectors will also be applied:

φ^[i](x_i)φ^[i](x_i)^TV 'V˜φ^[i](x_i),

where ˜V is anmi×mi product operational matrix for the vectorV. Using above equation and by the orthogonal property mentioned in (2), the elements{V˜_nm}^m_n,m=0ⁱ⁻¹ can be calculated from

V˜nm= 1 b^[i]_n

mi−1

∑

k=0

v_kg^[i]_nmk,

whereg^[i]_nmkis given by g^[i]_nmk=

Z

Ω_i

φn^[i](x_i)φ_m^[i](x_i)φ_k^[i](x_i)w_i(x_i)dx_i. Remark 3. g^[i]_nmk for the shifted Chebyshev and Leg- endre polynomials can be obtained as

i) shifted Chebyshev polynomials defined in[0,Li]:

g^[i]_nmk=L_ic_k

2 δn+m,k+δ|n−m|,k , c0=π, c_k=π

2, k=1,2,· · ·,mi−1. ii) shifted Legendre polynomials defined in[0,L_i]:

g^[i]_nmk=Li

2











dm−ldldn−l

(2n+2m−2l+1)d_n+m−l , k=n+m−2l, l=0,1,· · ·,m,

0, k6=n+m−2l,

l=0,1,· · ·,m, wherem≤nandd_l= (2l)!/2^l(l!)².

The product of two orthogonal vectors will also be applied as

Φ(x,t)Φ(x,t)^TA 'AΦ˜ (x,t),

where ˜Ais anN×Nproduct operational matrix for the vectorA. Using above equation and by the orthogonal property (4), the elements{A˜_{i j}}^N−1_i,_j=0can be calculated from

A˜i j= 1 c_i

N−1

∑

k=0

a_kg_{i jk},

whereg_{i jk}andciare given by g_{i jk}=

Z

Ω×Jφi(x,t)φ_j(x,t)φ_k(x,t)w(x,t)dxdt, g_{i jk}=g^[1]_i

1j1k1 g^[2]_i

2j2k2· · ·g^[d]_i

dj_dk_dg_{r pq}, ci=b^[1]_i

1 b^[2]_i

2 · · ·b^[d]_i

d br,

where i_s,j_s,k_s,r,p, and q are obtained by applying (6) – (8).

3. Numerical Examples

In this section, two examples of homogeneous and inhomogeneous two-dimensional parabolic equations are given to illustrate our results. Also, in the last example, we apply our combined method to the coupled Burgers equations. In all experiments, we consider the shifted Legendre polynomials defined in(0,π)as basis functions. It is noticeable that the results in shifted Chebyshev or other polynomials have no more differ- ence. Comparisons between present results and cor- responding analytical solutions are given. For these comparisons, the root mean square (RMS) error of the following form is applied:

RMS= s

∑^M_k=1 u(x_k,t_k)−u_N(x_k,t_k)2

M ,

whereu(x_k,t_k)andu_N(x_k,t_k)are achieved by the exact and the approximate solution on(x_k,t_k), andMis the number of test points.

3.1. The Inhomogeneous Two-Dimensional Heat Equation

Consider the following equation [12,13]:

f(x₁,x₂,t) =sin(x₁)·sin(x₂)·e^−t−4, g(x₁,x2) =sin(x₁)·sin(x₂) +x²₁+x²₂, 0<x1,x2<π, t>0,

with boundary conditions

u(0,x₂,t) =x²₂, u(x₁,0,t) =x₁², u(π,x₂,t) =x²₂+π², u(x₁,π,t) =x²₁+π². The exact solution of this problem is

(x₁,x2,t) =sin(x₁)·sin(x₂)·e^−t+x²₁+x²₂.

(5)

By applying the mentioned method, (1) is written as A^T

D_t−∇

Φ(x₁,x₂,t) =F^TΦ(x₁,x₂,t), where by using (5), we have f(x,t)'F^TΦ(x₁,x₂,t).

Now similar to the Tau method [18,20], one has EA=F, where E=

D_t−∇T

. (11)

Also the initial and boundary conditions are given as Im₁⊗I_m₂⊗Ψ(0)T

A=G, g(x₁,x₂)'G^TΦ^[1](x₁)⊗Φ^[2](x₂), Φ^[1](0)⊗I_m₂⊗I_sT

A=H₁, x²₂'H₁^TΦ^[2](x₂)⊗Ψ(t), I_m₁⊗Φ^[2](0)⊗I_sT

A=H₂, x²₁'H₂^TΦ^[1](x₁)⊗Ψ(t), Φ^[1](π)⊗I_m₂⊗I_sT

A=H₃, x²₂+π²'H₃^TΦ^[2](x₂)⊗Ψ(t), I_m₁⊗Φ^[2](π)⊗I_sT

A=H₄, x²₁+π²'H₄^TΦ^[1](x₁)⊗Ψ(t).

By substituting above equations in (11), one can obtain a system of algebraic equations and give a unique solution for the unknown coefficients{a_i}^N−1_i=0 . We solved this problem for different N and compared our result with the exact solution in Table1. The graphs of the absolute error functions form₁=7,m₂=7, ands=4 are shown in Figure1. These errors illustrate that the approximate solution is in good agreement with the exact solution.

3.2. The Homogeneous Two-Dimensional Heat Equation

In this example, we consider the two-dimensional linear homogeneous Burgers equation given by [12, 13]

g(x₁,x₂) =sin(x₁)·sin(x₂), 0<x1,x2<π, t>0,

Table 1. RMS error with some values ofN=m₁m₂sin Examples 1, 2, and 3.

Example m₁=m₂=s=4 m₁=m₂=s=5 m₁=m₂=s=6 m₁=m₂=s=7

1 2.85·10⁻³ 1.54·10⁻⁴ 6.73·10⁻⁵ 4.94·10⁻⁷

2 4.57·10⁻³ 3.15·10⁻⁴ 8.61·10⁻⁵ 5.45·10⁻⁷

3 2.67·10⁻³ 1.45·10⁻⁴ 5.36·10⁻⁵ 1.37·10⁻⁷

Fig. 1. Graph of the absolute error withN=196 for Exam- ple 1.

(6)

Fig. 2. Graph of the absolute error withN=196 for Exam- ple 2.

with boundary conditions

u(0,x₂,t) =0, u(x₁,0,t) =0, u(π,x2,t) =0, u(x₁,π,t) =0. The exact solution of this problem is

u(x₁,x2,t) =sin(x₁)·sin(x₂)·e^−2t.

Fig. 3. Graph of the absolute error for the coupled Burgers equation withN=64.

By applying the mentioned method, (1) is written as A^T

D_t−∇

Φ(x₁,x2,t) =0.

(7)

Now, similar to Tau method [18,20], one has EA=0, where E=

D_t−∇T

. (12)

Also the initial and boundary conditions are obtained as

I_m₁⊗I_m₂⊗Ψ(0)T

A=G, g(x₁,x₂)'G^TΦ^[1](x₁)⊗Φ^[2](x₂), Φ^[1](0)⊗I_m₂⊗I_sT

A=0, I_m₁⊗Φ^[2](0)⊗I_sT

A=0, Φ^[1](π)⊗I_m₂⊗I_sT

A=0, Im₁⊗Φ^[2](π)⊗IsT

A=0.

By substituting above equations in (12), one can obtain a system of algebraic equations and give a unique solution for the unknown coefficients{a_i}^N−1_i=0 . We solved this problem for differentN and compared it with the exact solution in Table1. The graphs of the absolute error functions form₁=7,m₂=7, ands=4 are shown in Figure2. These errors demonstrate that the approximate solution obtained by using this scheme is in good agreement with the exact solution.

3.3. The Coupled Burgers Equations

Finally, we consider the following nonlinear system of equations given by [12,13]:

u_t−∇u−2uu_x+ (uv)_x=0, (13) v_t−∇v−2vv_x+ (uv)_x=0, (14) subject to the initial and boundary conditions

u(x,0) =sin(x), v(x,0) =sin(x), u(0,t) =0, v(0,t) =0,

u(π,t) =0, v(π,t) =0. The exact solutions of this system are

u(x,t) =e^−tsin(x), v(x,t) =e^−tsin(x).

By assuming u(x,t) = A^TΦ(x,t) and v(x,t) = B^TΦ(x,t) and applying the proposed method, (13) and (14) are obtained as

A^T

D_t−∇−2D_xA˜+D_xB˜+V˜

Φ(x,t) =0, B^T

D_t−∇−2D_xB˜+D_xA˜+W˜

Φ(x,t) =0, whereV=D^T_xBandW=D^T_xA. Now, similar to the Tau method, we have

D_t−∇−2D_xA˜+D_xB˜+V˜T

A=0, (15) D_t−∇−2D_xB˜+D_xA˜+W˜T

B=0. (16) The conditions of the problem are obtained as

Im₁⊗Ψ(0)T

A=S,

Im₁⊗Ψ(0)T

B=S, Φ^[1](0)⊗I_sT

A=0,

Φ^[1](0)⊗I_sT

B=0, Φ^[1](π)⊗I_sT

A=0, big[Φ^[1](π)⊗I_sT

B=0, where sin(x) =S^TΦ^[1](x). By substituting above equations in (15) and (16), we can obtain a system of nonlinear algebraic equations and have a solution for the unknown coefficients{a_i,b_i}^N−1_i=0 . This problem is solved for different values ofNand the accomplished comparison with the exact solution is shown in Table1.

The graphs of the absolute error functions form₁=8 ands=8 are shown in Figure3. These errors reveal that the approximate solution is in good agreement with the exact solution.

4. Conclusions

In this paper, the Tau method for solving the coupled Burgers equation is applied. The orthogonal functions based on orthogonal polynomials for solving this time-depended equation by the means of the Kronecker product have been constructed. Also, a general formulation for the operational matrices of derivative and product has been derived. The achieved operational matrices along with the Tau method are used to reduce the problem to a system of algebraic equations.

In order to demonstrate the efficiency and reliability of the proposed technique, the root mean square (RMS) error are applied.

[1] J. M. Burgers, Adv. Appl. Mech.1, 171 (1948).

[2] J. D. Cole, Quart. Appl. Math.9, 225 (1951).

[3] S. E. Esipov, Phys. Rev. E.52, 3711 (1995).

[4] J. Nee and J. Duan, Appl. Math. Lett.11, 57 (1998).

[5] J. D. Fletcher, Int. J. Numer. Meth. Fluids 3, 213 (1983).

[6] P. C. Jain and D. N. Holla, Int. J. Numer. Meth. Eng.12, 213 (1978).

(8)

[7] A. R. Bahadir, Appl. Math. Comput.137, 131 (2003).

[8] C. J. Fletcher, J. Comput. Phys.51, 159 (1983).

[9] S. F. Radwan, J. Comput. Appl. Math.174, 383 (2004).

[10] M. A. Abdou and A. A. Soliman, J. Comput. Appl.

Math.181, 245 (2005).

[11] D. Kaya, Int. J. Math. Math. Sci.27, 675 (2001).

[12] M. Dehghan, A. Hamidi, and M. Shakourifar, Appl.

Math. Comput.189, 1034 (2007).

[13] A. Y. A. Kelleci, Numer. Meth. Part. D. E. 27, 982 (2011).

[14] C. Lanczos, J. Math. Phys.17, 123 (1938).

[15] M. Dehghan and A. Saadatmandi, Comput. Math.

Appl.52, 933 (2006).

[16] A. Saadatmandi and M. Dehghan, Commun. Numer.

Meth. Eng.24, 1467 (2008).

[17] T. Tajvidi, M. Razzaghi, and M. Dehghan, Chaos Soli- ton. Fract.35, 59 (2008).

[18] M. Razzaghi, S. Oppenheimer, and F. Ahmad, J. Quant.

Spectrosc. Ra.72, 439 (2002).

[19] A. Saadatmandi and M. Dehghan, Comput. Math.

Appl.59, 1326 (2010).

[20] A. Saadatmandi and M. Dehghan, Comput. Math.

Appl.62, 1135 (2011).

[21] E. H. Doha, A. H. Bhrawy, and S. S. Ezz-Eldien, Appl.

Math. Modell.35, 5662 (2011).