About this Journal Submit a Manuscript Table of Contents
Abstract and Applied Analysis
Volume 2014 (2014), Article ID 323260, 10 pages
http://dx.doi.org/10.1155/2014/323260
Research Article

A Fourth Order Finite Difference Method for the Good Boussinesq Equation

Department of Mathematics, College of Science, King Abdulaziz University, P.O. Box 80203, Jeddah 21589, Saudi Arabia

Received 24 October 2013; Accepted 16 January 2014; Published 26 February 2014

Academic Editor: Youyu Wang

Copyright © 2014 M. S. Ismail and Farida Mosally. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The “good” Boussinesq equation is transformed into a first order differential system. A fourth order finite difference scheme is derived for this system. The resulting scheme is analyzed for accuracy and stability. Newton’s method and linearization techniques are used to solve the resulting nonlinear system. The exact solution and the conserved quantity are used to assess the accuracy and the efficiency of the derived method. Head-on and overtaking interactions of two solitons are also considered. The numerical results reveal the good performance of the derived method.

1. Introduction

In recent years, remarkable developments have taken place in the study of nonlinear evolutionary partial differential equations. It is realized that many such equations possess special solutions in the form of pulses which retain their shapes and velocities after interacting with each other. Such solutions are called solitons. Many equations admitting soliton solutions are as follows: sine Gordon and double sine Gordon equations, Schrodinger equation, and KdV, MKdV, and complex modified KdV equations; many research works have been done on these equations [18]. Most of the current research is directed to solve coupled nonlinear systems analytically and numerically [921]. Solitons are of great interest in many physical areas, as, for example, in dislocation theory of crystals, plasma and fluid dynamics, magnetohydrodynamics, laser and fiber optics, and the study of the water waves.

In this work we will study numerically the “good” Boussinesq (GB) nonlinear equation which provides a balance between dispersion and nonlinearity, that leads to the existence of soliton solutions, similar to the Korteweg-de Vries (KdV) and cubic nonlinear Schrödinger equation [1, 8, 22].

The initial displacement associated with the partial differential equation given in (1) is assumed to take the form [2325] and the boundary conditions

The solitary wave solution of (1) is given by [3, 14] where is amplitude of the pulse and is an arbitrary parameter. In addition, the double soliton solution [3, 14, 24, 26] can be written as where is the velocity of the th soliton, is the corresponding amplitude, is an arbitrary parameter, and is the initial position. The exact solution (5) represents two solitary waves located initially at the positions and moving to the right or left according to their velocities.

Many research works on Boussinesq equation have been developed. Analytical solution of this equation was studied by many authors, such as the construction of -soliton solutions using the bilinear form [4], multiple soliton solutions for the GB equation using a simplified version of Hirota method [27] and Adomian decomposition method [28]. Construction of soliton solutions and periodic solution of Boussinesq equation by modified decomposition method are given in [29, 30]. A variational iteration method is developed for GB equation [31]. A solitary wave solution of the Boussinesq equation with power law nonlinearity is derived in [32].

Many numerical methods have been developed for solving the Boussinesq equation, such as Petrov-Galerkin method [19]. Bratsos et al. [23, 24, 33] have developed finite difference schemes and have considered the Bad Boussinesq (BB) and GB equations. El-Zhoheiry has developed an implicit finite difference scheme [34]. Aydin and Karasözen [35] constructed second order symplectic and multisymplectic integrators for the GB equation using the two-stage Lobatto IIIA-IIIB partitioned Runge-Kutta method. Daripa and Hua [36] developed finite difference method for BB equation of second order in time and space, where they have used the filtration and regularization techniques to control the growth of the errors arising from the instabilities and provide better approximate solutions of this equation. Ismail and Bratsos [25] have derived a predictor corrector scheme for the GB and BB equations; the scheme is fourth order in time and second order in space, and it is conditionally stable. Matsuo [37] also has derived conservative finite difference schemes for certain classes of nonlinear wave equations with some applications for the nonlinear Klein-Gordon and GB equations. Mohebbi and Asgari [26] also have solved the GB equation using a fourth order time stepping schemes with combination of discrete Fourier transform. Split step Fourier method is also used to solve Boussinesq-type equations. Dehghan and Salehi [38] have derived a mesh free method for the classical Boussinesq equation.

The paper is laid out as follows. In Section 2, the GB equation is transformed into a first order differential system in time; a finite difference scheme is derived for this system. The method is analyzed for accuracy and stability. In Section 3, a linearization technique is used to convert the nonlinear system obtained into a linear block tridiagonal system. Numerical tests are presented in Section 4. Concluding remarks are given in Section 5.

