Parametric envelope solitons in coupled Korteweg-de Vries equations

Parametric envelope solitons in coupled Korteweg-de Vries equations

10 March 1997 PHYSICS LETTERS A Physics Letters A 227 ( 1997) 47-54 ELSEVIER Parametric envelope solitons in coupled Korteweg-de Vries equations...

741KB Sizes 1 Downloads 72 Views

10 March 1997

PHYSICS

LETTERS

A

Physics Letters A 227 ( 1997) 47-54

ELSEVIER

Parametric envelope solitons in coupled Korteweg-de Vries equations Georg Gottwald, Roger Grimshaw, Boris Malomed’ Department of Mathematics,

Monash University, Clayton, Vie. 3168, Australia

Received 19 August 1996; revised manuscript received 17 December 1996; accepted for publication 28 December 1996 Communicated by A.R. Bishop

Abstract

We demonstrate that a system of linearly coupled Korteweg-de Vries equations, which inter alia is a general model of resonantly coupled internal waves in a stratified fluid, can give rise to broad envelope solitons produced by a double phaseand group-velocity resonance between the fundamental and second harmonics for certain wavenumbers. We derive asymptotic equations for the amplitudes of the two harmonics, which are identical to the second-harmonic-generation equations in a diffractive medium, that have recently attracted a lot of attention in nonlinear optics and give rise to the so-called parametric solitons. To check if the predicted solitons are close to exact solutions of the coupled Korteweg-de Vries equations, we perform direct numerical simulations, with initial conditions suggested by the above-mentioned parametric-sol&on solution to the asymptotic equations. Since the latter is known only in a numerical form, we use for them a recently developed analytical variational approximation. As a result, we observe very long-lived steadily propagating wave packets generated by these initial conditions. Thus we find a physical system that may allow experimental observation of propagating parametric solitons, while in nonlinear optics they are observed only as spatial solitons.

