Research Journal of Applied Sciences, Engineering and Technology 4(19): 3858-3864, 2012 ISSN: 2040-7467 © Maxwell Scientific Organization, 2012 Submitted: May 24, Accepted: June 21, 2012 Published: October 01, 2012

A Numerical Approach for Solving a General Nonlinear Wave Equation Zainal Abdul Aziz, Nazeeruddin Yaacob, Mohammadreza Askaripour Lahiji, Mahdi Ghanbari and C.C. Dennis Ling Department of Mathematics, Faculty of Science, Universiti Teknologi Malaysia, 81310 UTM Skudai, Johor, Malaysia Abstract: An analysis of various numerical schemes and boundary conditions on a general nonlinear wave equation is considered in this study. In particular, the Lax-Wendroff, Leapfrog and Iterated Crank Nicholson methods with Dirichlet boundary conditions are used to solve this nonlinear wave equation. The computation of the solution is made via the reduction of the nonlinear wave equation to the two variable and three variable systems. Keywords: Crank-Nicholson method, lax-wendroff method, leapfrog method, solution of nonlinear wave equation, three variable system, two variable system INTRODUCTION There are numerous methods to solve partial differential equations; especially when the concern is to generate numerical solutions (Morton and Mayers, 2005). Applications in numerical partial differential equations are particularly important in gravitationalwave research and numerical relativity (Barry et al., 1999; Kip, 1997; Will, 1999) and also various techniques have been devised over the years to solve such related equations. In this study, we proposed to discuss three numerical methods of relevance. Turkel (1974) studied second order schemes for hyperbolic systems and compared these with respect to stability properties and errors. Mansour et al. (1996) studied solutions using standard (Crank-Nicolson) vectorial beam propagation method for step-index waveguides containing small oscillations. Pelloni and Dougalis (2002) stated a fully discrete spectral method for the numerical solution of the initial and periodic boundary value problems. An alternating Crank-Nicolson method is proposed for the numerical solution of phase-field equations on a dynamically adaptive grid, which automatically leads to two decoupled algebraic subsystems, one is linear and the other is semi linear (Tan and Huang, 2008). Gilles et al. (2000) discussed the stability and accuracy of second order and fourth order accurate spatial central discretizations on staggered grids with a third order accurate spatial discretization on an unstaggered grid, combined with second order Leapfrog time integration scheme. A comparison is made with the standard Lax-Wendroff algorithm and application to obtaining steady-state solutions to the equations of inviscid flow by Kolibal (1992). Wang and Ruuth (2008) studied easily

implementable Variable Step-Size Implicit-Explicit (VSIMEX) linear multistep methods for timedependent PDEs. A relation between the truncation error and exact and approximate amplification factors is derived by Wesseling (1973). In addition, we view the wave equation as a system of first-order equation. There are many different options for this system. In this study, review of the finite difference method has been proposed. Moreover, the three-variable system and the two-variable system are considered. Then it is shown that the nonlinear partial differential equation is solved via the three methods proposed, with Dirichlet boundary conditions. Summary of our approach to designing a finite difference method: We can systematically create finite difference methods by separating the treatment of space and time derivatives, then designing a solver by choosing/testing:    

A time integrator A discretization for spatial derivatives Fourier's analysis for classification of the differential operator Writing code and testing SOME CLASSICAL FINITE DIFFERENCE METHODS

Leapfrog (space and time): We use centered differencing for both space x and time t (Bourchtein and Bourchtein, 2010): u u c t x

(1)

Corresponding Author: Zainal Abdul Aziz, Department of Mathematics, Faculty of Science, Universiti Teknologi Malaysia, 81310 UTM Skudai, Johor, Malaysia

3858

Res. J. Appl. Sci. Eng. Technol., 4(19): 3858-3864, 2012 u mn 1  u mn 1 u n  u mn 1  c ( m 1 ) 2t 2x

where, u  x, t  is the dependent variable, t and

replaced by the discretized centered differences. Consider the Taylor expansion:

(2)

x are

u ( x, t   )  u ( x, t )  

the independent variables representing time and onedimensional space respectively, c is a real positive constant. So we have:

u

n 1 m

u

n 1 m

ct n (um 1  umn 1 )  x