2. The Numerical Method

In order to derive a highly accurate method, we transform the GB equation (1) into the first order differential system in time as [19] with initial conditions and boundary conditions

By using the boundary conditions (10), the system (7) and (8) has the conserved quantity [2426, 38]

In order to derive a numerical method for solving the system (7) and (8), the region with its boundary consisting of the lines , and the axis is covered with a rectangular mesh of points with coordinates [25] where and are the space and time increments, respectively. We denote the theoretical solution of (7) and (8) at the typical mesh point by and , and the numerical solution by and . The fourth finite difference formula [39] is used to approximate the second order space derivative which appears in our system. By making use of (13) and the implicit midpoint rule, the finite difference nonlinear implicit scheme is obtained. The schemes (15) and (16) can be written as where . On expansion of the central differences using (14), this will lead us to the difference scheme By imposing the boundary conditions, a coupled nonlinear system is obtained, which can be solved by any iterative scheme. This system can be solved by Newton's method or fixed point method. Linearization method can also be used.

2.1. Accuracy of the Scheme

To study the accuracy of the proposed schemes (19) and (20), the numerical solutions and are replaced by the exact solutions and , respectively. By making use of the following expansions about the point : into (19) to obtain

By making use of the differential system (7) and (8), all terms inside the brackets in (22) are zero. Similar analysis can be done for (20). This will lead us to the conclusion that the derived scheme is of second order in time and fourth order in space, that is, . The numerical results confirm this.

2.2. Stability of the Scheme

To study the stability of the derived scheme, von Neumann stability analysis will be used. We consider the linear version of the proposed schemes (17) and (18) which can be given as where is a constant quantity defined by

To apply von Neumann method, we assume that the solutions of (22) and (23) are where , , and are constants. By substituting these solutions into (23) and (24), and after some manipulation, we get the system where . Equations (27) can be written in a matrix vector form as where

The von Neumann necessary condition for the stability of this system is the maximum eigenvalue of value of the amplification matrix in (28) is to be less than or equal to one. By direct calculation, the eigenvalues of the matrix are with modulus equal to one, and hence the proposed scheme is unconditionally stable.

3. Linearization Technique

In order to avoid the solution of the nonlinear systems (19) and (20), a linearization technique will be developed to overcome this difficulty. Using Taylor's series expansion of the nonlinear term about , we obtain and this will lead us to the following approximation: which preserves second order temporal accuracy. By replacing the exact solution by the approximating one, we get

By substituting (33) into (18), we obtain the linearly implicit finite difference scheme

By expanding the central difference operator in (34), we obtain the linearized scheme

By utilizing the boundary conditions, the finite difference scheme (35) produces a linear block tridiagonal system in the unknowns , which can be solved directly by Crout's method. The obtained difference scheme is still of second order in time and fourth order method in space, and it is unconditionally stable.

4. Numerical Results

In all experiments we use the following values: , , , , and . We study the accuracy of the proposed method by calculating the infinity error norm The numerical solution and trapezoidal rule are used to calculate the error and the conserved quantity [19, 25, 26, 37]. All numerical results in this section are obtained from the solution of the nonlinear schemes (19) and (20) using Newton's method. Linearization method is used for comparison purpose only.

4.1. Single Soliton

To test the derived method, we consider the initial condition where (37) represent soliton and kink solution, respectively. The parameters , , and are used in this test. The numerical results for the nonlinear and linearized schemes are given in Tables 1 and 2, respectively. The numerical results (amplitude and the conserved quantity) obtained display the high accuracy of the proposed method. The execution time required for the nonlinear scheme is  sec. compared to  sec. for the linearized scheme. The numerical results produced in Table 1 are more accurate than the one in Table 2, and this is due to the approximation of the nonlinear term in the linearization process. A comparison of some existing methods is given in Table 3, which indicates that our method is the most accurate one. In Figures 1 and 2, we display the numerical solutions and moving to the right for .

tab1
Table 1: Nonlinear scheme of single soliton .
tab2
Table 2: Linearization method of single soliton with .
tab3
Table 3: Comparison of numerical solutions with at .
323260.fig.001
Figure 1: Single soliton: numerical solution of .
323260.fig.002
Figure 2: Single kink: numerical solution of .

