- Email: [email protected]
Computing for pattern forming systems in nonlinear optics
___ k!!B
-
&g
ELSEVIER
ATI-IEMATICS AND COMPUTERS N SIMULATION Mathematics and Computers in Simulation 37 (1994) 385-403
Computing for pattern forming systems in nonlinear optics R.A. Indik *, A.C. Newell Department of Mathematics, University ofArizona, Tucson, AZ 85721, USA
Abstract Methods used to compute spectral method is presented systems where a modulational
solutions of nonlinear PDEs arising in nonlinear optics are described. In addition a for integrating dispersion dominated PDEs that works particularly well in integrating ansatz is appropriate.
1. Introduction
The human brain perceives the world through patterns. Indeed, sometimes our minds are so intent on discovering patterns that we may perceive them even when they are absent! We wish to see structure and regularity in the world, as that allows us to predict and categorize. Thus it is not surprising that the study of physical patterns has become very popular. Nor is it surprising that this study has been very fruitful. After all, patterns are the language which allows us to move from the particular to the general. This study has revealed seemingly universal behaviors of patterns from widely disparate areas of investigation. What is a pattern? Intuitively, a pattern is some regularity, some repetitive structure. The mathematical notion that can make this more precise is that of symmetry. Just any symmetry is not enough to make a pattern. We insist that the symmetry should allow us to deduce the structure of the whole from that of some small part. In terms of symmetry, we insist that the symmetries of the pattern can generate the entire pattern from some finite region. The simple patterns that we will restrict attention to will have length preserving symmetries only. All of the above assume a regularity more perfect that we really wish to insist on in physical phenomena. We must slightly relax our insistence on symmetry, and replace it with a local notion that holds almost everywhere. We must still have symmetry, but it may no longer be perfect, there may be “defects” in our patterns where the underlying symmetry fails (even locally). We are not aware of a completely general rigorous definition of patterns that satisfies these requirements, but as with pornography, we know it when we see it. * Corresponding
author.
Elsevier Science B.V. SSDI 0378-4754(94)00026-G
386
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
The hallmark of pattern formation is spontaneous symmetry breaking. One starts with some system with a very high degree of symmetry, typically the system is invariant under all affine transformations. There is some solution of the system which shares that full symmetry, but that symmetric solution is not stable for all parameter values. At some critical parameter value, the symmetric solution becomes unstable. Stable nonsymmetric solutions arise. Any time we have a solution of a symmetric system, we know that images of that solution under the action of the symmetries are also solutions. If the solution is not fixed by those symmetries there are many solutions to choose among. This “choice” that the systems must make has a great deal to do with the defects that appear in the patterns. Many times a system will “choose” differently in different spatial regions, and the interfaces between these choices must lead to defects in the patterns. In many cases the systems that are studied are only approximately or locally symmetric, and some higher order terms, or boundary conditions break the symmetry. In these cases one must understand how the patterns can fit boundary conditions, and where the higher order terms become important. This paper concerns itself primarily with the numerical experience we have had in studying one and two dimensional patterns in nonlinear optics. Nonlinear optics provides a particularly promising area for the study of patterns because of the wide range of time (and space) scales that naturally arise. Maxwell’s equations for optical fields suggest time scales of lo-i5 seconds and spatial scales of 10e7 meters, and even the envelope equations most frequently employed lead to time scales of order lO_” seconds. This should be compared to fluid convection, or chemical reaction diffusion problems, where the underlying time scale is such that often the experiments that yield patterns must be run for many hours. With nonlinear optical problems, we may expect that similar results might be available in microseconds. Indeed, often the slow scale in our pattern evolution is quite fast in practical terms, challenging the abilities of detectors. The study of nonlinear optical patterns is still at quite an early stage, and at this point the numerical results out-number the experimental results. Nonetheless, experiments demonstrating a variety of patterns have been performed [l-lo] and this is a very active area of investigation. The dominant role of diffraction in nonlinear optics require us to resolve different numerical issues than are faced in most other areas of pattern simulation. Most of those areas are well modeled by equations for which diffusion is quite central, dominating the spatial coupling which define the patterns. For many optical problems, diffusion may be entirely absent. This makes the similarities in the patterns observed even more striking, while making the numerical issues distinct. Some of the earliest work in the study of nonlinear optical pattern formation is in the context of passive media, and an excellent survey of work through 1989 in this field can be found in [ll]. In much of this work the pattern arises when transverse freedom is added to a bistable system [12], where the bistability is the result of feedback. Either a broad or (theoretically) infinite beam provides the light source, and the interaction is either unconfined, or confined in only one dimension. In this scenario, diffraction couples the behavior of the light field at nearby points. If the material in which the interaction takes place supports diffusion, there is diffusion as well. Stationary patterns including rolls and hexagons have been seen [14,13]. In addition to pattern forming systems based on bistability, there is also a physically simpler, though analytically more complex system which can form both hexagonal and square patterns [15,16]. In this system, two beams of light counterpropagate through a Kerr medium.
R.A. Indik, A.C. Newell/Mathematics
and Computers in Simulation 37 (1994) 385-403
387
Active media produce patterns as well, and the study of pattern formation in lasers has been lively. We [17-211 have made extensive theoretical and numerical studies of the spatio-temporal dynamics of the Maxwell-Bloch system describing a two level laser with one and two transverse dimensions, as well as of a model of a Raman laser, and we have seen a rich variety of pattern formation mechanisms reminiscent of the behavior of the Rayleigh-Bernard convection in fluids. Work in these areas relies heavily on numerical experiments. Ideally one might prefer that theorists perform analysis, and make predictions which can be motivated by, and checked against real experiments. What is the role of numerical experiment? We search for universality, and thus seek models and equations that are somewhat independent of the details of the physics and still produce the phenomena of interest. We must first predict the behavior of those equations, then we would like to compare that behavior to the behavior of physical systems. Physical experiments may be too real to start with. Particularly in studying symmetry effects, it is very hard to maintain the required high degree of symmetry in an experiment, and yet one needs to know what the behavior of the symmetric system is before one can understand the more realistic, but less symmetric experimental system. Thus one numerically studies artificial, idealized numerical systems, incorporating part of the physics of the problem, to compare the results of the predictions of theory to the numerical results incorporating the included portion of the physics. In addition, the theoretical predictions are often in the form of supposedly simpler equations that should describe the behavior of the system. These equations can rarely be understood without recourse to numerical integration. Thus we find that we must use numerical integration to both perform experiments and produce predictions from the theory. In what follows we shall discuss the methods that we chose to use, and the issues and experiences that have led to those choices for two particular pattern forming systems. The first is the Maxwell-Bloch [17,20,21] equations with transverse diffraction. These equations describe the evolution of a light field in an active medium (two level and Raman). They are the starting point for modeling lasers, and have proved quite rich, revealing structure similar to various convection problems in fluid mechanics. Figs. 1 and 2 show typical results of computations for these system. The second set of equations whose integration we shall discuss describe the behavior of counterpropagating light fields in a Kerr material. Physically, this is an extremely simple system whose key elements, counterpropagation and nonlinearity, are ubiquitous in optical devices. Because of the extra dimension introduced by the counterpropagation it is analytically and numerically more complex than the bistable systems mentioned earlier. In fact, the study of some of those physically more complex systems was motivated by a desire to understand the behavior of this system. It is a system that reveals a surprisingly rich pattern forming capability, showing hexagonal or square patterns, as well as developing secondary instabilities [15,16]. Fig. 3 shows some typical patterns that were computed.
2. Computational
techniques
In this section we will describe various numerical schemes that we have used in studying the Maxwell-Bloch system, and the counterpropagating beam problem in Kerr media. In the case of
388
R.A. Indik. A.C. Newell /Mathematics
and Computers
in Simulation
37 (1994) 385-403
t= 215
Fig. 1. Frames from an animation showing the evolution of the intensity field of the Raman laser equations and periodic boundary conditions. Solid and dashed lines are equiphases of the field.
for 0 > 0
the Maxwell-Bloch system we will discuss not only our experience integrating the full system, but also our experience in integrating an amplitude equation that arises from that system. The integration of the partial differential equations that arise in studying patterns in nonlinear optics present a significantly different numerical challenge than that typically rlcountered in problems arising from fluid dynamics, chemistry or condensed matter physics. in nearly all of those systems, the spatial coupling is dominated by diffusion. In optical systems in the paraxial approximation, the spatial coupling is due to diffraction. Where diffusion is dominant, a numerical scheme must deal with the presence of very strong damping of high wavenumbers, that is, it must be able to handle a stiff system. Diffraction, on the other hand, produces fast oscillation for those high wavenumbers. In the diffusive case, it is appropriate to choose schemes which may not accurately track the fastest dynamics, but are correct in the sense that they damp high wavenumbers. Thus unconditionally stable schemes of type L, are often preferred. In the diffractive case, we feel it is inappropriate to use a scheme that replaces oscillation with damping. Instead, it is reasonable to allow the scheme to produce fast oscillation (of the wrong speed perhaps) as its limiting behavior for such high wavenumbers.
R.A. Indik, A.C. Newell /Mathematics
389
and Computers in Simulation 37 (1994) 385-403
Fig. 2. Frames from an animation showing the evolution of the intensity field of the Raman and periodic boundary conditions. Solid and dashed lines are equiphases of the field.
laser equations
for 0 < 0
This choice can be justified by analogy with the choice made for L, stable schemes, where very strong damping is approximated by (perhaps not as strong) damping. Fig. 4 illustrates the difference between the two cases, by showing the distribution of eigenvalues for the behavior of finite difference discretizations of the Maxwell-Bloch system ae -iaV2e = -ue + op, at
$+(l+iO)p=(r-n)e, t
(1)
+ bn = Re(p*e)
near a stable traveling wave solution as described FitzHugh-Nagumo type as described in [22] A=+-+3-h+L h, = (4 - 2h)/lO
in [19] and a reaction
diffusion
system
of (2) (3)
390
R.A. Indik, A. C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
20
-20
-40
-1
I -0.5
I 0
I 0.5
Stress Parameter D Fig. 3. Diagram showing computed patterns as a function counterpropagating light fields in a Kerr Medium.
of stress parameter
and grating
parameter
for periodic
R.A. lndik, A.C. Newell /Mathematics
391
and Computers in Simulation 37 (1994) 385-403
: .
I
-101
-20 L--L-u -2.0 -1.5
L
-1 .o
-0.5
Damping
0.0
0.5
-15
-10
-5 Damping
0
Fig. 4. Real (growthrate) and imaginary (frequency) parts of the eigenvalues of linearized discretized Maxwell-Bloch equations near a traveling wave solution and (b) FitzHugh-Nagumo types equation front solution.
5
operator for (a) near a traveling
near a stable traveling wave solution. Note that though high wavenumber modes are damped in the Maxwell-Bloch case, their damping tends towards a constant value (a), while the imaginary parts of the eigenvalues tend towards +a~. In the case of the reaction diffusion system, the eigenvalues associated with high wavenumbers tend to --. It is worth pointing out that the order parameter equations that arise from optical problems are no different from those in other areas. Thus there are no new numerical issues in the treatment of these approximations. The spatial derivative terms in the amplitude equations come from the complex dispersion relation of the underlying problem. Thus the diffusion coefficient comes from the curvature of the growth rate curve near threshold, and the “diffraction” coefficient comes from the curvature of the imaginary part of the eigenvalue. Which of these coefficients is dominant does not directly depend on the form of the spatial derivatives in the original model equations. 2.1. Maxwell-Bloch We have invested considerable effort in studying the behavior of Maxwell-Bloch and related systems, and have used a variety of methods to integrate these systems. Our numerical investigations have been in both one and two dimensions using periodic, reflecting and absorbing boundary conditions. We have modeled both finite beams and (theoretical) infinite systems. Any numerical scheme must deal with the presence of high frequencies, as mentioned above. In addition, the multiple scales that were taken advantage of to introduce amplitude equations must be resolved by our numerics. Typically there is some basic underlying spatial scale (i.e. k, = m). In addition there may be a second, much larger spatial scale on which we expect the pattern to be modulated. In order to say that we are seeing a pattern in the solution, we must see repetition. In order to study how that pattern fails to be perfect, (defects and
392
R.A. Indik, A.C. Newell /Mathematics
and Computers
in Simulation
37 (1994) 385-403
modulation), that repetitiveness must be varied. Thus we must resolve a large number of periods of the pattern. Quite a lot of memory and computation can be required. Multiple time scales are also present. Fortunately, in the case of the Maxwell-Bloch equations, the frequencies of the traveling wave solutions are known, so that one may make a change of variables to remove that frequency from the computation. One therefore need not be concerned with resolving that fast frequency and so can concentrate on tracking the modulation dynamics occurring on a much slower time scale. The new equations are ae --
iRe - iaV2e = -ue
+
up,
at
ap
n)e,
z
+p
= (Y -
p
+ bn = Re(
(4)
p*e)
if R > 0 or de
-
aP :Iat
- -&f2
+
l+i---
- iaV2e = -ue n
l+u
p = (r 1
-
+
n)e,
up,
(5)
G+bn = Re( p*e) for 0 < 0. One issue that is troubling when working with systems where there are multiple time scales is that if one uses a method that takes steps that are large compared to the finest scales, those scales cannot be resolved accurately. We can justify using such a method if the details of the behavior at the fastest scales do not significantly influence the slow scale behavior. Certainly, a prerequisite for this is that the total intensity in the fast scales should remain small. If, in addition, the scales are well separated, then one may argue that the faster dynamics has only weak influence, since the slower modes will be driven well off resonance. If, however, the scales are distributed continuously, as is the case in diffractive or diffusive problems, there is a real possibility of significant coupling between all scales. It is fortunate that in the case of the Maxwell-Bloch system it is straightforward to show [17] that the higher modes are all essentially decoupled from the slower dynamics of the problem. We have found that the problem is,insensitive to the inaccuracies for high wavenumbers, if that decoupling is captured by the method. An implicit scheme of Crank-Nicolson type works c.f. [23] quite well for integrating the one dimensional MB system. Start by discretizing the partial differential equations (4) using second order central differences, and integrate the resulting ordinary differential equations using a
R.A. Indik, A.C. Newell /Mathematics
Pad& (1, 1) approximation we get
and Computers in Simulation 37 (1994) 385-403
for the propagation
e;_1-2e)+ej+,
(trapezoid
+ il2el + a( pj - e;)
A2 -2ej+ej+,
A2
+ iRe, + a( pi -
rule). Specifically
393
in the case 0 > 0
(6) ej)
(8)
Here the primes denote the values at t + h, while the unprimed quantities are the values at time t. This is a nonlinear equation, which can be solved using Newton-Raphson. The Newton steps require that one solve linear systems with block tridiagonal structure with blocksize five. Such a system can be solved in time proportional to the blocksize and the number of points. Typically only two or three iterates are required, so the method is not too costly, and remains stable providing the Newton iteration converges. In addition, the limiting behavior for large imaginary eigenvalues is oscillatory. As discussed above, we feel this is desirable for this diffraction dominated problem. We choose time steps reflecting the modulation dynamics since we first transform the equations to remove the frequency of the expected dominant traveling wave. Without that transformation, our time step is limited by that frequency. If we could have, we would have applied this method in two spatial dimensions as well, and indeed, the method makes sense in that setting. However, in two dimensions we must solve a sparse linear system, whose band structure is much more complex. A straight-forward approach leads to a cost and a storage requirement for solving the linear system which is proportional to Nm where N is the total number of points in the discretization. Another method that we have used quite extensively in both one and two spatial dimensions is a split-operator scheme, known as the beam propagation method [24,25]. We use this approach when we have either periodic or reflecting boundaries. We split the operator into two Maxwell-Bloch. The algorithm then involves the alternapieces, diffraction, and “plane-wave” tion of steps that propagate the electric field “in vacua” (diffraction) with steps that ignore transverse spatial coupling and calculate the evolution of many “plane-waves” simultaneously. If one likesone can consider this algorithm physically, as replacing the two level medium with a large number of very thin slabs of two level media separated by linear material, such that the average values of gain etc. are the same for the original two level medium, and the stack of thin plates. The great advantage of using this sort of splitting is that the diffractive propagation can be performed exactly in spectral coordinates (take the spatial Fourier transform). The planewave propagation can be performed numerically with either an explicit or an implicit scheme, depending on the stiffness of that problem which is controlled by the values of the parameters b and (T. Another very important advantage of this approach is the simplicity of implementation. The propagation under diffraction involves calculating an FFT, multiplying by a precalcu-
394
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
lated multiplier array and then calculating the inverse FFT. The plane-wave propagation requires only the solution of identical uncoupled blocks of five ODES. The disadvantage of this approach has to do with the step size. The limit on step size in this method comes from the splitting error. The splitting error comes from the failure of the two halves of the operator to commute. Thus, to control the size of the error, the change in the calculated quantities in each half of the splitting must be kept small. This is unfortunate. We are studying the system in quasi-equilibrium states. The solutions tend to be modulations of the underlying traveling wave solutions, and that modulation is taking place on a scales in space and time that are long and slow. Thus most often the changes in those two halves will nearly cancel each other. The diffraction step of the propagation will lead to phase shifts proportional to l/Ax2 in the highest Fourier modes. If we insist that those phase changes remain small, we must use a step size as small as those required by traditional explicit finite difference methods. For the counterpropagating beam problem discussed below, this limitation is required for stability in the defocusing case [26]. Fortunately, this is not required for the Maxwell-Bloch equations, providing there is not significant energy in the high Fourier modes. A significant improvement in efficiency was realized when we split the operator so that the detuning term 0 in equation (4) was solved together with the diffraction in the spectral half of the propagation, rather than in the plane wave half. This meant that the expected pure traveling wave solution was a fixed point of both halves of the operator. Thus in simulations where wavenumbers stayed close to that critical wavenumber, large steps could be used. These were precisely the cases where a modulational ansatz was appropriate. That is, one could consider the solution as being constructed from the traveling wave solutions, with parameters such as wavenumber and amplitude varying slowly in space and time. It is always desirable to split the operator so as to make each half as small as possible on expected solutions. In this case the analytic theory is sufficiently precise to provide the tools needed to achieve a reasonable balance. In two dimensional cases where the boundary conditions that we wished to model were incompatible with spectral methods, we used finite differences to discretize the PDE’s into ODES, and then used the SDRIVE package [27] to integrate using an Adams-Bashforth with adaptive step size and order. We took advantage of PDELIB, a package developed by Hyman et al. [28] to provide flexibility in specifying boundary conditions, and the order of the differencing. 2.2. Application to Newell- Whitehead-Segel The amplitude equations that describe the behavior of the Maxwell-Bloch system near threshold is a Newell-Whitehead-Segel equation aB (u+l)%+2a =a@--l)B
l-f
[
--iaV2B+
1
(IB12)
ua2 2[2ig-&+$]2B (a+l)
(9)
R.A. Indik, A.C. Newell/Mathematics
and Computers in Simulation 37 (1994) 385-403
for (fi > 0) and a Complex Ginzburg-Landau _
2iaQzA
=
(TXOA
395
equation
-
(10) r0
with r() = 1 + .n2/(1 + u)’
and
x0 = 1 - in/(1
+ a)
for 0 < 0. When integrating these amplitude equations, we no longer need to worry about multiple scales-these equations only include the coarse scale. So the computational challenge is not as great. However, at least in the case R > 0, we have fourth derivatives in y. The derivation of the amplitude equations is based on the assumption that the first derivative in x and the second derivative in y are coming in at the same scale. If we take this to heart in our numerical integration, and scale the domains accordingly, we will not see inordinately severe restrictions on step size from the higher derivatives. However, it is frequently of interest to explore the behavior of the amplitude equations beyond the regions where the original scaling assumptions are valid. In these cases, one can end up requiring finer steps to integrate the amplitude equation than were required for the original Maxwell-Bloch equations. Our method of choice for integrating these two dimensional equations has been the beam propagation method, splitting the operator into a part including all of the spatial derivatives, and a part containing the rest of the terms in the equation. Using periodic boundary conditions, the derivative part of the operator can be solved for exactly, while the remaining portion is solved numerically. We have also frequently made use of PDELIB in integrating these equations, using the adaptive order and timestep Adams-Bashforth in the SDRIVE package. In general, as would be hoped, the amplitude equations present less of a numerical challenge. 2.3. Countelpropagation in Kerr media Counterpropagating beams in a Kerr material can spontaneously form two dimensional patterns [15,16]. This is the case both for focusing and defocusing Kerr media. We have been studying this system for some time, and have performed considerable numerical modeling of those equations. The equations we study are g
+ g
=ipQ:F+iq[
(F12+GIB12]F,
(11)
-E+z=i/3Q:B
(12)
with boundary conditions F(x, y,z=O,
t)=F,
and
B(x,y,z=l,
t)=B,.
Analytically and numerically, this system differs from the propagation problems such as the Maxwell-Bloch system in that the plane wave version of the system is a set of hyperbolic PDEs, where for the Maxwell-Bloch system we have ODES. We are dealing with an initial boundary value problem rather than an initial value problem. There are two distinct characteristic
396
R.A. Indik, A.C. Newell/Mathematics
Fig. 5. Characteristics
and Computers in Simulation 37 (1994) 385-403
of the plane wave counterpropagating
beam system.
directions and the boundary values of the two fields are specified at different points. We chose to integrate along the characteristics of the system as shown in Fig. 5. The operator was split into diffraction and plane wave propagation. The diffraction propagation was performed spectrally using the fast Fourier transform to convert from spatial to spectral coordinates. The remaining operator was integrated along the characteristics of the system. Fig. 5 shows the characteristics of the system. In the plane wave system, only the phase of a field changes along its characteristic. Integrating for that phase change yields
g
=B,
-B
eiD~GIFIZ+If12
1
eiD1B,12heiDG(IF~,~2+4/F,~2+~F,~2)h/4
(14)
where the integrals are along the characteristics for F and B respectively. The phase is accurate to third order, and the amplitude is exact. The overall accuracy of the scheme is second order due to the splitting error. Before we arrived at this choice of integration scheme, we used several other schemes. Originally we did not take advantage of the invariance of the intensity along the characteristics. Instead we used a trapezoid rule. We found that when we started taking advantage of intensity invariance, we could use many fewer points. In hindsight, this is not at all surprising. Enforcing the preservation of the intensity along the characteristics, leads to a scheme that numerically conserves the transverse integrals of the intensity. This is a conserved quantity in the full system, and there has been considerable experience to support the expectation of improved faithfulness when a numerical scheme can be chosen which preserves the conservation laws of a system [29]. When we first chose to integrate the phase change along the characteristics, we chose a trapezoidal approximation of the integral since we felt there was no need to go beyond a second order method. What we discovered, however, was that the output of our simulations developed an oscillation in z at the grid scale, and that this oscillation seemed to be present no matter how fine a time step or spacing we chose. We felt that if it were a numerical instability, it would disappear with a sufficiently fine step. Finally we realized that the structure of our algorithm decoupled the system into two independent systems. We were initializing our system
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
397
close to an unstable fixed point, and studying the long term behavior of the system. The instability that we were studying was degenerate. The system had to break symmetry, choosing a direction in which to move away from the unstable fixed point. Since the numerical system was actually two decoupled systems, and since we were perturbing with a small amount of white noise, the two systems made different choices, and it looked to us as if we had a numerical instability. This seems to us a natural difficulty, and we believe others will find our experience instructive. As mentioned earlier, the disadvantage of beam propagation is that the individual halves of the propagation must stay small. This is not a limitation in the early transient evolution of the counterpropagating beam system, but once the system begins to settle down into a fixed pattern, the two halves of our operator nearly cancel. It would be better if we could use a method that allowed us to take steps on the slower scale on which the pattern is evolving. We have experimented with splitting the operator differently, but the other splittings seem to fair more poorly. We suspect this is related to the conservation law no longer being preserved. Also, the method of propagation along characteristics may not be a good choice in this system. As implemented, we must choose our time steps and longitudinal space steps to be equal. Even a stationary plane wave solution requires that the time step be small if the longitudinal structure of the solution is rapidly varying.
3. A promising technique In this section we will present a technique that promises calculations considerably. It is applicable in cases where spectral should be particularly effective for dispersion dominated systems. it to pattern forming systems, but our expectation is that it improvements. In the following we will present the application model equation with diffraction.
to improve efficiency in our methods can be applied, and As of yet we have not applied will lead to very significant of this method for a simple
3.1. A model problem Consider a partial differential u, =
equation of the form:
iu,, + g(u)
(15)
with periodic boundary conditions, and specified initial conditions U(X, 0) = U,(X). We can take advantage of the dispersion law for the linear diffraction operator to remove the fast scales from the numerics. Consider the Fourier transform of (15) tii, f -iw(k)u,
+ @(U))k
where w(k) = k’. At this point we make the key assumption that the term (g(u)jk varying. Given that assumption, it is reasonable to transform to new variables +=yk
e- io(k)t
(16) is slowly
(17)
398
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
to get jk = eiw@)‘(f(y)),
(18)
where f(y) = g(u). In practice the evaluation of the right hand side of (18) involves recovering uk from y,, taking the inverse Fourier transform (using the FFT), calculating g and taking its Fourier transform. Solving this system amounts to integrating a large system of coupled ODES most of which have highly oscillatory right hand sides. We consider that sort of problem in the next section. 3.2. Oscillatory ODES Consider the ODE j =
fP’f( t)
(19)
where we assume that f is slowly varying, but w may take any real value (not necessarily large). The function f(t) on the right-hand side may depend on t through a dependence on y and should be thought of as arising from (18). We wish to integrate this equation in such a way that we have second order accuracy for all values of w. There are many possibilities, but for specificity, we shall choose one based on the trapezoid rule. We approximate f(t) on the interval [t,, tl] by f(t) = (4 - s)f(t,) with s = (t - t,)/(tl y(tJ
-y&J
+ (+ + s)f(t,)
- to) -
= [l’f(t)
E
(tl -
(20)
l/2 for s E [ - l/2,
l/2]. Then we have
eiwt dt
t,)/l’*
((3 - S)f(to)
+ (i + S)f(tl))
eiW((f~+‘~)/2+s(t~-f~)) ds
-l/2
=: ctl
_
toj
f(tl)
+f(to)
e2iBs
eio((t,+f,)/2)
ds
2 +
ctl
_
to)
(tl)-f(to) 2
eio((to+tI)/2)
“*
s e2iOs
ds
/ -l/2
(21) where t9= o(t, - to)/2 and ,(,)=(~+;(cosB-~)) The two curves in Fig. 6 show the behavior of the real and imaginary parts of j3.
(22)
R.A. Indik, A.C. Newell/Mathematics
-0.5
and Computers in Simulation 37 (1994) 385-403
399
‘.:
-1.0.
0
IO 20 30 40 50 e
Fig. 6. Real (solid) and imaginary (dotted) parts of fi as a function of 13= w(tl - t,)/2.
Note that p is bounded for all 8 and becomes small for large o. This accurately captures the sort of decoupling that occurs when there are large frequency mismatches, without requiring those frequencies to be resolved. It is possible to prove the following error estimate for this scheme Theorem 1. Giuen that j = e’“‘f(t> with f(t) E C*([t,, tll> and f”
t1 - to = ei4(to+ll)/*) --+(toM*(@
y(tl) _y(to)
+f(G(V)
+ c(t, - to)3F min
where 8 = w(t, - to)/2 and p(e)
(y
=
+
8- y
J-(COS
))
(24)
By varying the sort of approximation that is used for the slowly varying function f, one produces a variety of schemes. Almost any of the traditional schemes for integrating ODES may be mixed with this method to produce a scheme for integrating oscillatory ODES with known frequency. Thus one can also produce explicit schemes, or schemes of higher order. 3.3. Application to PDE’s In practice, one wishes to use this method without transforming to the derived variables y,. If we rewrite equation (21) applied to the system of equation (18) in terms of uk = yke-io(k)f we get 4,)
-
=
U/&o)
e
-i4h
e
+t,bQ)
-i4kKtI-to)
t1 -
to
-+u(WW
+
(25)
e-i4’o+f1)/*)_
t1 2
to du(toM*(e)
(26)
This implicitly defines uk(t,) as the solution of a nonlinear equation. We do not wish to apply
400
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
SPATIAL-TEMPORAL
SPECTRUM - FREE RUNNING
+540
ETUNING (GW
-540
+ 4.4
-4.4 WAVE NUMBER (km-‘) BROAD AREA LASER (Index Guided)
Fig. 7. Intensity
as a function
of wavenumber
and frequency
for a free running
broad area semi-conductor
laser.
R.A. Indik, A.C. Newell /Mathematics and Computers in Simulation 37 (1994) 385-403
401
Newton-Raphson, as the equations are fully coupled, so the Jacobian matrix is dense. However, when we calculate the Jacobian (Jac) of the left-hand side of this equation we get 1
_
e-i4(t0+fl)/2)____
t1 2
to
P(e) Jack>.
(27)
Providing the step size times the Jacobian of g(u) is small, we can use the identity matrix I as a good approximation to the Jacobian of the nonlinear equation we must solve implicitly. Thus our iteration becomes:
where the term rhs, refers the the right hand side of (26). The cost of each iteration is dominated by the cost of evaluating g(u), which requires a pair of fast Fourier transforms. Since the discretization of time introduces an error of order (tl - t,)3, the iteration need not be performed beyond the point where that accuracy is obtained. As a starting point for the iteration we choose Uk(t,)
= Uk(to)
e-io(k)(fl-fd+
to
t, e-iw((ro+r~)/2)_g(u(to))(p(e)
+p*(O)).
This corresponds to using the oscillatory version of forward Euler to provide our first guess. We have tested this method for a few simple cases so far, and have been quite impressed with its effectiveness. In the case of the Nonlinear Schrodinger equation, the following table exemplifies our preliminary experience. The first column shows the step size, the second column the error in the L2 norm of the answer, and the last column, the number of evaluations of the function g, that were required. Step size
L2 error
evaluations
l/Q
0.0018 0.00066 0.00019 0.000048 0.000012 0.0000030 0.00000075 0.00000019 0.000000047
56 49 65 101 164 320 543 1036 2048
l/4 l/8 l/16 l/32 l/64 l/128 l/256 l/512
The problem was integrated using 1024 points in a domain of size 8~. One might think that in a case where one knows the nonlinear dispersion relation, even better results could be obtained. For nonlinear Schriidinger, this was not the case. This can be understood since the difference between the nonlinear, and linear dispersion relationships is bounded (in terms of the amplitude for NLS), and in our simulation was order 1. It is our hope that this method will improve our efficiency considerably. Fig. 7 illustrates how the linear dispersion relationship can dominate the behavior of a nonlinear system well beyond threshold. The figure is a plot showing intensity as a function of wavenumber and frequency in
402
R.A. Indik, A.C. Newell/Mathematics
and Computers in Simulation 37 (1994) 385-403
a simulation of a broad-area semiconductor laser. The equations that describe its behavior are not dissimilar to the Maxwell-Bloch equations. You can clearly see that the intensity is clustered around a parabolic dispersion. In order resolve that picture, we had to resolve the full range of frequencies. We expect that using this method, we will be able to resolve frequencies corresponding to the thickness of the intensity spread about that dispersion curve.
4. Conclusion In the preceding we have attempted to make clear some of the issues involved in numerical exploration of pattern forming systems in Nonlinear Optics. In order to do so we have presented our view of what a pattern is, and illustrated that point of view with a couple of specific examples we have worked with coming from Nonlinear Optics. The challenges of multiple time and space scales inherent in pattern forming systems, and the particular challenges of dispersion dominated equation in that regime determine which methods are most useful. In addition we have presented a very promising technique for improving efficiency in integration of these systems, by taking advantage of the dispersion, that we believe will be most useful.
Acknowledgements
The authors would like to thank the Arizona Center for Mathematical Sciences (ACMS) for support. ACMS is sponsored by AFOSR contract FQ8671-9000589 with the University Research Initiative Program at the University of Arizona. In addition, Robert Indik would like to acknowledge many useful discussions with J.B. Geddes, E. Meron, J.V. Moloney, W.J. Firth and N.B. Abraham.
References [l] G. Grynberg, M. Pinard and P. Verkerk, Spontaneous symmetry breaking in a ring four-wave mixing oscillator, Europhys. Lett. 9 (1989) 139. [2] G. Grynberg and L.A. Lugiato, Quantum properties of hexagonal patterns, Optics Comm. 101 (1993) 69. [3] M. Vallet, P. Verkerk and G. Grynberg, Instabilities in a ring cavity filled with a Kerr medium and submitted to a bimode field. [4] T. Honda, Hexagonal pattern formation due to counterpropagation in KNbO,, Optics Lett. 18 (1993) 598. [5] J. Pender and H. Hesselink, IEEE J. Quant. Electron 25 (1989) 395; Optics Comm. 75 (1990) 47. [6] J. Pender and H. Hesselink, Degenerate conical emissions in atomic Sodium vapor, Opt. Sot. Am. B 7 (1990) 1361. [7] N. Tan-No, H. Hoshimiya and H. Inaba, IEEE J. Quant. Electron. 16 (1980) 147. [8] F.T. Arecchi, S. Boccaletti, P.L. Ramazza and S. Residori, Transition from boundaryto bulk-controlled regimes in optical pattern formation, Phys. Rev. Lett. 70 (1993) 2277-2280. [9] R. Macdonald and H.J. Eichler, Spontaneous optical pattern formation in a nematic liquid crystal with feedback mirror, Optics Comm. 89 (1992) 289-295.
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
403
[IO] S.R. Liu and G. Indebetouw, Periodic and chaotic spatiotemporal states in a phase-conjugate resonator using a photorefractive BaTiO, phase conjugate mirror, J. Opt. Sot. Am. B 9 (1992) 1509-1520. [ll] N.B. Abraham and W.J. Firth, Overview of transverse effects in nonlinear optical systems J. Opt. Sot. Am. B 7 (1990) 9.51-962. [12] H.M. Gibbs, Optical Bistability-Controlling Light with Light (Academic Press, Orlando, FL, 1985). [13] W.J. Firth, private communication. [14] W.J. Firth, A. Scroggie, G.S. MacDonald and L.A. Lugiato, Hexagonal patterns in optical bistability, Phys. Rev. A. 46 (1992) R3609-12. [15] R. Chang, W.J. Firth, R.A. Indik, J.V. Moloney and E.M. Wright, Three-dimensional simulations of degenerate counterpropagating beam instabilities in a nonlinear medium, Optics Comm. 88 (1992) 167. [16] J.B. Geddes, R.A. Indik, W.J. Firth and J.V. Moloney, Self focusing media like hexagons, self defocusing media prefer squares, preprint. [17] P.K. Jakobsen, J.V. Moloney, A.C. Newell and R. Indik, Space-time dynamics of wide gain section lasers Phys. Rev. A. 45 (1992) 8129-8137. [18] A.C. Newell and J.V. Moloney, Nonlinear Optics (Addison-Wesley, Redwood City, CA, 1992). [ 191 Q. Feng, J.V. Moloney and A.C. Newell, Amplitude Instabilities of Transverse Traveling Waves in Lasers, Phys. Rev. Lett. 71 (1993) 1705. [20] P.K. Jacobsen, J. Lega, Q. Feng, M. Staley, J.V. Moloney and A.C. Newell, Nonlinear transverse modes of large aspect ratio homogeneously broadened lasers: I. Analysis and numerics, Phys. Rev. A., submitted. [21] P.K. Jacobsen, J. Lega, Q. Feng, M. Staley, J.V. Moloney and A.C. Newell, Nonlinear transverse modes of large aspect ratio homogeneously broadened lasers: I. Pattern analysis near and beyond threshold, Phys. Rev. A., submitted. [22] A. Hagberg and E. Meron, Domain walls in nonequilibrium systems and the emergence of persistent patterns, Phys. Rev. E. 48 (1993) 705. [23] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Wiley, New York, 1967, 2nd Ed.). [24] F.D. Tappert, The parabolic approximation method, in: J.B. Keller and J.S. Papadakis, eds., Wave Propagation and Underwater Acoustics, Lecture Notes in Physics 70 (Springer, Berlin). [25] J.A. Fleck, J.R. Morris and M.D. Feit, Appl. Phys. 10 (1976) 129. [26] W.J. Firth, C. Penman and C. Pare, Transverse counterpropagating instabilities in Kerr slices, Optics Comm. 75 (1990) 136-140. [27] D. Kahaner, C. Moler and S. Nash, Numerical Methods and Software (Prentice Hall, Englewood Cliffs, NJ, 1989). [28] J.M. Hyman and B. Larrouturou, The numerical differentiation of discrete functions using polynomial interpolation methods, in: J.F. Thompson, ed., Numerical Grid Generation of Partial Differential Equations (Elsevier, New York, 1982) 487-506. [29] C.W. Gear, Invariants and numerical methods for ODES, Physica D 60 (1992) 303-310.
-
&g
ELSEVIER
ATI-IEMATICS AND COMPUTERS N SIMULATION Mathematics and Computers in Simulation 37 (1994) 385-403
Computing for pattern forming systems in nonlinear optics R.A. Indik *, A.C. Newell Department of Mathematics, University ofArizona, Tucson, AZ 85721, USA
Abstract Methods used to compute spectral method is presented systems where a modulational
solutions of nonlinear PDEs arising in nonlinear optics are described. In addition a for integrating dispersion dominated PDEs that works particularly well in integrating ansatz is appropriate.
1. Introduction
The human brain perceives the world through patterns. Indeed, sometimes our minds are so intent on discovering patterns that we may perceive them even when they are absent! We wish to see structure and regularity in the world, as that allows us to predict and categorize. Thus it is not surprising that the study of physical patterns has become very popular. Nor is it surprising that this study has been very fruitful. After all, patterns are the language which allows us to move from the particular to the general. This study has revealed seemingly universal behaviors of patterns from widely disparate areas of investigation. What is a pattern? Intuitively, a pattern is some regularity, some repetitive structure. The mathematical notion that can make this more precise is that of symmetry. Just any symmetry is not enough to make a pattern. We insist that the symmetry should allow us to deduce the structure of the whole from that of some small part. In terms of symmetry, we insist that the symmetries of the pattern can generate the entire pattern from some finite region. The simple patterns that we will restrict attention to will have length preserving symmetries only. All of the above assume a regularity more perfect that we really wish to insist on in physical phenomena. We must slightly relax our insistence on symmetry, and replace it with a local notion that holds almost everywhere. We must still have symmetry, but it may no longer be perfect, there may be “defects” in our patterns where the underlying symmetry fails (even locally). We are not aware of a completely general rigorous definition of patterns that satisfies these requirements, but as with pornography, we know it when we see it. * Corresponding
author.
Elsevier Science B.V. SSDI 0378-4754(94)00026-G
386
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
The hallmark of pattern formation is spontaneous symmetry breaking. One starts with some system with a very high degree of symmetry, typically the system is invariant under all affine transformations. There is some solution of the system which shares that full symmetry, but that symmetric solution is not stable for all parameter values. At some critical parameter value, the symmetric solution becomes unstable. Stable nonsymmetric solutions arise. Any time we have a solution of a symmetric system, we know that images of that solution under the action of the symmetries are also solutions. If the solution is not fixed by those symmetries there are many solutions to choose among. This “choice” that the systems must make has a great deal to do with the defects that appear in the patterns. Many times a system will “choose” differently in different spatial regions, and the interfaces between these choices must lead to defects in the patterns. In many cases the systems that are studied are only approximately or locally symmetric, and some higher order terms, or boundary conditions break the symmetry. In these cases one must understand how the patterns can fit boundary conditions, and where the higher order terms become important. This paper concerns itself primarily with the numerical experience we have had in studying one and two dimensional patterns in nonlinear optics. Nonlinear optics provides a particularly promising area for the study of patterns because of the wide range of time (and space) scales that naturally arise. Maxwell’s equations for optical fields suggest time scales of lo-i5 seconds and spatial scales of 10e7 meters, and even the envelope equations most frequently employed lead to time scales of order lO_” seconds. This should be compared to fluid convection, or chemical reaction diffusion problems, where the underlying time scale is such that often the experiments that yield patterns must be run for many hours. With nonlinear optical problems, we may expect that similar results might be available in microseconds. Indeed, often the slow scale in our pattern evolution is quite fast in practical terms, challenging the abilities of detectors. The study of nonlinear optical patterns is still at quite an early stage, and at this point the numerical results out-number the experimental results. Nonetheless, experiments demonstrating a variety of patterns have been performed [l-lo] and this is a very active area of investigation. The dominant role of diffraction in nonlinear optics require us to resolve different numerical issues than are faced in most other areas of pattern simulation. Most of those areas are well modeled by equations for which diffusion is quite central, dominating the spatial coupling which define the patterns. For many optical problems, diffusion may be entirely absent. This makes the similarities in the patterns observed even more striking, while making the numerical issues distinct. Some of the earliest work in the study of nonlinear optical pattern formation is in the context of passive media, and an excellent survey of work through 1989 in this field can be found in [ll]. In much of this work the pattern arises when transverse freedom is added to a bistable system [12], where the bistability is the result of feedback. Either a broad or (theoretically) infinite beam provides the light source, and the interaction is either unconfined, or confined in only one dimension. In this scenario, diffraction couples the behavior of the light field at nearby points. If the material in which the interaction takes place supports diffusion, there is diffusion as well. Stationary patterns including rolls and hexagons have been seen [14,13]. In addition to pattern forming systems based on bistability, there is also a physically simpler, though analytically more complex system which can form both hexagonal and square patterns [15,16]. In this system, two beams of light counterpropagate through a Kerr medium.
R.A. Indik, A.C. Newell/Mathematics
and Computers in Simulation 37 (1994) 385-403
387
Active media produce patterns as well, and the study of pattern formation in lasers has been lively. We [17-211 have made extensive theoretical and numerical studies of the spatio-temporal dynamics of the Maxwell-Bloch system describing a two level laser with one and two transverse dimensions, as well as of a model of a Raman laser, and we have seen a rich variety of pattern formation mechanisms reminiscent of the behavior of the Rayleigh-Bernard convection in fluids. Work in these areas relies heavily on numerical experiments. Ideally one might prefer that theorists perform analysis, and make predictions which can be motivated by, and checked against real experiments. What is the role of numerical experiment? We search for universality, and thus seek models and equations that are somewhat independent of the details of the physics and still produce the phenomena of interest. We must first predict the behavior of those equations, then we would like to compare that behavior to the behavior of physical systems. Physical experiments may be too real to start with. Particularly in studying symmetry effects, it is very hard to maintain the required high degree of symmetry in an experiment, and yet one needs to know what the behavior of the symmetric system is before one can understand the more realistic, but less symmetric experimental system. Thus one numerically studies artificial, idealized numerical systems, incorporating part of the physics of the problem, to compare the results of the predictions of theory to the numerical results incorporating the included portion of the physics. In addition, the theoretical predictions are often in the form of supposedly simpler equations that should describe the behavior of the system. These equations can rarely be understood without recourse to numerical integration. Thus we find that we must use numerical integration to both perform experiments and produce predictions from the theory. In what follows we shall discuss the methods that we chose to use, and the issues and experiences that have led to those choices for two particular pattern forming systems. The first is the Maxwell-Bloch [17,20,21] equations with transverse diffraction. These equations describe the evolution of a light field in an active medium (two level and Raman). They are the starting point for modeling lasers, and have proved quite rich, revealing structure similar to various convection problems in fluid mechanics. Figs. 1 and 2 show typical results of computations for these system. The second set of equations whose integration we shall discuss describe the behavior of counterpropagating light fields in a Kerr material. Physically, this is an extremely simple system whose key elements, counterpropagation and nonlinearity, are ubiquitous in optical devices. Because of the extra dimension introduced by the counterpropagation it is analytically and numerically more complex than the bistable systems mentioned earlier. In fact, the study of some of those physically more complex systems was motivated by a desire to understand the behavior of this system. It is a system that reveals a surprisingly rich pattern forming capability, showing hexagonal or square patterns, as well as developing secondary instabilities [15,16]. Fig. 3 shows some typical patterns that were computed.
2. Computational
techniques
In this section we will describe various numerical schemes that we have used in studying the Maxwell-Bloch system, and the counterpropagating beam problem in Kerr media. In the case of
388
R.A. Indik. A.C. Newell /Mathematics
and Computers
in Simulation
37 (1994) 385-403
t= 215
Fig. 1. Frames from an animation showing the evolution of the intensity field of the Raman laser equations and periodic boundary conditions. Solid and dashed lines are equiphases of the field.
for 0 > 0
the Maxwell-Bloch system we will discuss not only our experience integrating the full system, but also our experience in integrating an amplitude equation that arises from that system. The integration of the partial differential equations that arise in studying patterns in nonlinear optics present a significantly different numerical challenge than that typically rlcountered in problems arising from fluid dynamics, chemistry or condensed matter physics. in nearly all of those systems, the spatial coupling is dominated by diffusion. In optical systems in the paraxial approximation, the spatial coupling is due to diffraction. Where diffusion is dominant, a numerical scheme must deal with the presence of very strong damping of high wavenumbers, that is, it must be able to handle a stiff system. Diffraction, on the other hand, produces fast oscillation for those high wavenumbers. In the diffusive case, it is appropriate to choose schemes which may not accurately track the fastest dynamics, but are correct in the sense that they damp high wavenumbers. Thus unconditionally stable schemes of type L, are often preferred. In the diffractive case, we feel it is inappropriate to use a scheme that replaces oscillation with damping. Instead, it is reasonable to allow the scheme to produce fast oscillation (of the wrong speed perhaps) as its limiting behavior for such high wavenumbers.
R.A. Indik, A.C. Newell /Mathematics
389
and Computers in Simulation 37 (1994) 385-403
Fig. 2. Frames from an animation showing the evolution of the intensity field of the Raman and periodic boundary conditions. Solid and dashed lines are equiphases of the field.
laser equations
for 0 < 0
This choice can be justified by analogy with the choice made for L, stable schemes, where very strong damping is approximated by (perhaps not as strong) damping. Fig. 4 illustrates the difference between the two cases, by showing the distribution of eigenvalues for the behavior of finite difference discretizations of the Maxwell-Bloch system ae -iaV2e = -ue + op, at
$+(l+iO)p=(r-n)e, t
(1)
+ bn = Re(p*e)
near a stable traveling wave solution as described FitzHugh-Nagumo type as described in [22] A=+-+3-h+L h, = (4 - 2h)/lO
in [19] and a reaction
diffusion
system
of (2) (3)
390
R.A. Indik, A. C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
20
-20
-40
-1
I -0.5
I 0
I 0.5
Stress Parameter D Fig. 3. Diagram showing computed patterns as a function counterpropagating light fields in a Kerr Medium.
of stress parameter
and grating
parameter
for periodic
R.A. lndik, A.C. Newell /Mathematics
391
and Computers in Simulation 37 (1994) 385-403
: .
I
-101
-20 L--L-u -2.0 -1.5
L
-1 .o
-0.5
Damping
0.0
0.5
-15
-10
-5 Damping
0
Fig. 4. Real (growthrate) and imaginary (frequency) parts of the eigenvalues of linearized discretized Maxwell-Bloch equations near a traveling wave solution and (b) FitzHugh-Nagumo types equation front solution.
5
operator for (a) near a traveling
near a stable traveling wave solution. Note that though high wavenumber modes are damped in the Maxwell-Bloch case, their damping tends towards a constant value (a), while the imaginary parts of the eigenvalues tend towards +a~. In the case of the reaction diffusion system, the eigenvalues associated with high wavenumbers tend to --. It is worth pointing out that the order parameter equations that arise from optical problems are no different from those in other areas. Thus there are no new numerical issues in the treatment of these approximations. The spatial derivative terms in the amplitude equations come from the complex dispersion relation of the underlying problem. Thus the diffusion coefficient comes from the curvature of the growth rate curve near threshold, and the “diffraction” coefficient comes from the curvature of the imaginary part of the eigenvalue. Which of these coefficients is dominant does not directly depend on the form of the spatial derivatives in the original model equations. 2.1. Maxwell-Bloch We have invested considerable effort in studying the behavior of Maxwell-Bloch and related systems, and have used a variety of methods to integrate these systems. Our numerical investigations have been in both one and two dimensions using periodic, reflecting and absorbing boundary conditions. We have modeled both finite beams and (theoretical) infinite systems. Any numerical scheme must deal with the presence of high frequencies, as mentioned above. In addition, the multiple scales that were taken advantage of to introduce amplitude equations must be resolved by our numerics. Typically there is some basic underlying spatial scale (i.e. k, = m). In addition there may be a second, much larger spatial scale on which we expect the pattern to be modulated. In order to say that we are seeing a pattern in the solution, we must see repetition. In order to study how that pattern fails to be perfect, (defects and
392
R.A. Indik, A.C. Newell /Mathematics
and Computers
in Simulation
37 (1994) 385-403
modulation), that repetitiveness must be varied. Thus we must resolve a large number of periods of the pattern. Quite a lot of memory and computation can be required. Multiple time scales are also present. Fortunately, in the case of the Maxwell-Bloch equations, the frequencies of the traveling wave solutions are known, so that one may make a change of variables to remove that frequency from the computation. One therefore need not be concerned with resolving that fast frequency and so can concentrate on tracking the modulation dynamics occurring on a much slower time scale. The new equations are ae --
iRe - iaV2e = -ue
+
up,
at
ap
n)e,
z
+p
= (Y -
p
+ bn = Re(
(4)
p*e)
if R > 0 or de
-
aP :Iat
- -&f2
+
l+i---
- iaV2e = -ue n
l+u
p = (r 1
-
+
n)e,
up,
(5)
G+bn = Re( p*e) for 0 < 0. One issue that is troubling when working with systems where there are multiple time scales is that if one uses a method that takes steps that are large compared to the finest scales, those scales cannot be resolved accurately. We can justify using such a method if the details of the behavior at the fastest scales do not significantly influence the slow scale behavior. Certainly, a prerequisite for this is that the total intensity in the fast scales should remain small. If, in addition, the scales are well separated, then one may argue that the faster dynamics has only weak influence, since the slower modes will be driven well off resonance. If, however, the scales are distributed continuously, as is the case in diffractive or diffusive problems, there is a real possibility of significant coupling between all scales. It is fortunate that in the case of the Maxwell-Bloch system it is straightforward to show [17] that the higher modes are all essentially decoupled from the slower dynamics of the problem. We have found that the problem is,insensitive to the inaccuracies for high wavenumbers, if that decoupling is captured by the method. An implicit scheme of Crank-Nicolson type works c.f. [23] quite well for integrating the one dimensional MB system. Start by discretizing the partial differential equations (4) using second order central differences, and integrate the resulting ordinary differential equations using a
R.A. Indik, A.C. Newell /Mathematics
Pad& (1, 1) approximation we get
and Computers in Simulation 37 (1994) 385-403
for the propagation
e;_1-2e)+ej+,
(trapezoid
+ il2el + a( pj - e;)
A2 -2ej+ej+,
A2
+ iRe, + a( pi -
rule). Specifically
393
in the case 0 > 0
(6) ej)
(8)
Here the primes denote the values at t + h, while the unprimed quantities are the values at time t. This is a nonlinear equation, which can be solved using Newton-Raphson. The Newton steps require that one solve linear systems with block tridiagonal structure with blocksize five. Such a system can be solved in time proportional to the blocksize and the number of points. Typically only two or three iterates are required, so the method is not too costly, and remains stable providing the Newton iteration converges. In addition, the limiting behavior for large imaginary eigenvalues is oscillatory. As discussed above, we feel this is desirable for this diffraction dominated problem. We choose time steps reflecting the modulation dynamics since we first transform the equations to remove the frequency of the expected dominant traveling wave. Without that transformation, our time step is limited by that frequency. If we could have, we would have applied this method in two spatial dimensions as well, and indeed, the method makes sense in that setting. However, in two dimensions we must solve a sparse linear system, whose band structure is much more complex. A straight-forward approach leads to a cost and a storage requirement for solving the linear system which is proportional to Nm where N is the total number of points in the discretization. Another method that we have used quite extensively in both one and two spatial dimensions is a split-operator scheme, known as the beam propagation method [24,25]. We use this approach when we have either periodic or reflecting boundaries. We split the operator into two Maxwell-Bloch. The algorithm then involves the alternapieces, diffraction, and “plane-wave” tion of steps that propagate the electric field “in vacua” (diffraction) with steps that ignore transverse spatial coupling and calculate the evolution of many “plane-waves” simultaneously. If one likesone can consider this algorithm physically, as replacing the two level medium with a large number of very thin slabs of two level media separated by linear material, such that the average values of gain etc. are the same for the original two level medium, and the stack of thin plates. The great advantage of using this sort of splitting is that the diffractive propagation can be performed exactly in spectral coordinates (take the spatial Fourier transform). The planewave propagation can be performed numerically with either an explicit or an implicit scheme, depending on the stiffness of that problem which is controlled by the values of the parameters b and (T. Another very important advantage of this approach is the simplicity of implementation. The propagation under diffraction involves calculating an FFT, multiplying by a precalcu-
394
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
lated multiplier array and then calculating the inverse FFT. The plane-wave propagation requires only the solution of identical uncoupled blocks of five ODES. The disadvantage of this approach has to do with the step size. The limit on step size in this method comes from the splitting error. The splitting error comes from the failure of the two halves of the operator to commute. Thus, to control the size of the error, the change in the calculated quantities in each half of the splitting must be kept small. This is unfortunate. We are studying the system in quasi-equilibrium states. The solutions tend to be modulations of the underlying traveling wave solutions, and that modulation is taking place on a scales in space and time that are long and slow. Thus most often the changes in those two halves will nearly cancel each other. The diffraction step of the propagation will lead to phase shifts proportional to l/Ax2 in the highest Fourier modes. If we insist that those phase changes remain small, we must use a step size as small as those required by traditional explicit finite difference methods. For the counterpropagating beam problem discussed below, this limitation is required for stability in the defocusing case [26]. Fortunately, this is not required for the Maxwell-Bloch equations, providing there is not significant energy in the high Fourier modes. A significant improvement in efficiency was realized when we split the operator so that the detuning term 0 in equation (4) was solved together with the diffraction in the spectral half of the propagation, rather than in the plane wave half. This meant that the expected pure traveling wave solution was a fixed point of both halves of the operator. Thus in simulations where wavenumbers stayed close to that critical wavenumber, large steps could be used. These were precisely the cases where a modulational ansatz was appropriate. That is, one could consider the solution as being constructed from the traveling wave solutions, with parameters such as wavenumber and amplitude varying slowly in space and time. It is always desirable to split the operator so as to make each half as small as possible on expected solutions. In this case the analytic theory is sufficiently precise to provide the tools needed to achieve a reasonable balance. In two dimensional cases where the boundary conditions that we wished to model were incompatible with spectral methods, we used finite differences to discretize the PDE’s into ODES, and then used the SDRIVE package [27] to integrate using an Adams-Bashforth with adaptive step size and order. We took advantage of PDELIB, a package developed by Hyman et al. [28] to provide flexibility in specifying boundary conditions, and the order of the differencing. 2.2. Application to Newell- Whitehead-Segel The amplitude equations that describe the behavior of the Maxwell-Bloch system near threshold is a Newell-Whitehead-Segel equation aB (u+l)%+2a =a@--l)B
l-f
[
--iaV2B+
1
(IB12)
ua2 2[2ig-&+$]2B (a+l)
(9)
R.A. Indik, A.C. Newell/Mathematics
and Computers in Simulation 37 (1994) 385-403
for (fi > 0) and a Complex Ginzburg-Landau _
2iaQzA
=
(TXOA
395
equation
-
(10) r0
with r() = 1 + .n2/(1 + u)’
and
x0 = 1 - in/(1
+ a)
for 0 < 0. When integrating these amplitude equations, we no longer need to worry about multiple scales-these equations only include the coarse scale. So the computational challenge is not as great. However, at least in the case R > 0, we have fourth derivatives in y. The derivation of the amplitude equations is based on the assumption that the first derivative in x and the second derivative in y are coming in at the same scale. If we take this to heart in our numerical integration, and scale the domains accordingly, we will not see inordinately severe restrictions on step size from the higher derivatives. However, it is frequently of interest to explore the behavior of the amplitude equations beyond the regions where the original scaling assumptions are valid. In these cases, one can end up requiring finer steps to integrate the amplitude equation than were required for the original Maxwell-Bloch equations. Our method of choice for integrating these two dimensional equations has been the beam propagation method, splitting the operator into a part including all of the spatial derivatives, and a part containing the rest of the terms in the equation. Using periodic boundary conditions, the derivative part of the operator can be solved for exactly, while the remaining portion is solved numerically. We have also frequently made use of PDELIB in integrating these equations, using the adaptive order and timestep Adams-Bashforth in the SDRIVE package. In general, as would be hoped, the amplitude equations present less of a numerical challenge. 2.3. Countelpropagation in Kerr media Counterpropagating beams in a Kerr material can spontaneously form two dimensional patterns [15,16]. This is the case both for focusing and defocusing Kerr media. We have been studying this system for some time, and have performed considerable numerical modeling of those equations. The equations we study are g
+ g
=ipQ:F+iq[
(F12+GIB12]F,
(11)
-E+z=i/3Q:B
(12)
with boundary conditions F(x, y,z=O,
t)=F,
and
B(x,y,z=l,
t)=B,.
Analytically and numerically, this system differs from the propagation problems such as the Maxwell-Bloch system in that the plane wave version of the system is a set of hyperbolic PDEs, where for the Maxwell-Bloch system we have ODES. We are dealing with an initial boundary value problem rather than an initial value problem. There are two distinct characteristic
396
R.A. Indik, A.C. Newell/Mathematics
Fig. 5. Characteristics
and Computers in Simulation 37 (1994) 385-403
of the plane wave counterpropagating
beam system.
directions and the boundary values of the two fields are specified at different points. We chose to integrate along the characteristics of the system as shown in Fig. 5. The operator was split into diffraction and plane wave propagation. The diffraction propagation was performed spectrally using the fast Fourier transform to convert from spatial to spectral coordinates. The remaining operator was integrated along the characteristics of the system. Fig. 5 shows the characteristics of the system. In the plane wave system, only the phase of a field changes along its characteristic. Integrating for that phase change yields
g
=B,
-B
eiD~GIFIZ+If12
1
eiD1B,12heiDG(IF~,~2+4/F,~2+~F,~2)h/4
(14)
where the integrals are along the characteristics for F and B respectively. The phase is accurate to third order, and the amplitude is exact. The overall accuracy of the scheme is second order due to the splitting error. Before we arrived at this choice of integration scheme, we used several other schemes. Originally we did not take advantage of the invariance of the intensity along the characteristics. Instead we used a trapezoid rule. We found that when we started taking advantage of intensity invariance, we could use many fewer points. In hindsight, this is not at all surprising. Enforcing the preservation of the intensity along the characteristics, leads to a scheme that numerically conserves the transverse integrals of the intensity. This is a conserved quantity in the full system, and there has been considerable experience to support the expectation of improved faithfulness when a numerical scheme can be chosen which preserves the conservation laws of a system [29]. When we first chose to integrate the phase change along the characteristics, we chose a trapezoidal approximation of the integral since we felt there was no need to go beyond a second order method. What we discovered, however, was that the output of our simulations developed an oscillation in z at the grid scale, and that this oscillation seemed to be present no matter how fine a time step or spacing we chose. We felt that if it were a numerical instability, it would disappear with a sufficiently fine step. Finally we realized that the structure of our algorithm decoupled the system into two independent systems. We were initializing our system
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
397
close to an unstable fixed point, and studying the long term behavior of the system. The instability that we were studying was degenerate. The system had to break symmetry, choosing a direction in which to move away from the unstable fixed point. Since the numerical system was actually two decoupled systems, and since we were perturbing with a small amount of white noise, the two systems made different choices, and it looked to us as if we had a numerical instability. This seems to us a natural difficulty, and we believe others will find our experience instructive. As mentioned earlier, the disadvantage of beam propagation is that the individual halves of the propagation must stay small. This is not a limitation in the early transient evolution of the counterpropagating beam system, but once the system begins to settle down into a fixed pattern, the two halves of our operator nearly cancel. It would be better if we could use a method that allowed us to take steps on the slower scale on which the pattern is evolving. We have experimented with splitting the operator differently, but the other splittings seem to fair more poorly. We suspect this is related to the conservation law no longer being preserved. Also, the method of propagation along characteristics may not be a good choice in this system. As implemented, we must choose our time steps and longitudinal space steps to be equal. Even a stationary plane wave solution requires that the time step be small if the longitudinal structure of the solution is rapidly varying.
3. A promising technique In this section we will present a technique that promises calculations considerably. It is applicable in cases where spectral should be particularly effective for dispersion dominated systems. it to pattern forming systems, but our expectation is that it improvements. In the following we will present the application model equation with diffraction.
to improve efficiency in our methods can be applied, and As of yet we have not applied will lead to very significant of this method for a simple
3.1. A model problem Consider a partial differential u, =
equation of the form:
iu,, + g(u)
(15)
with periodic boundary conditions, and specified initial conditions U(X, 0) = U,(X). We can take advantage of the dispersion law for the linear diffraction operator to remove the fast scales from the numerics. Consider the Fourier transform of (15) tii, f -iw(k)u,
+ @(U))k
where w(k) = k’. At this point we make the key assumption that the term (g(u)jk varying. Given that assumption, it is reasonable to transform to new variables +=yk
e- io(k)t
(16) is slowly
(17)
398
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
to get jk = eiw@)‘(f(y)),
(18)
where f(y) = g(u). In practice the evaluation of the right hand side of (18) involves recovering uk from y,, taking the inverse Fourier transform (using the FFT), calculating g and taking its Fourier transform. Solving this system amounts to integrating a large system of coupled ODES most of which have highly oscillatory right hand sides. We consider that sort of problem in the next section. 3.2. Oscillatory ODES Consider the ODE j =
fP’f( t)
(19)
where we assume that f is slowly varying, but w may take any real value (not necessarily large). The function f(t) on the right-hand side may depend on t through a dependence on y and should be thought of as arising from (18). We wish to integrate this equation in such a way that we have second order accuracy for all values of w. There are many possibilities, but for specificity, we shall choose one based on the trapezoid rule. We approximate f(t) on the interval [t,, tl] by f(t) = (4 - s)f(t,) with s = (t - t,)/(tl y(tJ
-y&J
+ (+ + s)f(t,)
- to) -
= [l’f(t)
E
(tl -
(20)
l/2 for s E [ - l/2,
l/2]. Then we have
eiwt dt
t,)/l’*
((3 - S)f(to)
+ (i + S)f(tl))
eiW((f~+‘~)/2+s(t~-f~)) ds
-l/2
=: ctl
_
toj
f(tl)
+f(to)
e2iBs
eio((t,+f,)/2)
ds
2 +
ctl
_
to)
(tl)-f(to) 2
eio((to+tI)/2)
“*
s e2iOs
ds
/ -l/2
(21) where t9= o(t, - to)/2 and ,(,)=(~+;(cosB-~)) The two curves in Fig. 6 show the behavior of the real and imaginary parts of j3.
(22)
R.A. Indik, A.C. Newell/Mathematics
-0.5
and Computers in Simulation 37 (1994) 385-403
399
‘.:
-1.0.
0
IO 20 30 40 50 e
Fig. 6. Real (solid) and imaginary (dotted) parts of fi as a function of 13= w(tl - t,)/2.
Note that p is bounded for all 8 and becomes small for large o. This accurately captures the sort of decoupling that occurs when there are large frequency mismatches, without requiring those frequencies to be resolved. It is possible to prove the following error estimate for this scheme Theorem 1. Giuen that j = e’“‘f(t> with f(t) E C*([t,, tll> and f”
t1 - to = ei4(to+ll)/*) --+(toM*(@
y(tl) _y(to)
+f(G(V)
+ c(t, - to)3F min
where 8 = w(t, - to)/2 and p(e)
(y
=
+
8- y
J-(COS
))
(24)
By varying the sort of approximation that is used for the slowly varying function f, one produces a variety of schemes. Almost any of the traditional schemes for integrating ODES may be mixed with this method to produce a scheme for integrating oscillatory ODES with known frequency. Thus one can also produce explicit schemes, or schemes of higher order. 3.3. Application to PDE’s In practice, one wishes to use this method without transforming to the derived variables y,. If we rewrite equation (21) applied to the system of equation (18) in terms of uk = yke-io(k)f we get 4,)
-
=
U/&o)
e
-i4h
e
+t,bQ)
-i4kKtI-to)
t1 -
to
-+u(WW
+
(25)
e-i4’o+f1)/*)_
t1 2
to du(toM*(e)
(26)
This implicitly defines uk(t,) as the solution of a nonlinear equation. We do not wish to apply
400
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
SPATIAL-TEMPORAL
SPECTRUM - FREE RUNNING
+540
ETUNING (GW
-540
+ 4.4
-4.4 WAVE NUMBER (km-‘) BROAD AREA LASER (Index Guided)
Fig. 7. Intensity
as a function
of wavenumber
and frequency
for a free running
broad area semi-conductor
laser.
R.A. Indik, A.C. Newell /Mathematics and Computers in Simulation 37 (1994) 385-403
401
Newton-Raphson, as the equations are fully coupled, so the Jacobian matrix is dense. However, when we calculate the Jacobian (Jac) of the left-hand side of this equation we get 1
_
e-i4(t0+fl)/2)____
t1 2
to
P(e) Jack>.
(27)
Providing the step size times the Jacobian of g(u) is small, we can use the identity matrix I as a good approximation to the Jacobian of the nonlinear equation we must solve implicitly. Thus our iteration becomes:
where the term rhs, refers the the right hand side of (26). The cost of each iteration is dominated by the cost of evaluating g(u), which requires a pair of fast Fourier transforms. Since the discretization of time introduces an error of order (tl - t,)3, the iteration need not be performed beyond the point where that accuracy is obtained. As a starting point for the iteration we choose Uk(t,)
= Uk(to)
e-io(k)(fl-fd+
to
t, e-iw((ro+r~)/2)_g(u(to))(p(e)
+p*(O)).
This corresponds to using the oscillatory version of forward Euler to provide our first guess. We have tested this method for a few simple cases so far, and have been quite impressed with its effectiveness. In the case of the Nonlinear Schrodinger equation, the following table exemplifies our preliminary experience. The first column shows the step size, the second column the error in the L2 norm of the answer, and the last column, the number of evaluations of the function g, that were required. Step size
L2 error
evaluations
l/Q
0.0018 0.00066 0.00019 0.000048 0.000012 0.0000030 0.00000075 0.00000019 0.000000047
56 49 65 101 164 320 543 1036 2048
l/4 l/8 l/16 l/32 l/64 l/128 l/256 l/512
The problem was integrated using 1024 points in a domain of size 8~. One might think that in a case where one knows the nonlinear dispersion relation, even better results could be obtained. For nonlinear Schriidinger, this was not the case. This can be understood since the difference between the nonlinear, and linear dispersion relationships is bounded (in terms of the amplitude for NLS), and in our simulation was order 1. It is our hope that this method will improve our efficiency considerably. Fig. 7 illustrates how the linear dispersion relationship can dominate the behavior of a nonlinear system well beyond threshold. The figure is a plot showing intensity as a function of wavenumber and frequency in
402
R.A. Indik, A.C. Newell/Mathematics
and Computers in Simulation 37 (1994) 385-403
a simulation of a broad-area semiconductor laser. The equations that describe its behavior are not dissimilar to the Maxwell-Bloch equations. You can clearly see that the intensity is clustered around a parabolic dispersion. In order resolve that picture, we had to resolve the full range of frequencies. We expect that using this method, we will be able to resolve frequencies corresponding to the thickness of the intensity spread about that dispersion curve.
4. Conclusion In the preceding we have attempted to make clear some of the issues involved in numerical exploration of pattern forming systems in Nonlinear Optics. In order to do so we have presented our view of what a pattern is, and illustrated that point of view with a couple of specific examples we have worked with coming from Nonlinear Optics. The challenges of multiple time and space scales inherent in pattern forming systems, and the particular challenges of dispersion dominated equation in that regime determine which methods are most useful. In addition we have presented a very promising technique for improving efficiency in integration of these systems, by taking advantage of the dispersion, that we believe will be most useful.
Acknowledgements
The authors would like to thank the Arizona Center for Mathematical Sciences (ACMS) for support. ACMS is sponsored by AFOSR contract FQ8671-9000589 with the University Research Initiative Program at the University of Arizona. In addition, Robert Indik would like to acknowledge many useful discussions with J.B. Geddes, E. Meron, J.V. Moloney, W.J. Firth and N.B. Abraham.
References [l] G. Grynberg, M. Pinard and P. Verkerk, Spontaneous symmetry breaking in a ring four-wave mixing oscillator, Europhys. Lett. 9 (1989) 139. [2] G. Grynberg and L.A. Lugiato, Quantum properties of hexagonal patterns, Optics Comm. 101 (1993) 69. [3] M. Vallet, P. Verkerk and G. Grynberg, Instabilities in a ring cavity filled with a Kerr medium and submitted to a bimode field. [4] T. Honda, Hexagonal pattern formation due to counterpropagation in KNbO,, Optics Lett. 18 (1993) 598. [5] J. Pender and H. Hesselink, IEEE J. Quant. Electron 25 (1989) 395; Optics Comm. 75 (1990) 47. [6] J. Pender and H. Hesselink, Degenerate conical emissions in atomic Sodium vapor, Opt. Sot. Am. B 7 (1990) 1361. [7] N. Tan-No, H. Hoshimiya and H. Inaba, IEEE J. Quant. Electron. 16 (1980) 147. [8] F.T. Arecchi, S. Boccaletti, P.L. Ramazza and S. Residori, Transition from boundaryto bulk-controlled regimes in optical pattern formation, Phys. Rev. Lett. 70 (1993) 2277-2280. [9] R. Macdonald and H.J. Eichler, Spontaneous optical pattern formation in a nematic liquid crystal with feedback mirror, Optics Comm. 89 (1992) 289-295.
R.A. Indik, A.C. Newell /Mathematics
and Computers in Simulation 37 (1994) 385-403
403
[IO] S.R. Liu and G. Indebetouw, Periodic and chaotic spatiotemporal states in a phase-conjugate resonator using a photorefractive BaTiO, phase conjugate mirror, J. Opt. Sot. Am. B 9 (1992) 1509-1520. [ll] N.B. Abraham and W.J. Firth, Overview of transverse effects in nonlinear optical systems J. Opt. Sot. Am. B 7 (1990) 9.51-962. [12] H.M. Gibbs, Optical Bistability-Controlling Light with Light (Academic Press, Orlando, FL, 1985). [13] W.J. Firth, private communication. [14] W.J. Firth, A. Scroggie, G.S. MacDonald and L.A. Lugiato, Hexagonal patterns in optical bistability, Phys. Rev. A. 46 (1992) R3609-12. [15] R. Chang, W.J. Firth, R.A. Indik, J.V. Moloney and E.M. Wright, Three-dimensional simulations of degenerate counterpropagating beam instabilities in a nonlinear medium, Optics Comm. 88 (1992) 167. [16] J.B. Geddes, R.A. Indik, W.J. Firth and J.V. Moloney, Self focusing media like hexagons, self defocusing media prefer squares, preprint. [17] P.K. Jakobsen, J.V. Moloney, A.C. Newell and R. Indik, Space-time dynamics of wide gain section lasers Phys. Rev. A. 45 (1992) 8129-8137. [18] A.C. Newell and J.V. Moloney, Nonlinear Optics (Addison-Wesley, Redwood City, CA, 1992). [ 191 Q. Feng, J.V. Moloney and A.C. Newell, Amplitude Instabilities of Transverse Traveling Waves in Lasers, Phys. Rev. Lett. 71 (1993) 1705. [20] P.K. Jacobsen, J. Lega, Q. Feng, M. Staley, J.V. Moloney and A.C. Newell, Nonlinear transverse modes of large aspect ratio homogeneously broadened lasers: I. Analysis and numerics, Phys. Rev. A., submitted. [21] P.K. Jacobsen, J. Lega, Q. Feng, M. Staley, J.V. Moloney and A.C. Newell, Nonlinear transverse modes of large aspect ratio homogeneously broadened lasers: I. Pattern analysis near and beyond threshold, Phys. Rev. A., submitted. [22] A. Hagberg and E. Meron, Domain walls in nonequilibrium systems and the emergence of persistent patterns, Phys. Rev. E. 48 (1993) 705. [23] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Wiley, New York, 1967, 2nd Ed.). [24] F.D. Tappert, The parabolic approximation method, in: J.B. Keller and J.S. Papadakis, eds., Wave Propagation and Underwater Acoustics, Lecture Notes in Physics 70 (Springer, Berlin). [25] J.A. Fleck, J.R. Morris and M.D. Feit, Appl. Phys. 10 (1976) 129. [26] W.J. Firth, C. Penman and C. Pare, Transverse counterpropagating instabilities in Kerr slices, Optics Comm. 75 (1990) 136-140. [27] D. Kahaner, C. Moler and S. Nash, Numerical Methods and Software (Prentice Hall, Englewood Cliffs, NJ, 1989). [28] J.M. Hyman and B. Larrouturou, The numerical differentiation of discrete functions using polynomial interpolation methods, in: J.F. Thompson, ed., Numerical Grid Generation of Partial Differential Equations (Elsevier, New York, 1982) 487-506. [29] C.W. Gear, Invariants and numerical methods for ODES, Physica D 60 (1992) 303-310.