u  2  2 u   O ( 3 ) 2 t 2 t

(9)

We consider the equation:

u u  c t x

(3)

(10)

and It is interesting to note that the leapfrog method is a three-level method. To find the value of the function at one time step, it is important to know the value of the function at the previous two-time steps. Similarly, we know that the leapfrog time stepping method is only stable for operators with eigen values in the range (Hesthaven and Warburton, 2008):

 t   i   1, 1

x c

(5)

(6)

u ( xm , tn 1 )  u ( xm , tn 1 )  u ( xm1 , tn )  u ( xm1 , tn )   ct   2 2x   t  O(t 3 )  c O(x3 ) x

Tmn  O(t 3 )  cO(x3 )

 u c 2 2  2 u  x 2 x 2

(12)

c n c 2 2 n (u i 1  u in1 )  (u i 1  2u in  u in1 ) 2h 2h 2

(13)

This is the approximation of the Eq. (12). This is a second-order difference equation (in space) which differentiates with the second-order difference equation (in time), in series approximation. Crank-Nicholson method: In numerical analysis, the Crank-Nicholson method is a finite difference method used for numerically solving the heat equation and similar partial differential equations. It is a secondorder method, implicit in time and is numerically stable (Teukolsky, 2000; Hansen et al., 2004; Bourchtein and Bourchtein, 2010). The Crank-Nicholson method is based on central difference in space and the trapezoidal rule in time, giving second order convergence in time. Equivalently, it is the average of (7) forward Euler and backward Euler in time. For example, in one dimension, if the partial differential equation is:

Tmn 

x  t so:

u ( x , t   )  u ( x , t )  c

u in 1  u in 

We can perform a full truncation analysis (in space and time):

We know that if

So, to second order:

and thus:

where, m is the number of data , m  1, 2,3 and with a condition :

t 

(11)

(4)

where,  is the eigenvalue and i   1 . The centered difference derivative matrix is a skew symmetric matrix with eigenvalues: ic 2 m s in ( ) x M

 2u  u  2u  (c )  c2 2 t t x x 2

(7)

(8)

Lax-wendroff method: The Lax-Wendroff is a secondorder difference method in both time and space. For this method, we will approximate our function by the Taylor expansions where the time derivatives are