To test whether the proposed numerical scheme exhibits the expected convergence rate in space, we perform some numerical experiments for various values of and fixed value of . In these experiments we choose to ensure that the temporal error is negligible. The rate of convergence for the scheme is calculated using the formula where is the infinity error norm. The errors and the rate of convergence for are given in Table 4.

tab4
Table 4: Rate of convergence with .
4.2. Head-On Soliton Interaction

To study the head-on collision of two solitons, we choose the initial conditions where

In [19] it was reported that solution blew up numerically when . By choosing the parameters , , and . The initial conditions represent two solitons and two kinks located at and . In Figures 3 and 4, we display the interaction scenario for to ; we can easily observe that the two solitons and the two kinks have been separated completely after the interaction and restored their original shapes and velocities. The velocity and the conserved quantity are given in Table 5. It is noted that when the two waves overlap (Figures 3 and 4), the joint amplitude is greater than the sum of individual amplitudes; this is in full agreement with [37]. By choosing , the situation dramatically changes; both Figures 5 and 6 end around ~. Newton's method fails to find the solution there. This agrees with the blowup results in [19, 37].

tab5
Table 5: Interaction of two solitons .
323260.fig.003
Figure 3: Head-on interaction: numerical solution of   with .
323260.fig.004
Figure 4: Head-on interaction: numerical solution of   with .
323260.fig.005
Figure 5: Blowup: numerical solution of with .
323260.fig.006
Figure 6: Blowup: numerical solution of with .
4.3. Overtaking Soliton Interaction

In this test we will study the interaction of two solitons moving in the same direction. In this case we choose the initial conditions where

The following parameters are selected:

The initial conditions represent two solitons and two kinks are initially located at and , moving in the same direction with velocities and . The usual interaction has taken place and the faster wave interacted and separated from the slower one and left it behind; the interaction scenario is given in Figures 7 and 8. In Table 6, we track the amplitude and the conserved quantity during the interaction scenario. In Figure 9, we display the contours of the numerical solution . Many numerical tests have been conducted for different amplitudes. We have noticed that the interaction occurred for , , which is completely different from the head-on interaction case, where the solution blew up numerically when , . This type of interaction for the GB equation seems to be reported in the literature for the first time as far as the authors know. We have noticed that when the two waves overlap (Figures 7 and 8), the joint amplitude is less than the sum of individual amplitudes.

tab6
Table 6: Overtaking interaction with .
323260.fig.007
Figure 7: Overtaking soliton interaction: numerical solution of with , .
323260.fig.008
Figure 8: Overtaking soliton interaction: numerical solution of with , .
323260.fig.009
Figure 9: Contours of overtaking soliton interaction, where , .
4.4. Birth of Solitons

To study the evolution of arbitrary pulse [19], we choose the initial conditions

This test is carried out for values of , mainly in the vicinity of , the theoretical value of the amplitude for . We considered two cases and . In the first case we choose ; we have noticed that as time evolves, the initial pulse split into two solitons moving in the opposite directions with equal amplitudes and velocities ; see Table 7. In Figures 10 and 11, we display the splitting scenario and some dispersive oscillations which occurred after splitting. In the second case we choose ; as time evolves, the initial pulse is split into two solitons moving in the opposite direction with amplitudes ; see Table 8. Figures 12 and 13 display the splitting scenario; no dispersive oscillation is observed in this case (linearized scheme failed to produce the numerical solution in this case). Many numerical tests have been conducted for different amplitudes; we have noticed that the numerical solution blew up if .

tab7
Table 7: Birth of solitons with .
tab8
Table 8: Birth of solitons with ,  .
323260.fig.0010
Figure 10: Birth of solitons, where .
323260.fig.0011
Figure 11: Birth of solitons, where .
323260.fig.0012
Figure 12: Birth of solitons, where .
323260.fig.0013
Figure 13: Birth of solitons, where .

5. Concluding Remarks

In this paper, a highly accurate finite difference method to solve GB equation is proposed. The main idea is to transform the GB equation into a first order differential system in time. The method is fourth order in space and second order in time, and it is unconditionally stable. The method produced a coupled nonlinear system; Newton's and linearization methods are used to solve it. The nonlinear scheme produced more accurate results than the linearized scheme. Numerical results for single soliton, head-on, and overtaking solitons interactions are given; the most interesting phenomena we observed were when two solitons collide; if both solitons are small enough, they pass through each other like the other usual solitons do, but when they exceed some limit, the solution blows up at the collision, even if both amplitudes are smaller than for being stable solitons [14, 15]. The overtaking solitons interaction is also discussed as a new numerical test. Elastic interaction occurred for large amplitudes (, , and this is contrary to the head-on interaction, where blowup is observed for moderate values of the amplitudes (, . Birth of solitons is also observed for initial pulse with amplitude . A blowup was observed when and this is in full agreement with the existing methods [14, 15].

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors would like to express their sincere thanks to the referees for their valuable comments on the improvement of the paper.

