A Decomposition Method for Solving the Nonlinear Klein–Gordon Equation

JOURNAL OF COMPUTATIONAL PHYSICS ARTICLE NO.

124, 442–448 (1996)

0071

A Decomposition Method for Solving the Nonlinear Klein–Gordon Equation E. Y. DEEBA

AND

S. A. KHURI

Department of Computer and Mathematical Sciences, University of Houston-Downtown, Houston, Texas 77002 Received May 15, 1995; revised September 6, 1995

and b is real, g is a given nonlinear function, and f is a known function. Equations (1.1)–(1.2) is one of the important mathematical models in quantum mechanics [23–24] and it occurs in relativistic physics [20–21] as a model of dispersive phenomena. There are numerous papers dealing with the existence, uniqueness of the smooth and weak solutions of (1.1)–(1.2) (see [18]), and with the numerical solutions using finite difference, finite element, or collocation methods such as in [17, 19, 22].

In this paper, Adomian’s decomposition scheme is presented as an alternate method for solving the nonlinear Klein–Gordon equation. The method is demonstrated by several examples. Comparing the scheme with existing collocation, finite difference and finite element techniques shows that the present approach is highly accurate and converges rapidly. Q 1996 Academic Press, Inc.

1. INTRODUCTION

Recently, Adomian’s decomposition scheme is emerging as an alternate method for solving a wide range of problems whose mathematical models yield equations or system of equations involving algebraic, differential, integral, integro-differential, or differential-delay terms (for example, see [1–16]). It has been shown that the method yields rapidly convergent series solutions to linear and nonlinear deterministic and stochastic equations. In some instances, if the problem is nonlinear, then the linearization of the perturbation method yields a set of partial differential equations that needs to be solved. On the other hand, the decomposition method provides a direct scheme for solving the underlying problem and does not require linearization. We shall, among other things, present this in the paper. Our attention will focus on the nonlinear Klein–Gordon equation of the form

2. DECOMPOSITION METHOD AND KLEIN–GORDON EQUATION

In this section we shall describe the main algorithm of Adomian’s decomposition method as it applies to a general nonlinear equation of the form u 2 N(u) 5 f,