u u  2 u , )  F (u , x , t , t x x 2

(14)

Then, letting u ( i  x , n  t )  u n , the Crank-Nicholson i method is the average of the forward Euler method at n and the backward Euler method at n  1 (Fig. 1) The Crank-Nicholson stencil for a 1D problem is given by: Forward Euler:

3859 

Res. J. Appl. Sci. Eng. Technol., 4(19): 3858-3864, 2012 u  2u  2u  a( 2  2 ) t x y

(22)

can be solved with Crank-Nicholson discretization of: uin, j 1  uin, j  Fig. 1: Crank-Nicholson method

1 a t 2 ( x ) 2

(uin11, j  uin11, j  uin, j 11  uin, j 11  4uin, j 1 )    n n n n n   (ui 1, j  ui 1, j  ui , j 1  ui , j 1  4ui , j )  (23)

u

n 1 i

u t

n i

 Fi n (u , x , t ,

u  u , ) x x 2 2

(15)

Assuming that a square grid is used, so that x  y , then this equation can be simplified by rearranging them and using:

Backward Euler: u in  1  u in u  2u  F i n 1 ( u , x , t , , ) t x x 2

 

(16)

u in 1  u in 1  n 1 u  2 u u  2 u    Fi (u , x , t , , 2 )  Fi n (u , x , t , , 2 )  t t x x x  2 (17)



(1 2)uin,j1  (uin1,1j  uin1,1j  uin,j11  uin,j11)  2



The function F should be discretized spatially with a central difference. Note that this is an implicit method to get the next value of u in time and a system of algebraic equations must be solved. If the partial differential equation is nonlinear, the discretization will be nonlinear so that advancing time will involve the solution of a system of nonlinear algebraic equations.

(1 2)uin, j  (uin1, j  uin1, j  uin, j1  uin, j1) 2 (25)

THREE VARIABLES SYSTEM Proposition: For some functions u  x, t  , v  x, t  , w  x, t  , the partial differential equation u  c 2 u is equivalent tt

u  2u  a t x 2

(18)

As an example, for linear diffusion, whose CrankNicholson discretization is given by: uin 1  uin a (uin11  2uin 1  uin11 )  (uin1  2uin  uin1 )   t 2( x ) 2

at 2( x) 2

(26)

vt  c 2 w x

(27)

wt  v x

 (28)

w (x, 0) = u x (x, 0)

(29)

Proof: First, assume u  c 2u for some function tt xx Set v  ut and w  u x . Then:

(20)

vt  utt  c 2u xx  c 2 wx

 ruin11  (1  2r )uin1  ruin11  ruin1  (1  2r )uin  ruin1 (21) n i

ut  v

where,

(19)

Letting:

xx

to the system of equation:

Example: 1D diffusion:

This is a tridiagonal problem, so that u may be solved by using tridiagonal matrix algorithm. Example: 2D dimensional The two-dimensional heat equation:

(24)

For Crank-Nicholson numerical scheme, thus we can write the scheme as:

Crank-Nicholson:

r 

at ( x)2

u.

(30)

and

wt  u xt  u tx  v x

(31)

Thus, a solution of the wave equation leads to a solution of the three variables system. Now assume 3860 

Res. J. Appl. Sci. Eng. Technol., 4(19): 3858-3864, 2012

wt  v x , vt  c 2 wx and u t  v , with initial condition

ut  vin

w( x,0)  u x ( x,0) . It is easy to see:

u tt  v t  c w x  c u x x 2

2

at time t  0 . Now, we need to show w  u x for all t , That is:

t

v in  1  v in   v t 

 vx  utx . Hence:

0  u tx  u xt 

(38.b)

Considering equations (36.c) and (36.d), the following difference equation is obtained:

 ( w  ux ) t  0

We know w

w in 1  w in 1  Vin 2h

u tt 

(32) 

    

(38.a)

 2

v tt (39)

where,

u x (w  u x ) w   t t t

(33)

Nonlinear system: Let us consider a nonlinear wave equation:

utt  uxx  V (u)  0

w in 1  w in 1  V i n 2h

vt 

v tt 

(34)

where V (u) is some smooth (potential) function. Define three-variable system and write up:

ut  v

(35.a)

vt  wx  V (u )

(35.b)

ux  w

w in  1  w in   w t 

utt  wx  V (u)

(36.b)

v t  w x  V  (u )

(36.c)

 w tt 2

(41)

where, w

Lax-wendroff method: Lax-Wendroff approximates u by the second-order Taylor series, where the derivatives are replaced by 3 equivalent difference equation: (36.a)

(40.b)

Considering equations (36.e) and (36.f), the following difference equation is obtained:

(35.c)

ut  v

v in 1  2 v in  v in 1  v in V i  n h2

(40.a)

t



w tt 

v in 1  v in 1 2 h

w in 1  2 w in  w in 1  w inV i  n h2

(42.a)

(42.b)

Leapfrog method: We will use half-steps for nonlinear wave equation. In fact, we will use the same grid system as in the three-variable wave equation that is, our lattice will be constructed. So if we consider the values of u at the lattice points  i, n  , the difference

vtt  vxx  vV (u )

equation will be:

(36.d)

wt  vx

(36.e)

wtt  wxx  wV (u )

(36.f)

u

n  1 i

 u

n i

 u

 u

n 1 i

 2  ( v in  0 . 5 )

(43.a)

and also the values of v at the lattice points  i, n  5 , the difference equation will be:

Considering equations (36.a) and (36.b), the following difference equation is obtained:

u

n 1 i

t



 2

u

(37) tt

win 0.5  win 0.5 (43.b)  Vi n0.5 ) h and then the values of w at the lattice points  i, n  5  , the difference equation will be: vin  0.5  vin  0.5  2 (

where, 3861 

Res. J. Appl. Sci. Eng. Technol., 4(19): 3858-3864, 2012

w in01.5  w in01.5  2 (

v in 0 .5  v in 0 .5 ) h

(43.c)

Iterative Crank-Nicholson method: We will consider each term in the second term on the right-hand side of the FTCS (Forward Time, Centered Space) method and replaced it by the average of that value and the value of the function at the next time step. Thus, our difference equations are:

u in  1  u v in  1  v in 





n i

 w in11  w in11 2

(

(44.a)

( v in  1  v in )

2

2h

 V in 1 

w in 1  w in 1  V in ) 2h (44.b)

w in  1  w in 

 v in11  v in11 2

(

2



v in 1  v in 1 ) 2

(44.c)

Example: A specific of Eq. (34) can be considered as: V (u )  1  cos( u )

Fig. 2: Crank-Nicholson on sine-gordon equation:

(45)

u ( x,0)  0  u i0  0

which results in the sine-Gordon equation:

utt  u xx  sin(u )

u t ( x ,0 )  0 

(46)

To simulate the periodic boundary condition as follows: u ( 0 , t )  u (1, t )

(47)

(48.a)

u t ( x ,0 )  0

(48.b)

Proposition: For some function, u (x. t) and v(x, t), the partial differential equation u  c 2 u is equivalent to tt

 win11  win11

vin1  vin  ( 2

2h

( v in 1  v in )

 sin(uin1) 

w in 1  w in 

2

(

2

vt  cu x

(51.a)

ut  cv x

(51.b)

(49.a)

Proof: If u  cv and v  cu for functions u , v, then: t x t x

win1  win1  sin(uin )) 2h (49.b)

 v in11  v in11

xx

the system of equations:

We use the Crank-Nicholson method (44) to discretize (46) and to solve the system of equations:

2

(50.b)

TWO-VARIABLE SYSTEM

u ( x ,0 )  0





( u i0  u i1 )  0  u i1  u i0

Note that, the black points are the known data that can be used to evaluate and therefore, the values of ui1 , vi1 , wi1 for i  1, 2 ,..., n can be computed (Fig. 2).

and the initial conditions are given:

u in 1  u in 1 

1

(50.a)

v n  v in1 (49.c)  i 1 ) 2

utt  cvxt  cvtx  c 2uxx

(52)

Now assume u  c 2 u for some function u . There tt xx exists a family of functions v such that vt  cu x . Then each v has the form:

To predict the values two initial conditions are required which they are:

v ( x, t )   cu x dt  f ( x )

3862 

Res. J. Appl. Sci. Eng. Technol., 4(19): 3858-3864, 2012 where, for some function f . Now we choose f so that cv  u at t  0 . Therefore, we x t need only to show that cv  u always, that is x

vt 

uin1  2uin  uin1 Vin 2 h

vtt 

vin1  vin1 n n  vi Vi 2h

t

. Since:  ( cv x  u t )  0 t

 cv x  cv xt  cv tx  c 2 u xx t

Leapfrog method: In the method, we have a second derivative. So we should evaluate both u and v on whole steps, thus our difference equations are: (54)

uin 1  uin 1  2 (vin )

Thus the wave equation and the system are equivalent. Nonlinear system: The nonlinear wave Eq. (34) is considered. It is defined two-variable system and can be written as:

ut  v

(55.a)

vt  uxx  V (u )

(55.b)

Lax-Wendroff method : Lax-Wendrof approximates u by the second-order Taylor series, where the derivatives are replaced by two (56) equivalent difference equation:

ut  v

(56.a)

utt  uxx  V (u)

(56.b) (56.c)

vtt  vxx  vV (u)

(56.d)

2

u tt

u in  1  u in 

 2

( v in  1  v in )

2

(

 Vi  n 1 

u

n i 1

h2  2 u in  u in1  Vi  n ) h2

V (u )  1  cos(u )

 2

v tt

(63)

which results in the sine-Gordon equation: (64)

To simulate the periodic boundary condition as: u ( 0 , t )  u (1, t )

(58.b)

Considering equations (56.c) and (56.d), the following difference equation is obtained:

v in  1  v in   v t 

(62.b)

(57)

(58.a)

u in1  2uin  u in1  Vin 2 h

(62.a)

 u in11  2 u in 1  u in11

utt  u xx  sin(u)

u t  v in

(61.b)

Iterative Crank-Nicholson method: We will consider each term in the second term on the right-hand side of the FTCS method (Forward Time, Centered Space) and replace it by the average of that value and the value of the function at the next time step. Thus, our difference equations are:

where,

u tt 

uin1  2uin  uin1  Vin ) 2 h