References

  1. M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform, vol. 4 of SIAM Studies in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 1981. View at MathSciNet
  2. A. Biswas and K. R. Raslan, “Numerical simulation of the modified Korteweg-de Vries equation,” Physics of Wave Phenomena, vol. 19, no. 2, pp. 142–147, 2011. View at Publisher · View at Google Scholar · View at Scopus
  3. G. Domairry, A. G. Davodi, and A. G. Davodi, “Solutions for the double sine-Gordon equation by exp-function, tanh, and extended tanh methods,” Numerical Methods for Partial Differential Equations, vol. 26, no. 2, pp. 384–398, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  4. R. Hirota, “Exact N-soliton solutions of the wave equation of long waves in shallow-water and in nonlinear lattices,” Journal of Mathematical Physics, vol. 14, pp. 810–814, 1973. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  5. H. Leblond, H. Triki, and D. Mihalache, “Derivation of a generalized double sine-Gordon equation describing ultrashort soliton propagation in optical media composed of multilevel atoms,” Physical Review A, vol. 86, Article ID 063825, 9 pages, 2012. View at Publisher · View at Google Scholar
  6. H. Liu and J. Li, “Symmetry reductions, dynamical behavior and exact explicit solutions to the Gordon types of equations,” Journal of Computational and Applied Mathematics, vol. 257, pp. 144–156, 2014. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  7. R. Mokhtari and M. Mohseni, “A meshless method for solving mKdV equation,” Computer Physics Communications, vol. 183, no. 6, pp. 1259–1268, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  8. G. B. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, New York, NY, USA, 1974. View at MathSciNet
  9. A. Aydın and B. Karasözen, “Lobatto IIIA-IIIB discretization of the strongly coupled nonlinear Schrödinger equation,” Journal of Computational and Applied Mathematics, vol. 235, no. 16, pp. 4770–4779, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  10. D. Bai and J. Wang, “The time-splitting Fourier spectral method for the coupled Schrödinger-Boussinesq equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 3, pp. 1201–1210, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  11. M. S. Ismail and T. R. Taha, “Parallel methods and higher dimensional NLS equations,” Abstract and Applied Analysis, vol. 2013, Article ID 497439, 9 pages, 2013. View at Publisher · View at Google Scholar · View at MathSciNet
  12. M. S. Ismail and A. Biswas, “1-soliton solution of the Klein-Gordon-Zakharov equation with power law nonlinearity,” Applied Mathematics and Computation, vol. 217, no. 8, pp. 4186–4196, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  13. M. S. Ismail, “Numerical solution of complex modified Korteweg-de Vries equation by collocation method,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 3, pp. 749–759, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  14. M. S. Ismail, “Numerical solution of a coupled Korteweg-de Vries equations by collocation method,” Numerical Methods for Partial Differential Equations, vol. 25, no. 2, pp. 275–291, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  15. M. S. Ismail, “A fourth-order explicit schemes for the coupled nonlinear Schrödinger equation,” Applied Mathematics and Computation, vol. 196, no. 1, pp. 273–284, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  16. M. S. Ismail, “Numerical solution of coupled nonlinear Schrödinger equation by Galerkin method,” Mathematics and Computers in Simulation, vol. 78, no. 4, pp. 532–547, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  17. M. S. Ismail and T. R. Taha, “A linearly implicit conservative scheme for the coupled nonlinear Schrödinger equation,” Mathematics and Computers in Simulation, vol. 74, no. 4-5, pp. 302–311, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  18. S. Kutluay and Y. Ucar, “Numerical solution of coupled modied Korteweg-devries equation by Galerkin method using quadratic spline,” International Journal of Computer Mathematics, vol. 90, no. 11, pp. 2353–2371, 2013. View at Publisher · View at Google Scholar
  19. V. S. Manoranjan, A. R. Mitchell, and J. Ll. Morris, “Numerical solutions of the good Boussinesq equation,” Society for Industrial and Applied Mathematics, vol. 5, no. 4, pp. 946–957, 1984. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  20. S. S. Ray, “Soliton solutions for time fractional coupled modified KdV equations using new coupled fractional reduced differential transform method,” Journal of Mathematical Chemistry, vol. 51, no. 8, pp. 2214–2229, 2013. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  21. S. Zhou and X. Cheng, “Numerical solution to coupled nonlinear Schrödinger equations on unbounded domains,” Mathematics and Computers in Simulation, vol. 80, no. 12, pp. 2362–2373, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  22. H. Hu, L. Liu, and L. Zhang, “New analytical positon, negaton and complexiton solutions of a coupled KdV-mKdV system,” Applied Mathematics and Computation, vol. 219, no. 11, pp. 5743–5749, 2013. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  23. A. G. Bratsos, “A second order numerical scheme for the solution of the one-dimensional Boussinesq equation,” Numerical Algorithms, vol. 46, no. 1, pp. 45–58, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  24. A. G. Bratsos, Ch. Tsitouras, and D. G. Natsis, “Linearized numerical schemes for the Boussinesq equation,” Applied Numerical Analysis and Computational Mathematics, vol. 2, no. 1, pp. 34–53, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  25. M. S. Ismail and A. G. Bratsos, “A predictor-corrector scheme for the numerical solution of the Boussinesq equation,” Journal of Applied Mathematics & Computing, vol. 13, no. 1-2, pp. 11–27, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  26. A. Mohebbi and Z. Asgari, “Efficient numerical algorithms for the solution of “good” Boussinesq equation in water wave propagation,” Computer Physics Communications, vol. 182, no. 12, pp. 2464–2470, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  27. A.-M. Wazwaz, “Multiple-soliton solutions for the Boussinesq equation,” Applied Mathematics and Computation, vol. 192, no. 2, pp. 479–486, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  28. B. S. Attili, “The Adomian decomposition method for solving the Boussinesq equation arising in water wave propagation,” Numerical Methods for Partial Differential Equations, vol. 22, no. 6, pp. 1337–1347, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  29. A.-M. Wazwaz, “New travelling wave solutions to the Boussinesq and the Klein-Gordon equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 13, no. 5, pp. 889–901, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  30. A. M. Wazwaz, “Construction of soliton solutions and periodic solutions of the Boussinesq equation by the modified decomposition method,” Chaos, Solitons and Fractals, vol. 12, no. 8, pp. 1549–1556, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  31. A. Yildirim and S. T. Mohyud-Din, “A variational approach for soliton solutions of good Boussinesq equation,” Journal of King Saud University, vol. 22, no. 4, pp. 205–208, 2010. View at Publisher · View at Google Scholar · View at Scopus
  32. A. Biswas, D. Milovic, and A. Ranasinghe, “Solitary waves of Boussinesq equation in a power law media,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 11, pp. 3738–3742, 2009. View at Publisher · View at Google Scholar · View at Scopus
  33. A. G. Bratsos, “A parametric scheme for the numerical solution of the Boussinesq equation,” The Korean Journal of Computational & Applied Mathematics, vol. 8, no. 1, pp. 45–57, 2001. View at Zentralblatt MATH · View at MathSciNet
  34. H. El-Zoheiry, “Numerical investigation for the solitary waves interaction of the “good” Boussinesq equation,” Applied Numerical Mathematics, vol. 45, no. 2-3, pp. 161–173, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  35. A. Aydin and B. Karasözen, “Symplectic and multisymplectic Lobatto methods for the “good” Boussinesq equation,” Journal of Mathematical Physics, vol. 49, no. 8, Article ID 083509, 2008. View at Publisher · View at Google Scholar · View at Scopus
  36. P. Daripa and W. Hua, “A numerical study of an ill-posed Boussinesq equation arising in water waves and nonlinear lattices: filtering and regularization techniques,” Applied Mathematics and Computation, vol. 101, no. 2-3, pp. 159–207, 1999. View at Publisher · View at Google Scholar · View at MathSciNet
  37. T. Matsuo, “New conservative schemes with discrete variational derivatives for nonlinear wave equations,” Journal of Computational and Applied Mathematics, vol. 203, no. 1, pp. 32–56, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  38. M. Dehghan and R. Salehi, “A meshless based numerical technique for traveling solitary wave solution of Boussinesq equation,” Applied Mathematical Modelling, vol. 36, no. 5, pp. 1939–1956, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  39. A. R. Mitchell and D. Gritthes, Finite Difference Method, Wiley, New York, NY, USA, 1980.