where N is a nonlinear operator on a Hilbert space H, f is a known element of H, and we are seeking u [ H satisfying (2.1). We assume that for every f [ H, Eq. (2.1) has a unique solution. The decomposition technique consists of representing the solution of (2.1) as a series

­ 2u(x, t) 2 Du(x, t) 1 bu(x, t) 1 g(u(x, t)) 5 f(x, t) (1.1) ­t 2 u(x, 0) 5 a0 (x),

­u (x, 0) 5 a1 (x) ­t

(2.1)

u5

(1.2)

Ou y

(2.2)

n

n50

and the nonlinear operator N is decomposed as

with

N(u) 5

OA , y

(2.3)

n

n50

x 5 (x1 , x2 , ..., xm) [ Rm, t [ (0, T ], where the An’s are Adomian’s polynomials of u0 , ..., un given by

where D5

O ­x­ m

j50

2 2 j

An 5

(1.3) 442

0021-9991/96 $18.00 Copyright  1996 by Academic Press, Inc. All rights of reproduction in any form reserved.

F SO DG

1 dn N n! dln

y

i50

liui

, n 5 0, 1, ... .

l5 0

(2.4)

443

SOLVING THE NONLINEAR KLEIN–GORDON EQUATION

Upon substituting Eq. (2.2) and (2.3) into the functional equation (2.1) yields

O u 2 O A 5 f. y

The An9 s can finally be written in the following convenient way

y

n

(2.5)

n

n50

An 5

n50

O c(n, n)g n

(n)

n5 1

(u0).

(2.8)

The result is The convergence of the series in (2.5) will yield A0 5 g(u0) u0 5 f

A1 5 u1

u1 5 A0 u2 5 A1

(2.6)

A3 5 u3

_ Thus, one can recurrently determine every term of the y series on50 un . The convergence of this series has been established (see [13, 14]). In [14] a proof of convergence is established using fixed point theorems. In [13] the hypotheses for proving convergence are less restrictive and are generally satisfied in physical problems. The two hypotheses that are necessary for proving convergence of Adomian’s technique as given in [13] are:

To illustrate the scheme, let N(u) be a nonlinear function of u, say g(u), where

u31 d3 d d2 g(u0) 1 u1u2 2 g(u0) 1 g(u0). du0 du0 3! du30

How do we interpret and solve the Klein–Gordon equation in this setting? Following the Adomian decomposition analysis [4-5], define the linear operators Lt 5

­2 ­2 Lxi 5 2 , i 5 1, 2, ..., m. 2, ­t ­xi

(2.10)

Consequently, (1.1) can be rewritten in terms of the operators (2.10) as

1. The nonlinear functional equation (2.1) has a series y y solution on50 un such that on50 (1 1 «)nuunu , y, where « . 0 may be very small. 2. The nonlinear operator N(u) can be developed in y series according to u: N(u) 5 on50 anun.

(2.9)

u21 d 2 d A2 5 u2 g(u0) 1 g(u0) du0 2! du20

_ un 5 An21

d g(u0) du0

Ltu 5

O L u 2 bu 2 g(u) 1 f. m

i51

(2.11)

xi

It was shown in [18] that Eq. (1.1) with conditions (1.2) possesses a unique solution. Thus the inverse operator of Lt , namely Lt21, exists and is the twofold indefinite integral; i.e., [Lt21 f ](t) :5

u 5 u0 1 lu1 1 l u2 1 ... 2

E du E dv f (v). t

u

0

0

(2.12)

Operating on both sides of (2.11) with Lt21 yields then the first four Adomian’s polynomials An are given by Lt21Lt u 5 A0 5 g(u( l ))ul50 5 g(u0) A1 5 (dg/du)(du/dl )ul50 A2 5 As[(d h/du )(du/dl ) 1 (dg/du)(d u/dl )]ul50 2

2

2

2

2

(2.7)

A3 5 Ah[(d 3g/du3)(du/dl )3 1 2(d 2g/du2)(du/dl )(d 2u/dl2) 1 (d 2g/du2)(d 2u/dl2)(du/dl ) 1 (dg/du)(d 3u/dl3)]ul50 .

OL m

i50

21 t Lxiu

2 bLt21u 2 Lt21(g(u)) 1 Lt21f (2.13)

from which it follows, upon using the initial conditions given in (1.2), u(x, t) 5 a0(x) 1 a1(x)t 1

OL m

i50

21 t Lxiu

2 bLt u 2 Lt (g(u)) 1 Lt f. 21

21

21

(2.14)

444

DEEBA AND KHURI

The Adomian decomposition method yields the solution in the form given in (2.2). The first term u0 can be determined as u0(x, t) 5 a0(x) 1 a1(x)t 1 Lt ( f (x, t)). 21

y

un11 5

OL

t

1 Lt21(22 sin x sin t).

21

y

u0 5 t sin x 1 Lt21(22 sin x sin t)

21

Hence, using (2.9), the first three iterates are given by

5 2 sin x sin t 2 t sin x

m

u1 5 Lt21Lxu0 1 2Lt21u0 5 Lt21(Lx(2 sin x sin t 2 t sin x))

u2 5 oi50 Lt21Lxiu1 2 bLt21u1 2 Lt21(u1g(u0)) m

S

2 Lt21 u2

1 2Lt21(2 sin x sin t 2 t sin x)

(2.17)

m

D

5 22 sin t sin x 1 2t sin x 2

For later numerical computation, let the expression

O u (x, t)

1 3 t sin x. 3!

In a like manner, we find

n21

i

(3.4)

5 Lt21(2 sin x sin t 2 t sin x)

u21 d2 d g(u0) 1 g(u0) . du0 2! du20

fn 5

(3.3)

and

u1 5 oi50 Lt21Lxiu0 2 bLt21u0 2 Lt21(g(u0))

u3 5 oi50 Lt21Lxiu2 2 bLt21u2

(3.2)

If u(x, t) 5 on50 un , then (2.15) and (2.17) imply that the various iterates can be determined as

Lxiun 2 bLt un 2 Lt An , n $ 0. (2.16)

21

i50

u(x, t) 5 t sin x 1 Lt21Lxu 1 2Lt21u

(2.15)

If we set N(u) 5 g(u) 5 on50 An , then the next iterates are determined via the recursive relation m

Equation (2.14) implies that

u2 5 Lt21Lxu1 1 2Lt21u1

(2.18)

i 50

t5 5 2 sin x sin t 2 2t sin x 1 Ad t sin x 2 sin x 5!

(3.5)

3

denote the n-term approximation to u. Since the series n21 converges very rapidly, the sum fn 5 oi50 ui can serve as a practical solution. We will show through several examples, that the number of terms required to obtain an accurate computable solution is very small.

and u3 5 Lt21Lxu2 1 2Lt21u2 5 22 sin x sin t 1 2t sin x 2 Ad t 3 sin x

3. ILLUSTRATION OF THE METHOD

1

In this section we shall consider four examples where in (1.1) g is assumed to be either u, u2, u3, or sin u, respectively. The outcome of Adomian’s decomposition method shall be compared with any known solution to the underlying Klein–Gordon equation. The solutions obtained are generated by using MAPLE. EXAMPLE 1. Consider the Klein–Gordon equation of the form utt 2 uxx 2 2u 5 22 sin x sin t u(x, 0) 5 0, ut(x, 0) 5 sin x.

(3.1)

The term 2 sin x sin t will be shown to be a noise term.

(3.6)

7

2 5 t t sin x 2 sin x. 5! 7!

Also, u4 5 Lt21Lxu3 1 2Lt21u3 5 2 sin x sin t 2 2t sin x 1 Ad t 3 sin x 2

(3.7)

2 1 2 5 t sin x 1 t7 sin x 2 t 9 sin x. 5! 7! 9!

Upon summing these iterates, we observe that

f4 5

O u 5 sin x St 2 3!t 1 5!t 2 7!t D 3

3

i

i50

5

7

(3.8)

445

SOLVING THE NONLINEAR KLEIN–GORDON EQUATION

with 0 , « ! 1 and k is a specified constant. When w 5 «u, the above equation reduces to

and

f5 5

O u 5 2 sin x sin t 4

utt 2 c2uxx 1 c2u 2 «2su3 5 0

i

i 50

S

2 sin x t 2

3

5

7

9

(3.9)

D

t t t t 1 2 1 . 3! 5! 7! 9!

u(x, 0) 5 cos ks, ut(x, 0) 5 0.

Equation (3.13) is an example of (1.1) with g being a cubic function and f (x, t) 5 0. Equation (2.14) implies that

This explains the phenomena that 2 sin x sin t is the self-canceling ‘‘noise’’ term. Further, canceling the noise term we obtain inductively the exact solution to (3.1) given by u(x, t) 5 sin x sin t.

(3.10)

Table I shows the errors obtained upon solving the Klein– Gordon equation (3.1) using only three terms of the decomposition method. In Table I and thereafter, error is understood to be 5 utrue value 2 approximate valueu. Hence, in Table I error 5 usin x sin t 2 f3u. It is noted that only three terms were needed to obtain an error of less than 1%. The overall errors can be made much smaller by adding new terms of the decomposition. Our next example (see [20, 24]) focuses on a nonlinear problem that cannot be solved explicitly. The standard method that is used to handle such an equation is the perturbation method which reduces the underlying equation to a class of linear partial differential equations. We will present the decomposition method as an alternative for solving the equation.

u(x, t) 5 cos ks 1 c2Lt21Lxu 2 c2Lt21u 1 «2sLt21(u3).

wtt 2 c wxx 1 c w 2 sw 5 0, 2

3

(3.14)

For convenience later we set g2 5 c2k2 1 c2, which is the dispersion relation for the linear Klein–Gordon equation. y The iterates u(x, t) 5 on50 un , upon using (2.15) and (2.17), are given by u0 5 cos kx

(3.15)

and u1 5 c2Lt21Lxu0 2 c2Lt21u0 1 «2sLt21(u30)

(3.16) 1 22 2 2 5 2 g t cos kx 1 « s t (Dk cos kx 1 Ak cos 3kx) 1 O(«3). 2!

Similarly, u2 5 c 2Lt21Lxu1 2 c 2Lt21u1 1 «2sLt21(3u20 u1) 5

EXAMPLE 2. Consider the hyperbolic equation 2

(3.13)

1 1 44 g t cos kx 2 «2st 4[3g2 cos kx 4! 4!

(3.17)

1 (3k2c2 1 c2) cos 3kx] 1 O(«3) (3.11) and

where c, c, s are appropriate physical constants, with the initial conditions

u3 5 c2Lt21Lxu2 2 c 2Lt21u2 1 «2sLt21(3u20u2 1 3u0u21) 52

w(x, 0) 5 « cos kx, wt(x, 0) 5 0, 2y , x , y (3.12)

«2s 6 1 66 g t cos kx 1 t [75g4 cos kx 6! 4.6!

(3.18)

1 (129c4k4 1 90c2c2k2 1 25c4) cos 3kx] 1 O(«3). TABLE I

Finally,

Error Obtained Using Decomposition Method with Three Iterations x

t 5 0.2

t 5 0.4

t 5 0.6

0.2 0.4 0.6

2.6 3 1024 2.10 3 1023 7.03 3 1023

5.2 3 1024 4.12 3 1023 1.377 3 1022

7.5 3 1024 5.98 3 1023 1.9964 3 1022

u4 5 c2Lt21Lxu3 2 c2Lt21u3 1 «2sLt21(3u20 u3 1 6u0 u1 u2 1 u31) 5

1 44 «2s 8 g t cos kx 2 t [39g6 cos kx 8! 2.7! 1 (13c6 1 93c2c4k4 1 54c4c2k2 1 84c6k6) cos 3kx] 1 O(«3).

(3.19)

446

DEEBA AND KHURI

Combining these iterates yields

S

12

with the following initial conditions:

D

1 22 1 44 1 66 1 88 g t 1 g t 2 g t 1 g t cos kx 2! 4! 6! 8

F F

1 «2s 1 «2s 1 52

G

3 2 3 24 75 4 6 39 6 8 t 2 gt 1 gt 2 g t cos kx 8 4! 4.6! 2.7! 1 2 1 t 2 (3k2c2 1 c 2)t 4 8 4!

1 (129c4k4 1 90c2k2c2 1 25c4)t6 4.6!

G

1 84c6k6)t 8 cos 3kx 1 O(«3). Inductively, the first term of the series is cos gt cos kx while the other terms, upon comparing them with the results of the perturbation method, yield

1

F

9s t sin gt 32g

G

3s (cos gt 2 cos 3gt) cos kx 128g2

1 «2

F

utt 2 uxx 2 sin u 5 0

G

where

l2 5 9c2k2 1 c2. Adomian’s method led to the same approximate solution obtained by the perturbation method. While in the perturbation method, one needs to solve a class of linear partial differential equations and match the boundaries to obtain the solution; in the decomposition method, once the problem is properly set, one only needs to differentiate and integrate to obtain Adomian’s polynomials An and consequently u.

(3.24)

which is a model for the nonlinear meson fields with periodic properties for the unified description of mesons and their particle sources (see [21]). Simple solutions of (3.24) may be found in which u is a function of j 5 (x 2 vt)/ Ï1 2 v2. In this case, the field equation (3.24) reduces to the pendulum-like equation d2u/dj 2 5 sin u. For a real physical velocity v , 1, there are implicit solutions given by sin As u 5 6sech( j 2 j0) which have particle number N 5 61 and total energy E 5 (2/f)(1/ Ï1 2 v2). These may be interpreted as the fields associated with a particle of mass 2/f centered at j 5 j0 and moving with velocity v [21]. We will show, using Adomian’s decomposition method, how to obtain solutions that coincide with the implicit solutions of Eqs. (3.22)–(3.23) given by

(3.21) 3s s (cos gt 2 cos lt) 1 128c2k2 128c2

(cos lt 2 cos 3gt) cos 3kx 1 O(«3),

(3.23)

Equation (3.22) arises from the general sine–Gordon equation

(3.20)

1 (13c 6 1 93c 2c4k4 1 54c 4c2k2 2.7!

u(x, t) 5 cos gt cos kx 1 «2

du (0) 5 22. dt

u(0) 5 f,

sin As u 5 sech t.

(3.25)

We will first describe the decomposition scheme as it applies to the general sine–Gordon equation (3.24). If we set Lt 5 ­2 /­t2, Lx 5 ­2 /­x2, N(u) 5 sin u, then Eq. (2.14) implies that the solution of (3.22) can be expressed in the operator form u(x, t) 2 u(x, 0) 2 tut(x, 0) 5 Lt21Lxu 1 Lt21N(u). (3.26) y

Thus writing u(x, t) 5 on50 un , and N(u) 5 sin u 5 y on50 An then, using (3.26) and (2.16), the various iterates are given by un11 5 Lt21Lxun 1 Lt21An , n $ 0

(3.27)

u0 5 u(x, 0) 1 tut(x, 0).

(3.28)

Our third example involves the sine–Gordon equation. EXAMPLE 3. Consider the pendulum-like equation d 2u 5 sin u dt 2

(3.22)

with

447

SOLVING THE NONLINEAR KLEIN–GORDON EQUATION

For N(u) 5 f (u) 5 sin u and upon using (2.9), we have

In a like manner, we find upon using the identity cos(f 2 2t) 5 2cos 2t,

A0 5 sin u0 u2 5 Lt21A1 5 Lt21(u1 cos u0)

A1 5 u1 cos u0 A2 5 u2 cos u0 2

5 2Lt21(u1 cos 2t)

u21 sin u0 2!

u31 A3 5 u3 cos u0 2 u2u1 sin u0 2 cos u0 3!

(3.34)

5 dsG t 2 Ak sin 2t 1 Ak t cos 2t 2 hfA sin 2t cos 2t

(3.29) and

. .

S

.

u3 5 Lt21A2 5 Lt21 u2 cos u0 2

Therefore, since u0 is known, Eqs. (3.27)–(3.28) provide y the series solution on50 un , where

D

u21 sin u0 2!

5 Lt21h(2cos 2t)u2 1 (2As sin 2t)u21j 5 hfG t 2 gjh FJ sin 2t 1 ask AD t cos 2t

(3.35)

2 hfA sin 2t cos 2t 1 dsA t sin 2t 2

u0 5 u(x, 0) 1 tut(x, 0)

1 hfA t cos2 2t 2 sd;f A cos2 2t sin 2t 1 aags A sin3 2t,

u1 5 Lt Lxu0 1 Lt A0 21

21

u2 5 Lt21Lxu1 1 Lt21A1 .

(3.30)

. un11 5 Lt21Lxun 1 Lt21An with the A9ns given in (3.29). We will adapt the above algorithm in (3.30) to the pendulum-like equation given in (3.22) with initial conditions (3.23), where u is only a function of t. For this special case, the polynomials An given in (3.29) are dependent on t only. The solution (3.27) thus reduces to un11 5 Lt21An , n $ 0,

(3.31)

where in (3.35) the identities cos(f 2 2t) 5 2cos 2t and sin(f 2 2t) 5 sin 2t are used. Higher iterates can be determined similarly. Upon combining the first six iterates and expanding in Taylor’s series around t 5 0, we obtain

u 5 f 2 2t 1

2 3 10 5 61 7 t 2 t 1 t 3! 5! 7!

2770 9 103058 11 2 t 1 t 1 ... 9! 11!

(3.36)

which coincide with Taylor expansion of (3.25). We observed that each time an iterate is added, the Taylor expansion coincides to the next higher term.

where Our final example deals with a nonhomogeneous Klein– Gordon equation having a quadratic nonlinear term. u0 5 u(0) 1 t

du (0) 5 f 2 2t. dt

(3.32)

EXAMPLE 4. Consider the Klein–Gordon equation of the form

Using (3.30) and (3.31)–(3.32) the various iterates are given as utt 2 uxx 1 u1 5 Lt21A0 5 Lt21(sin u0) 5 As t 2 Af sin 2t.

(3.33)

f f2 u 1 u2 5 x2 sin2 t 4 2

f u(x, 0) 5 0, ut(x, 0) 5 x. 2

(3.37)

448

DEEBA AND KHURI

Equation (2.14) implies that

TABLE III Errors Obtained Using Decomposition Method with Four Iterations

f2 u(x, t) 5 As fxt 1 Lt21Lxu 2 Lt21u 4

S

(3.38)

D

f 2 Lt21(u2) 1 Lt21 x2sin2 t . 2 y

Thus if u(x, t) 5 on50 un , then using (2.15)–(2.17) the first three terms are given by

S

f u0 5 As fxt 1 Lt21 x2 sin2 t 2 5 As fxt 2 As

S

2

D

(3.39)

D

x 2 1 Af t2 1 2 cos ft x2 f2 f

x 5 0.1

t 0.1 0.2 0.3 0.4 0.5

2.8 3.7 2.3 4.3 4.2

3 3 3 3 3

10214 10211 1029 1028 1027

x 5 0.2 7.0 2.5 1.8 3.6 3.6

3 3 3 3 3

10215 10211 1029 1028 1027

x 5 0.3 1.7 8.6 1.0 2.5 2.7

3 3 3 3 3

10214 10212 1029 1028 1027

x 5 0.4 4.6 1.2 4.2 8.1 1.2

3 3 3 3 3

10214 10211 10211 1029 1027

x 5 0.5 7.8 3.7 1.3 1.4 8.0

3 3 3 3 3

10214 10211 1029 1028 1028

decreases as the number of iterations increase. Convergence is rapid and the errors are extremely small. The presence of self-canceling noise terms such as 2As(x2 /2f 2) and Af t 2x2 in Eq. (3.39) affects the monotonicity of the errors (see [3]).

and ACKNOWLEDGMENTS

f2 u1 5 Lt Lxu0 2 Lt21u0 2 Lt21(u20). 4 21

(3.40)

REFERENCES

In a like manner, u2 5 Lt21Lxu1 2 u3 5 Lt21Lxu2 2

f 21 Lt u1 2 Lt21(2u0u1) 4 2

(3.41)

f 2 21 Lt u2 2 Lt21(u21 1 2u0u2). (3.42) 4

The exact solution of (3.37) is given by f u(x, t) 5 x sin t. 2

(3.43)

The exact solution was compared with the approximate solution using Adomian’s method. Tables II and III show the errors obtained by using the approximation f3 and f4 , as defined in (2.18), respectively. It is evident that the error

TABLE II Errors Obtained Using Decomposition Method with Three Iterations t 0.1 0.2 0.3 0.4 0.5

x 5 0.1 1.2 2.9 1.3 1.5 9.9

3 3 3 3 3

10211 1029 1027 1026 1026

x 5 0.2 6.3 4.2 3.3 4.9 4.7

The authors thank the referees for their valuable suggestions.

3 3 3 3 3

10211 1029 1029 1027 1026

x 5 0.3 1.2 1.2 1.6 7.7 1.8

3 3 3 3 3

10210 1028 1027 1027 1026

x 5 0.4 1.8 2.1 3.3 2.3 9.8

3 3 3 3 3

10210 1028 1027 1026 1026

x 5 0.5 2.4 3.1 5.3 4.0 1.9

3 3 3 3 3

10210 1028 1027 1026 1025

1. G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method (Kluwer, Dordrecht, 1994). 2. G. Adomian and R. Rach, J. Math. Anal. Appl. 174, 118 (1993). 3. G. Adomian and R. Rach, Math. Appl. 24(11), 61 (1992). 4. G. Adomian, Comput. Math. Appl. 21(5), 101 (1991). 5. G. Adomian, Comput. Math. Appl. 21(5), 133 (1991). 6. G. Adomian and R. Rach, Comput. Math. Appl. 10(12), (1990). 7. G. Adomian, Nonlinear Stochastic Systems Theory and Applications to Physics (Kluwer, Dordrecht, 1989). 8. G. Adomian, J. Math. Anal. Appl. 115, 233 (1986). 9. G. Adomian, J. Math. Anal. Appl. 113, 202 (1986). 10. G. Adomian, J. Math. Anal. Appl. 102, 420 (1984). 11. G. Adomian, Comput. Appl. Math. 11(2), 225 (1984). 12. R. E. Bellman and G. Adomian, Partial Differential Equations— Methods for their Treatment and Solutions (Reidel, Dordrecht, 1984). 13. Y. Cherruault, G. Saccomandi, and B. Some, Math. Comput. Modeling 16(2), 85 (1992). 14. Y. Cherruault, Kybernetes 18(2), 31 (1989). 15. B. K. Datta, Comput. Math. Appl. 20(1), 61 (1990). 16. B. K. Datta, J. Math. Anal. Appl. 142, 6 (1989). 17. P. Kuo and L. Vazquez, J. Appl. Sci. 1, 25 (1983). 18. J. L. Lions, Quelques Me´thods de Re´solution des Proble`mes aux Limites Non Line´aires (Dunod, Paris, 1969). 19. W. Cao and B. Guo, J. Comput. Phys. 108, 296 (1993). 20. J. K. Perring and T. H. R. Skyrme, Nucl. Phys. 31, 550 (1962). 21. L. I. Schiff, Phys. Rev. (2) 84(1), 1 (1951). 22. W. Strauss and L. Vazquez, J. Comput. Phys. 28, 271 (1978). 23. G. B. Whitham, Linear and Nonlinear Waves (Wiley, New York, 1974). 24. E. Zauderer, Partial Differential Equations of Applied Mathematics (Wiley, New York, 1983).