(61.a)

Example: A specific of equation (34) can be considered as:

Considering equations (56.a) and (56.b), the following difference equation is obtained:



vin 1  vin 1  2 (

vin 1  vin 

vt  uxx  V (u )

u in  1  u in   u t 

(60.b)

(53)

We have:  (cv x  u t )  c 2 u xx  u tt  0 t

(60.a)

(65)

and initial conditions are given:

(59)

3863 

u ( x ,0 )  0

(66)

u t ( x ,0 )  v

(67)

Res. J. Appl. Sci. Eng. Technol., 4(19): 3858-3864, 2012

Fig. 3: Leapfrog for sine-gordon

We use the Leapfrog method (62) to discretize (65), so solve the system of equations:

u in 1  u in 1  2 ( vin ) vin 1  vin 1  2 (

u in1  2u in  u in1  sin( u in )) h2

(68.a)

(68.b)

Note that, the black points are the known date that can be used to evaluate. Therefore, the values of u1 , v1 i i for i  1,2,..., n can be computed (Fig. 3). CONCLUSION In this study, we have used Lax-Wendroff, Leapfrog and Iterated-Crank-Nicholson methods to solve numerically certain generalized nonlinear wave equation in terms of two-variable and three-variable systems. In addition, provided with the knowledge of using the various numerical methods and imposing the various boundary conditions, we consider a nonlinear system without general solution. REFERENCES Barry, C., W. Barish and W. Rainer, 1999. LIGO and the detection of gravitational waves. Phys. Today, 52: 44-50. Bourchtein, A. and L. Bourchtein, 2010. On iterated Crank-Nicholson methods for hyperbolic and parabolic equations. Comp. Phys. Commun., 181: 1242-1250. Gilles, L., S.C. Hangness and L. Vazques, 2000. Comparison between staggered and unstaggered finite difference time domain grids for few-cycle temporal optical soliton propagation. J. Comput. Phys., 161: 379-400.

Hansen, J., A. Khokhlov and I. Novikov, 2004. Properties of four numerical schemes applied to a scalar nonlinear scalar wave equation with a GRType nonlinearity. Int. J. Mod. Phys. D., 13: 961-982. Hesthaven, J.S. and T. Warburton, 2008. Nodal discontinuous Galerkin Methods: Algorithms, analysis and application. Springer Texts in Applied Mathematics. Springer Verlag, New York, Vol. 54. Kip, S.T., 1997. Probing Black Holes and Relativistic Stars with Gravitational Waves. Retrieved from: arXiv: gr-qc/9706079. Kolibal, J., 1992. Extensions of Lax-Wendroff to multicolor schemes. Ulam Quaterly., 1: 12-29. Mansour, I., A. Daniele and C. Rosa, 1996. Noniterative vectorial beam propagation method with a smoothing digital filter. J. Lightwave Technol., 4: 908-913. Morton, K.W. and D.F. Mayers, 2005. Numerical Solution of Partial Differential Equation. Cambridge University Press, United Kingdom, pp: 102-140. Pelloni, B. and V.A. Dougalis, 2002. Error estimates for a fully discrete spectral scheme for a class of nonlinear, nonlocal dispersive wave equation. Appl. Numer. Math., 37: 95-107. Tan, Z. and Y. Huang, 2008. An alternating CrankNicolson method the numerical solution of the phase-field equations using adaptive moving meshes. Int. J. Numer. Meth. Fluids., 56: 1673-1693. Teukolsky, S.A., 2000. Stability of the iterated CrankNicholson method in numerical relativity. Phys. Rev. D., 61: 087501. Turkel, E., 1974. Phase error and stability of second order methods for hyperbolic problems. J. Comput. Phys., 15: 226-250. Wang, D. and S.J. Ruuth, 2008. Variable step-size implicit-explicit linear multistep methods for timedependent partial differential equations. J. Comput. Math., 26: 838-855. Wesseling, P., 1973. On the construction of accurate difference schemes for hyperbolic partial differential equations. J. Eng. Math., 7: 19-31. Will, G.M., 1999. Gravitational radiation and the validity of general relativity. Phys. Today., 52: 485-495.

3864