1. Introduction Envelope solitons supported by parametric interactions [ 11, i.e., quadratic nonlinear interactions coupling three underlying harmonic waves, or only two waves in the case of second-harmonic generation (SHG) in diffractive media, have recently attracted a lot of attention in nonlinear optics [2-4], after being introduced long ago in the pioneering paper by Karamzin and Sukhorukov [S] . These solitons were observed experimentally in the form of stationary self-supported localized beams in a three- or two-dimensional bulk (in the former case, the beam I Permanent address: Department of Interdisciplinary Studies, Faculty of Engineering, Tel Aviv University,Tel Aviv 69978, Israel.

is cylindrical) [6]. However, they have never been observed in the temporal domain, i.e., as propagating solitons. A principal difficulty is that it is virtually impossible to fabricate a sufficiently long optical fibre with a quadratic nonlinearity, which would be a natural medium to support a temporal-domain soliton. In this work we aim to demonstrate that such parametric solitons in the time-domain can be easily launched in other physical systems, where they may be interesting objects in their own right. The model which will be considered in this work is a system of two coupled Korteweg-de Vries (KdV) equations, ut + uxxx - 6uux = ux,

0375-9601/97/$17.00 Copyright @ 1997 Published by Elsevier Science B.V. All rights reserved PII SO375-9601(97)00021-2

(1.1)

48

4 +

G. Gottwald et al/Physics

+ vu, ~~x,x

- 6,uvv, =

KU,.,

(1.2)

where v is the long wave phase speed mismatch between the two components, S and p are relative dispersion and nonlinear coefficients, respectively, and the coupling constant K is normalized to unity in the sequel. This system appears in a theory describing nearresonant interaction of internal waves in a stratified fluid [7], and also for the analogous near-resonant interaction of planetary waves [8]. The relative dispersion S plays a crucial role. In principle, it may be either positive or negative. When S < 0, there is a gap in the spectrum of the phase speeds of the linearized system, in which envelope gap-solitons, very different from the classical KdV solitary waves, may stably exist [ 91. It is also noteworthy that for the case of negative relative dispersion (6 < 0) but with opposite sign of the linear coupling (i.e. K = - 1) intersection of the dispersion curves of the corresponding uncoupled KdV subsystems may give rise, instead of the above-mentioned gap in the phase-speed spectrum, to a gap in the wavenumber spectrum, which actually implies a novel fundamental instability in the coupled-wave system [ lo]. If, on the other hand, the relative dispersion is positive (which is the more typical case for applications to internal waves), the gap is not opened, and strictly speaking solitons do not exist at all in this model, since, whatever the velocity of a wave packet, one can always find a linear wave with the same phase velocity, which can be expected to give rise to the emission of these linear waves by the nonlinear wave packet. Nevertheless, we will show below that this model can support, through the SHG process, very long-lived approximate envelope solitons of a new type. They can exist as stationary solutions within the framework of an asymptotic expression for the envelope amplitudes, very slowly decaying into linear waves. Before proceeding to analysis of the SHG resonance in ( 1.1) and ( 1.2), it is relevant to note that the same system can generate a sort of a parametric soliton in a simpler way, viz., through the interaction between a fundamental and zero (rather than second) harmonic. Indeed, let us take, for simplicity, a purely symmetric system with 6 = ,CL= 1, 77 = 0 and we recall that K = 1. In this case, the spectrum of the linearized system is very simple, c = fl - k2, where c is the phase

Letters A 227 (1997) 47-54

velocity, k is the wavenumber, and the alternate signs correspond to two different branches of the dispersion curve. The corresponding group velocity as a function of k is v = ~!~tl- 3k2.

(1.3)

If we now start the analysis with a harmonic with wavenumber k, its quadratic self-interaction generates the zero harmonic. For the latter one, the group velocity is Z!Z~(see ( 1.3) ) . Picking up the lower sign, one notices that it coincides with the group velocity of the fundamental harmonic belonging to the other branch at k2 = i. Thus, the fundamental harmonic at this wavenumber can be resonantly coupled to the zero harmonic (“mean field”). After a straightforward analysis, the system of asymptotic equations for the slowly varying complex envelope 111of the fundamental and the real-valued mean field VO take the form i(Vl),

+ ~(VI>..

(Uo)t + dI~I12L

+

UOVI

= 0,

=

0,

(1.4) (1.5)

where the variables x and t are scaled versions of the original variables in a frame moving with the common group velocity, and (+ is a constant. These equations are exactly integrable and support an obvious soliton solution similar to the soliton of the nonlinear Schrijdinger equation [ 111. In this work we will consider within the framework of the general system ( 1.l ) and ( 1.2)) another way to support a parametric soliton, which is based on the resonance between the fundamental and second harmonics (SH). In the lowest-order approximation, this resonance demands that the phase and group velocities of the fundamental and second harmonics coincide for certain wavenumbers, which proves to be possible in general, but incidentally is not possible in the above-mentioned symmetric system. The rest of the paper is organized as follows. In Section 2, we analyze the dispersion curves for the general system and find the wavenumbers where this double resonance can occur. In Section 3, which is the technical core of this work, we derive by means of a standard multiscale technique, an asymptotic system for the fundamental and secondary envelopes (similar to ( 1.4) and ( 1.5) for the resonance with the zero harmonic). In their eventual form, these equations coin-

49

C. Gottwald et al. /Physics Letters A 227 (I 997) 47-54

tide with “canonical” equations that give rise to parametric solitons in nonlinear optics [ 2,3]. Therefore, we can expect existence of the known one-parametric family of stationary soliton solutions [ 3,4 1. However, as we mentioned above, the coupled KdV equations with dispersions of the same sign (i.e. 6 > 0), which is the case under consideration, cannot strictly speaking have exact soliton solutions because there is no gap in the spectrum of the linear waves. Usually, the solitons obtained as exact solutions to the asymptotic amplitude equations exist for very long times as almost exact solitons to the underlying system (here the coupled KdV equations ( 1.1) and ( 1.2)), suffering from an extremely slow radiative decay (see, e.g., Ref. [ 91) . In order to understand the meaning of these solutions, in Section 4 we report results of direct numerical simulations of these coupled KdV equations, ( 1.1) and ( 1.2). As an initial state, we take wave forms suggested by the soliton solutions of the asymptotic equations. The result is that, over fairly long times, these initial wave forms undergo practically no evolution except for translation at a constant velocity (an arbitrary localized wave packet completely decays on the same time scale). An important practical problem in doing these simulations is representing the soliton solutions to the asymptotic amplitude equations, because they are known in an exact analytical form for only a single parameter value [ 51, while at all other parameter values they were found numerically. However, in Ref. [ 41 it was demonstrated that one can approximate the general soliton solutions to a very reasonable accuracy by analytical expressions based on a variational approximation. This is exactly the form in which we take solutions to the asymptotic equations in order to generate the initial states for our direct simulations. The final result is that we indeed obtain “almost genuine” envelope solitons governed by the linearly coupled KdV equations ( 1.1) and ( 1.2). Concluding remarks are collected in Section 5.

Wk,)

(2.1)

= 4X,),

where k, is the resonant wavenumber. Here w(k) is the dispersion relation, generally multi-valued with two or more branches, which may depend on several parameters. In order to guarantee a sufficiently strong interaction the group velocities must also coincide so that $(k,)

= $(2k,).

(2.2)

In general, the two resonance conditions impose a certain relation between the parameters. For the coupled Korteweg-de Vries equations ( 1.l > and ( 1.2) the linear dispersion relation is (w+k3)(w+cYk3

-vk)

-k’=O.

(2.31

The two solution branches are given by 2w=-(1+S)k3+77kfk&l-6)k2+77]*+4, (2.4) where we recall that here K = 1 in the system ( 1.1) and ( 1.2). We shall assume that 6 > 0, so that each KdV component in the coupled system ( 1.1) and ( 1.2) has the same-signed dispersion. The resonance conditions (2.1) and (2.2) can be met if the first and second harmonics belong to the lower and upper branches respectively. We find that the resonant wavenumber k, is given by

k= r

5 7?‘+4 1)’ J MS-

(2.5)

provided that the parameters

satisfy the condition

, 3 s+1 516=

d

l-

16 v* -25q2+4’

(2.6)

and we require that v( 6 - 1) > 0, which only for a nonzero value of v, and for typical example of the dispersion curves Fig. 1 at the same parameter values which in Section 4 for the numerical simulations.

is possible S + 1. A is shown in we will use

2. Dispersion curves and the double resonance 3. Derivation of the amplitude The condition for second-harmonic resonance (i.e. equality of the phase-velocities of the fundamental and second harmonics) requires that

equations

To investigate the dynamics of the interaction of a fundamental harmonic with the wavenumber k, with

G. Gottwald et al./Physics Letters A 227 (1997) 47-54

50

Thus, at the lowest order, 0( e2), (3.2) is trivially satisfied by the linear dispersion relation for each mode, yielding ‘... ,/

*-

\

,.:

a_,Q

&(k,)

j

‘:.

,,,’

‘:

= -+,

h(kd

= -~

‘..,

::

w2 + k; k:!



,,I ‘.,, -12

-

/

‘.,

/’ -1,

‘.,

:/

,:’ .18

‘$,.

..:

At the next order, 0( e3), linearity again allows us to separate the two modes, and we obtain for the first harmonic component

i

“8.2

.,5

45

.1

f

I

0.1

(5

2

Fig. 1. Dispersion curves, c as a function of k, where c = o/k, for the parameter values ,u = 7) = 1 and the corresponding resonant value of the relative dispersion S = 4.6. The resonant wavenumber is kr = 0.9.

its second harmonic with wavenumber 2k, in the system ( 1.1) and ( 1.2), we will perform a multiscale perturbation analysis for the resonant waves. Thus, introducing a small parameter E and setting X = EX, Tt = Et, T2 = e’t, we write u = e2A(X,T)eio1+E2B(X,T)ei82+e3u3+e4u4+c.c., u = e2,f1A(X,T)eieI

+

(3.1)

where tYl = k,x - ult,

E

=Jz/(u,lJ),

(3.2)

0 where the linear operator matrix is LU = a, +a,,,,

L,, =

-a,,

Lu = -a,, L,,= a,+ 7)a, + 6a,,, and the nonlinear N,(u,u)

=6uux,

dt

k, zAx

>>

e iBI

-AT, - $(kl)Ax)e”‘,

(3.3)

and

(3.4)

below. There are similar expressions for the second harmonic on the right-hand side. Since the eigenfunctions of the homogeneous adjoint operator are proportional to eisr and eiB2,we can set u3 = 03 = 0, to obtain

e2 = k2x - 02t,

with the carrier-wavenumbers kl , k2 and oi = w( ki), i = 1,2 and T = (Tl,Tz,...). Recall that we choose the lower (upper) branch, respectively, for i = 1 and i = 2. Here kl s k, + dk and k2 3 2(k, + EAk) = 2kI, where Ak allows a small deviation from the exact resonance-condition (2.2). The system ( 1.1) and ( 1.2) can be written as e

=

-AT~+ #Ax + 61Ax +

Note that here we are regarding h( ki) as an operator, with ki ----)ki-ieax, for reasons which we shall discuss

l2&B(X,T)ei”2

+ E3U3 + 64tJ4 + C.C.,

=

right-hand-side

aw1 AT, + ak

= -ts(k,)AkAx

+O(e2), (3.5)

&-, +

2

(2kr) Bx = -E-j-ga2w2 (2kr)2AkBx

+ O(e2>.

(3.6) The corrections on the right-hand side will be absorbed into higher-order terms of the expansion. We note that the ansatz ( 3.1) , with ZJi( ki ) interpreted as an operator as described above ensures vanishing of the right-hand sides of (3.3) and (3.4) at the leading order, and so ~3 = ~3 = 0. Here, of course, the form (3.1) implicitly includes an 0( e3) correction due to [i( ki). Next, the form of (3.5) and (3.6) implies that we make the transformation

is

N~(u,v) =~,uuv,~.

(k,)Ax

where

G. Gottwald

vG

5

%(k,)

et al. /Physics

= X do’ (2k,),

because at the leading order, A and B are functions of X’ and T2. Henceforth we omit the primes on X’. At the next order, O( l 4), separation of the two waves is no longer possible because of the nonlinear terms. Thus, we have to take into account the nonlinearity-induced phase-mismatch

51

Letters A 227 (1997) 47-54

The compatibility conditions for solution of the coupled system of Eqs. (3.8) and (3.9) yields a pair of evolution equations for the amplitudes A and B, iAr2 i- @Axx

+ 2iQlAkAx

i&, +@zBxx

+4i@zAkBx+

where the coefficients

+ NIA*BeiX

= 0,

N2A2eeiX = 0,

(3.10)

are defined by (3.1 I)

2 -(Ak)’

(3.12) (3.13)

Using the explicit form of 5 as described obtains for the resonant components,

=

above, one (3.14) The phase-mismatch of the nonlinear eliminated by the transformation

-A,+~-gy-i “W’ (k,)Axx

A = A’exp( -iAkX

SW, - z(k,)AkAx

eiBl .

+

-&

+

-; 3

-

a2w2 z(2kr)2AkBx

+

6ik

B = B’exp(-2iAkX

2 Wr)

- iaiT2), - iazT2).

Azeie2-ix

(3.15) (3.16)

Substituting this transformation into the system (3.10)) the phase-mismatch can be removed by choosing

Bxx

eiQ2 2c-q - (+2 = (4@2 - 2@,) (Ak)‘.

r

terms can be

+

(jik,A*Beiel+iX,

(3.17)

(3.8) Finally, we obtain, omitting

where the asterisk stands for the complex conjugate. Analogously, one obtains the second equation

the primes,

iAT2+@lAxx+S,A+NIA*B=0, i&z + QizBxx + S2B + N2A2 = 0,

(3.18)

where S1 = g1 + Qlt (Ak)‘, S, = ~2 + 4@2(Ak)’

+ 6pi k&f ( k,) A2eie2-ix + 6,uikJj

( k,)~2(2k,)A*BeiB1fiX.

(3.9)

= 2Sl.

(3.19)

We now see that we can set Si = SZ = 0 by choosing = -4@2( Ak)2, which satisfies (+l = -@l(Ak)2,(72 the constraint (3.17). It is readily seen that in terms of the new variables, the effective phases are 01 = i 02 = k, - w ( k,) t and are evaluated at precisely the resonant wavenumber, satisfying both resonance conditions (2.1) and (2.2). However, it is useful to retain St and S2 in (3.18)) which is equivalent to adding the

52

G. Gottwald et al/Physics

terms -SITZ and -S2T2 to the phases 81 and f32, respectively. Such terms in effect replace the resonant frequencies wi( kr) with wi( k,) + e2Si, i = 1,2. Thus, each Si can be regarded as an 0(e2) frequency detuning term, for which the parameters in the coupled KdV equations ( 1.1) and ( 1.1) fail to satisfy (2.6) by terms of 0( l 2) . Note that (3.19) implies that we can either satisfy the resonance condition (2.1) exactly, i.e. we can choose Ak = 0 while introducing a frequency detuning via the remaining free parameter ~1, or we can satisfy (2.2) exactly introducing the detuning via the wavenumber mismatch Ak. We will choose the frequency detuning and set Ak = 0. It is important to mention that stationary-wave solutions of (3.18) only exist in the detuned case St # 0. In the exactly resonant case we can obtain oscillatory solitary wave solutions, for which the phase of the oscillations introduces a necessary term balancing the dispersion in the asymptotic limit. To obtain stationary solitary wave solutions of (3.18)) we must assume that @t St < 0 and @2S2 < 0. Then, using the transformations

A=

d

w2

A’

,

N1

(3.20)

+ Bxx - .cJB+ iA2 = 0,

where

For the system ( 1.1) and ( 1.2) we have @i < 0 and @2 < 0 and so 5 > 0. In the stationary case a particular solitary wave solution of (3.20) for 5 = 1 was found in Ref. [ 51 to be

=f

where the parameters a2 = (p+

(3.21)

were found to be

1)(2p+y)(y+5)

b= (P+wmT,

y_

4P2 . 1-P

42 The parameter p is determined of the cubic equation

as a real positive root

(3.22)

2op3 + (4 - 35)p2 + 45p = 5.

In the parameter range of interest, there is only one real root, which is positive, i.e. p is uniquely determined.

3x4 2 cosh’( x/2)



B,(X)

=

3 2 cosh2 (x/2)

simulations

In this section we address the question of the validity of our theoretical results; to be specific we investigate whether the coupled system of equations (3.20) and their approximate solutions (3.2 1) is an appropriate representation for an SHG resonance in the full coupled KdV equations ( 1.1) and ( 1.1) . To show this, we integrate the coupled KdV equations numerically using the approximate solution (3.21) as an initial condition, suitably transformed back into the original variables of the KdV system. We use a pseudo-spectral code, where the linear terms are treated with a semiimplicit Crank-Nicholson scheme, and the nonlinear terms with an explicit leapfrog. There are only two free parameters, /J and r], since 6 is given through (2.6) as

T2 = +,T’,

iAT+Axx-A+A*B=O,

A,(X)

B = be-Yxz,

A = ae-Px2,

4. Numerical



where r = f 1 for Nt N&Q @2 > 0 or < 0 respectively, and omitting the primes, the system (3.18) is cast into the final form,

;iJBr

For 5 > 1 localized solutions of the system (3.20) can be obtained using an asymptotic expansion around 1/l (see, e.g., Ref. [ 21) . For arbitrary 4’ # 1, Steblina et al. [4] have found good agreement of the numerically obtained solution with a Gaussian ansatz, whose parameters were determined by means of a variational approximation. We will use their result for the ansatz

BJ?!!p

2N1 N2@1

x=

Letters A 227 (1997) 47-54



s = -3 + J(97?2 3 - J(9$

+ lOO)/($

+ 4)

+ lOO)/(?+

+ 4) ’

(4.1)

which implies 6 > 4. Note that ,u = 0 is also a possible case, i.e. one equation of the system ( 1.1) and ( 1.1) can be linear. It is pertinent to mention that when the parameters take 0( 1 )-values the extension of the

G. Gortwald et al./Physics

Letters A 227 (1997) 47-54

Fig. 2. Left panel: initial condition. right panel: evolved wavepacket at

t = 20. The upper pictures refer to

53

U, the lower to 6. The parameters

are chosen as p = 71 = 1, and the corresponding S = 4.6. The resonant wavenumber is k, = 0.9.

wavepacket is very large, whereas the amplitude is very small to satisfy EAk < k,. Fig. 2 shows the initial condition determined as described above, and its evolution. During the evolution, the wavepacket readjusts slightly to get closer to an exact solution. It has been checked that over this time an arbitrary initial wave packet not satisfying Eqs. (3.20) or the system ( 1.1) and ( 1.1) strongly disperses (see Fig. 3). Similar numerical experiments have been performed at other parameter values, all revealing that the soliton solution to (3.20) approximates very closely a solution of the original coupled KdV system.

5. Conclusion In this work we have analyzed a resonance between the fundamental and second harmonics in the linearly coupled KdV system ( 1.1) and ( 1.2), which is a general model of resonant interaction of internal waves in stratified fluids. We have shown that the resonance may indeed take place. By means of a mul-

tiscale expansion, we have derived asymptotic envelope amplitude equations which coincide with the standard second-harmonic generation equations in nonlinear optics. We have then performed direct simulations of the coupled KdV equation, with initial states corresponding to the parametric solitons of the asymptotic equations. The latter were taken in the form of a Gaussian ansatz provided by the variational approximation. The result is that these initial states give rise to almost exact envelope solitons of the coupled KdV equation. This provides the possibility to directly produce parametric solitons in the temporal domain, which is not possible in available nonlinear optical media. Finally, we note that in realising the physical applicability of the theoretical results obtained here, we should recall that the coupled KdV equations ( 1.1) and ( 1.2) have been derived under a long-wave hypothesis. Eqs. ( 1.l ) and ( 1.2) are in nondimensional form, but, for example their derivation in the internal wave context requires that kdh < 1 where kd is the dimensional wavenumber, and h is a length scale for the vertical stratification (e.g. the total fluid depth).

54

G. Gottwald et al. /Physics Letters A 227 (1997) 47-54

I

“5w

I -PO

.2a

-ram

0

rca

a0

3c.l

4m

I

“~‘L

-x.3

4m

.x0

‘:,

-,m

0

,a,

I 2.m

Fig. 3. Left panel: initial condition, right panel: evolved wavepacket at r = 20. The upper pictures refer to u, the lower to u. Same parameters as in Fig. 2, but with a disturbed p.

The condition for second-harmonic resonance (2.5) requires that the nondimensional wavenumber k, be of order unity, which implies that the dimensional resonant wavenumber krd should be finite, but k,d < 1. This condition can readily be achieved in oceanic, atmospheric or laboratory situations. Acknowledgement

One of the authors, G. Gottwald, is deeply indebted to the “Deutscher Akademischer Austauschdienst”, which sponsored this work through the HSP II/AUPE, and to the CRC-Southern Hemisphere Meteorology for their funding. References [ I] Y.R. Shen, The principles of nonlinear optics (Wiley, New York, 1984). [2] G.I. Stegeman, M. Sheik-Bahae, E. Van Stryland and G. Assanto, Opt. Lea. 18 (1993) 13; R. Schick, J. Opt. Sot. Am. B 10 (1993) 1848;

M.J. Werner and PD. Drummond, J. Opt. Sot. Am. B 10 (1993) 2390; L. Tamer, C.R. Menyuk and Cl.1 Stegeman, Opt. Lett. 19 (1994) 1617; C.R. Menyuk, R. Schick and L. Tamer, J. Opt. Sot. Am. B 11 ( 1994) 2434; C. Etrich, U. Peschel, E Lederer and B.A. Malomed, Phys. Rev. A 52 (1995) R3444. [3] A.V. Buryak and Y.S. Kivshar, Opt. Lett. 19 (1994) 1612: Phys. Len. A 197 (1995) 407. [4] V. Steblina, Y.S. Kivshar, M. Lisak and B.A. Malomed, Opt. Commun. 118 (1995) 345. [ 5 ] Y.N. Karamzin and Al? Sukhomkov, Pis’ma Zh. Eksp. Teor. Fiz. 20 ( 1974) 734 [JETP Lett. 20 ( 1974) 3391. [ 61 W.E. Tomrellas, Z. Wang, D.J. Hagan, E.W. Van Stryland, G.I. Stegeman, L. Tomer and CR. Menyuk, Phys. Rev. Lett. 74 (1995) 5036; R. Schick, Y. Baek and G.I. Stegeman, Phys. Rev. E 53 (1996) 1138. [7] J. Gear and R. Grimshaw, Stud. Appl. Math. 70 ( 1984) 235 [ 81 H. Mitsudera, J. Atmos. Sci. 51 3137; G. Gottwald and R. Grimshaw, in preparation. [9] R. Grimshaw and B.A. Malomed, Phys. Rev. Len. ( 1993); R. Grimshaw, B.A. Malomed and Xin Tian, Phys. Lett. A 201 (1995) 285. [lo] R. Grimshaw, J. He and B.A. Malomed, in preparation. [ 111 Y.C. Ma, Stud. Appl. Math 59 (1978) 201.