- Email: [email protected]
Polygonal approximation to the flow on the critical surface for the bistable equation
C o m p u t e r s Math. Applic. Vol. 25, No. 1, pp. 45--51, 1993 Printed in Great Britain. All rights reserved
0097-4943/93 $5.00 + 0.00 Copyright(~ 1992 Pergamon Press Ltd
P O L Y G O N A L A P P R O X I M A T I O N TO T H E F L O W O N T H E CRITICAL SURFACE FOR THE BISTABLE EQUATION V. MOLL Department of Mathematics, Tulane University New Orleans, LA 70118, U.S.A. (Received November 1991) A b s t r a c t - - W e present a method to follow the evolution of the flow on the stable manifold of the stationary wave for the bistable equation. Previous calculations indicated that a large class of initial data (those whose support above threshold is kept away from 0 and oo) on this surface, are characterized by a constant value of their area. The results presented here show that this is not sufficient to characterize data on threshold. We have observed that the areas converge smoothly under the flow to the area of the standing wave. We conjecture that the critical surface J~ is divided into three distinct regions. In the first the critical multiplier is minimal, in the second data has constant area, and the third is a region containing the standing wave. The procedure is an adaptation of Euler's polygonal method for the solution of ordinary differential equations. 1. I N T R O D U C T I O N T h e s t u d y of s y s t e m s of e q u a t i o n s of t h e f o r m ut = ux~ + f ( u , v),
(1)
v, = ~ g ( u , v)
h a s often b e e n g u i d e d b y r e s u l t s o b t a i n e d for specific choices o f f a n d g. T h e F i t z h u g h - N a g u m o system: ut = u ~ + u(1 - u ) ( u - a) - v,
v, = e(u - ~v)
(2)
t h a t first a p p e a r e d as a simplified version of H o d g k i n - H u x l e y ' s e q u a t i o n m o d e l l i n g n e r v e c o n d u c t i o n a n d its piecewise linear version ut = u~: - u + H ( u -
a)-
v,
vt = ~(u - 7 v )
(3)
i n t r o d u c e d b y M c K e a n in [1] are f r e q u e n t l y used to discover g e n e r a l results for (1). (Here, H is t h e H e a v i s i d e s t e p f u n c t i o n . I t is a s s u m e d t h a t 0 < a < 1, e > 0, 7 > 0.) A c e n t r a l p r o b l e m a s s o c i a t e d w i t h (1) is to classify t r a v e l l i n g wave solutions, a n a l y z e t h e i r s t a b i l i t y p r o p e r t i e s a n d prove t h a t stabilization o f initial data to these w a v e s occurs. Here, we c o n t i n u e t h e s t u d y i n i t i a t e d in [2] for t h e single e q u a t i o n o b t a i n e d f r o m (2) or (3) b y d r o p p i n g t h e v a r i a b l e v: ut = u ~ + u(1 - u ) ( u - a), (4) and u t = ux~ -
u + H(u
-
a).
(~)
Partially supported by the Louisiana Board of Regents LEQSF RD-A-31.
Typeset by JIA/tS-TEX 45
46
V. MOLL
Both equations admit monotone fronts, travelling waves ~o(z), z = x + kt, with ~o(-oo) = 0 and ~o(oo) = 1. In the cubic case, the speed is k = v~(½ - a) and in the piecewise linear case it is given implicitly by
~o°° ~- ~'ooQ(x, s) dz ds = a. Here, -
1
e_(t+x2/4,)
is the kernel of the linear part of equation (5). Fife and McLeod [3] proved that in the cubic case the front is globally stable. This was proven by McKean [4] in the piecewise linear case by studying the free boundary re(t) defined by u(m(t),t) = a. These equations also admit a standing wave w(x) that is symmetric about its peak and vanishes at -t-co. The infinitessimal analysis shows that this wave is unstable with a single unstable mode. Thus, this equilibrium is locally a saddle point. This has been verified in the large in [5] for the piecewise linear case. The global stable manifold A4 for the standing wave w was constructed using a shooting method along rays of the form Cuo, C E R +, where the initial data u0 is normalized by u0(0) = 1. The result is that each ray contains a unique value where convergence to the wave takes place. This is the critical multiplier. This construction was extended by Flores [6] to Nagumo's equation (4). , ~ is a smooth surface of codimension 1 in the space of initial data and it consists of all functions attracted to the standing wave w under the evolution described by (4) and (5). It separates data that collapse to 0 from those that expand to 1. In this sense, the critical manifold .M is a threshold surface. Data converging to 1 are above Iv/and are referred as super-threshold, data collapsing to 0 are below A4 and are called sub-threshold. Calculations for the critical multiplier in the piecewise linear case were described in [2]. The main finding of the numerical experiments was that critical data have the same area, provided the set where the initial data is above the threshold parameter a is kept away from 0 and oo. These results were confirmed by an asymptotic expansion of the free boundary re(t) that controls the evolution in the piecewise linear case. It was quite surprising that this common area is not the area of the standing wave. The critical multiplier attains its minimum value a for the class of data considered by Terman [7]. He has proven the existence of a constant M(a) with tl~,e following property: if uo(x) is above the threshold level a on {x G R :lxl _< M(a)}, then u0 is superthreshold. Therefore, the critical multiplier for this type of data must be a. The area for data in this class can be made arbitrarily large by considering examples that are above threshold in larger and larger domains. The main objective of this paper is to expand and clarify some conjectures stated in [2]. In particular, we will describe an Euler-type method used to approximate the evolution of the partial differential equations (4) and (5) for initial data on A4. These calculations will show that the area is not always constant and moreover that along a polygonal path on A4 one obtains convergence to the area of the wave in a smooth manner. We conjecture that the critical surface M is formed by three distinct regions that are determined by simple parameters. The first one A41 consists of the data described by Terman's result, there the critical multiplier has its minimum value a. The second region fl42 consists of the data of constant area described in [2]. Finally, the third region A~3 is a neighborhood of the standing wave. We propose that the region A42 occupies a large portion of A4, but the evolution of initial data in this region is exponentially attracted to .bt3. 2. T H E C R I T I C A L S U R F A C E The initial value problem
U¢=Uxx-i-f(u),
(6)
is solved for initial data in A'={Uo:U0(Z)=Uo(-X),Uo(Oo)=O,u~(x)
forx>O}.
Polygonal approximation
47
The nonlinear term f(u) is one of the two functions appearing in (4) or (5). The solution u(z,t) starting at u0 remains in X. In particular, there is a unique point z E R + such that u(z, t) = a; this value z = re(t) is called the median. The construction of the critical manifold A4 in [5,6] is based on a shooting method along rays of the form CUl, where ul is a normalized initial data (Ul(0) = 1). The problem (6) is solved with Uo(Z) = C u l ( z ) , where C is the shooting parameter. If C is too small, then u(z, t) converges to 0 exponentially fast. If C is too large, then u(z, t) propagates to 1 uniformly on compact sets with tails that move away like two separating fronts. These two cases comprise all positive multipliers with one exception: this is called the critical multiplier C[ul]. The normalization on ul is dropped by extending C to the ray as a homogeneous function of degree - 1,
C[ u0] =
C[u0].
(8)
Note that C[ul] >__a and equality holds for the data described by Terman. Therefore, the stable manifold of the standing wave w can be described as a level set for the functional C: (9)
A4 = {u0 e X : C[u0] = I} = C-1({I}).
The function C of initial data is smooth with nonvanishing gradient (see [5,6] for details). Therefore, A,[ is a smooth manifold of eodimension 1 in X. The expressions C[Ul] = sup{C E R + : Cul collapses to 0},
C[Ul] = inf{C e R + : Cul expands to 1}
(10)
will be used in the next section to calculate the critical multiplier. The results described in [2] have motivated the study of the set
A4A =
{
uo e X :
F
uo(z) d x = A(a)
oo
}
,
where A(a) is a constant that depends only upon a. Note that A4a is a convex surface. Moreover, as in the case of the critical surface A4 we can define the multiplier
with the property that for any u0 E X we have CA[uo] Uo E A4A. Therefore,
•A
= 6]1(1).
The results of [2] suggest that M N A4A is nontrivial. The analysis of this set will be subject of future work. 3. C A L C U L A T I O N S F O R T H E P I E C E W I S E L I N E A R C A S E Consider the initial value problem us = u~:~ -
u(x, o) = uo(x), u~:(O,t) = O,
u + g(u
- a),
(11)
ux(L,t) = O. We discretize (11) in a standard way (implicitly in the linear part and explicitly in the nonlinear term). The first boundary condition is needed for symmetry and the second one is used to simulate the behavior at co. This is used to follow the evolution and compute the critical multiplier for a fixed value of L using bisection on the scaling of the initial data. The results shown were obtained with a minimal L which does not influence the cMculations. Details are provided in [2]. C,q4~ L~sI-0
48
V. MOLL
T h e value L = 2 is sufficient in the piecewise linear case where we have used A z = 0.01. This relatively small value for L suffices to simulate the half-line problem because of the discontinuity in the nonlinear term: the evolution pushes data away from the level u = a in an exponential way. A much larger value of L is required for smooth nonlinearities. The speed of the calculations in the piecewise linear case is also improved by the use of the middle stationary solution WL to (5) on [0, L] as a comparison function. If u ( z , t) becomes larger than W L then convergence to the constant solution u = 1 is guaranteed. Similarly, if u(z, t) becomes smaller than wL, then collapse occurs. In this case, wL can be computed explicitly:
{
1 +[ ~[ ~ - ]
WL(X)
=
-
½] e ' - ~
m
tle O
,
a
where D =
+ [ -a- e- °B/ " -
m
and m is given by the quadratic equation
e m -4- e 2 L - m
(12)
e2Le - 4 m -- (1 -- 2 a ) ( e ~L -- 1)e - 2 m -- 1 = O.
Figure 1 represents the height of a critical rectangle as a function of the base computed for a = 0.20 and L = 2.00 in log-log scale. Observe that the slope of the linear part is - 1 . This calculation suggests that the area of a critical rectangle is constant provided the support is bounded away from 0 and cx~. This phenomenon was analyzed in [2] for rectangles as well as other classes of initial data. The result persists. In [2], we conjectured that for a large class of initial data the critical area is constant A ( a ) depending only upon a. We do not claim that area is the only factor determining the criticality of u0. In particular, the constant A ( a ) is not the area of the standing wave. One may conjecture that the area of critical data is always A ( a ) with the only exception being the standing wave. This is not the case. Consider the initial data: ttO(ig)
f w(x), t e2,
if I=1 > el, otherwise.
This represents a one parameter perturbation of the standing wave. The value e2 is adjusted to make u0 critical. Observe that in this case we have good apriori bounds on the location of the critical value for c2 : W(el) < e2 < w(0). Figure 2 shows the critical value of ¢2 as a function of el. The critical multiplier for the standing wave for these values of a and L is 0.24. As el ~ 0 we see convergence of the ~2 to this critical multiplier. This example shows that the wave is not a singularity with respect to the area. Height
Critical HcKean c~se. I n i t i a l 1". 0 0
.
.
.
.
epsi
Critical
d~t~t Rectangles
perturb&lions o f
.
Ion2 the wave
0+. 3 0
0.26
0.40
¢ -~ - 0 . 2 0 -C C~
Z
z
O~ 0 . 2 2
g
@
Z
Z Z
g
Z
z
Z ZZZZZZCZZZ*X
.c c~ o
-0.80
m 0.1B
IZZzI lz I
-1.40
-2.002.00.. ,
0.14
-1 I.4@
,
- 0 .i 8 0
,
_0i. 20 ,.
i , 0.40
I .00
log(basel Figure 1. Height of a critical rectangle. Piecewise linear equation, a = 0.20, L = 2.00.
0.100.00
,
01.20 ,
0 i. 4 0
,
e p s i I on
i 0.60
,
0 i. 8 0
,
.00
I
Figure 2. Critical multiplier for a perturbation of the wave.
Polygonal approximation
49
In the next section, we will present a polygonal approximation to the flow on 2~4 that will verify the absence of sharp transitions at the wave. 4. P O L Y G O N A L A P P R O X I M A T I O N
T O T H E F L O W O N A4
The scaling of u0 to the surface .M is now used to implement an Euler type method to approximate the flow on Ad. We start with any data u0 E X. The bisection method explained in Section 2 is used to compute the critical multiplier for u0 : Ul = C[u0] u0. This is the starting point of the flow. The evolution equation (11) is next solved with initial data ul and this evolution is followed for a short time At. The exact calculation would produce a curve completely contained in .M. Instead, errors in the discretization will produce a curve tangent to .£4 at ul. The bisection method can now be used again to find a scalar multiple of u(x, At) that is critical: u2(x) = C[u(x, At)] u(x, At). This forms the second point of the polygonal approximation. This procedure is iterated until convergence is observed. This procedure is the analogue of Euler's polygonal approximation method to solve initial value problems for ordinary differential equations. Note that the vertices of the polygon are always on A/[. Also we have some control on the propagation of the error. At each step the critical multiplier is computed within a prescribed accuracy. Therefore the error is controlled by the size of the gradient of C at points on the critical surface A4. Unfortunately, the expression for this gradient given in [5, Section 9] forbids any reasonable estimation. Observe that if At is sufficiently small then one has good information on the size of the next critical multiplier. This was the point where the code used in [2] spent most of its time: getting an initial approximation for the critical multiplier. This makes the polygon calculations relatively fast. We have produced a sequence of data lying on the critical surface f14. Figure 3 shows a sequence starting with rectangular initial data of base 0.30 and consisting of 20 points on the polygon. In this case one observes convergence to the standing wave w. Figure 4 shows the evolution of the critical multipliers and Figure 5 shows the corresponding areas. Recall that because of the normalization of initial data, the critical multiplier of the wave is not 1 but its value at the origin. Finally Figure 6 shows the relative difference with respect to the wave measured in the L 2 norm. Calculations were performed with other types of initial data and similar results were obtained. Data
along
the
of
r'e¢{~XngJes.
F'vo|u~lon
0.50
,
,
,
,
,
polygon ,
Hultlpliers
BaseaO.30
,
,
0.50
0.40
.~
X0 " 3 0
D
Evolution
along of
~he
rectangles.
polygon
Bame-0.30
0.40
0.30
0.20
-
0. •- 0.20
0.10
-
0.10
0.00
.-.
0.00
0.30
0.60
0.90
1 .20
o 0"00.00
1 .50
i
1.00
I
2.00
I
3.00
,
4.00
5.00
t*10x~(-31 F i g u r e 3. E v o l u t i o n o n t h e c r i t i c a l s u r f a c e s t a r t i n g from a rectangle.
5. C O M P A R I S O N
F i g u r e 4. C r i t i c M m u l t i p l i e r s a l o n g t h e p o l y g o n M a p p r o x i m a t i o n t o t h e flow.
WITH ASYMPTOTICS
IN T H E P I E C E W I S E
LINEAR CASE
The evolution of the area can be described explicitly in terms of the free boundary re(t). Integrating the equation u,
=
u~,: -
u + H(u - a)
(13)
50
V. Areas 0.50
a|ong the Evolution of r e c t a n g l e l .
MOLL
polygon
Errors &long the polygon
Ba.e=0.30
Evolution of r e c t a n g l e l .
0.20
0.40
Ba$e=0,38
0.1b
0.30
0.12 L 0 C
o
t. a) 0.08
L 0.20
CJ _J
0.10
8"000.00
0.04
2100'4100
6100
0.0 00100
8!00 t 10.00
~.~o t
2.00 '
'
3 I00
'
4.0o '
'
s.oo
t~K10~K ( - 3 )
Figure 5. Critical a r e a s a l o n g t h e p o l y g o n a l a p p r o x i m a t i o n t o t h e flow.
Figure 6. Convergenceto the wave measured in the L2 norm.
over R and introducing A ( t ) = f~_o~ u ( x , t ) dx we obtain + 2re(t).
A'(t) = -A(t)
Therefore, A ( t ) = A(O) e - t + 2 f0 t re(s) e - ( t - ' )
ds.
(14)
Let ~ = lim re(t) (this limit exists for critical initial data) and define ml (t) as re(t) = r=n + m l (t) e - '
(15)
and replace it in (14) to get
A(t) -
= [.4(0) - 2
le-' + 2 e - '
f
mx(s) d,.
(16)
This shows convergence of A(t) to the area of the wave ( = 2 ~ ) , provided that the integral term is small. The result of [2] that m l ( t ) = O(t -1) as t ---, eo it is not sufficient to guarantee this, and an improvement of this error is required. Nevertheless, we compare the approximation
.4(0 =
+ [A(0) - 2,
le-t
(17)
with the areas computed by the polygonal method . The results are shown in Figure 7. An improved estimate of the remainder term in (15) will produce a better agreement. 6. E X A M P L E
FOR NAGUMO'S
EQUATION
In this section, we briefly present corresponding results for Nagumo's equation u , = u~,~ + u ( 1 -
u)(u
-
a).
(18)
As before, the problem is discretized on [0, L] with Neumann boundary conditions at the endpoints. The calculations were performed with several types of initial data for L = 12.00 and a = 0.20. Figure 8 shows the critical height for rectangular initial data in log-log scale. Note that the slope of the linear part is - 1. Calculations with other classes of initial data also reveal the constant critical area phenomenon as in the piecewise linear case. As before, this constant area is not the area of the standing wave. Polygonal approximation to the flow on the critical surface does show however, that there is smooth convergence to the area of the wave.
Polygonal approximation
Comparison 0,.60
Critical
asymptotlcs
with
Z , computed ,, , =
x ,
51
Height
Nagumo case. I n i t i a l
asymptotic
data
Rectangle=
1.00
i
x 0.60
0,4U zzzzzz* 0.36
xx=xxxxxxxxXx
z
z
-w 0 . 2 0
zzzzzzz xx~xxXx xxX z = xxxxx
z z
~t~xxxx
,0
= z
0.24
zzzxzz z
o
•
0"000.00
zZ=Zz I
-0.60
0.12
J
0:.24 '
0:.48 '
0h.72 '
0L.96
'
1.20
t,10~*(-2} Figure 7. Comparison with the asymptotic approximation.
-1.00
-1,00
zz
_0=.68 i
-0L.36 '
-0=.04 '
0=.28 '
0.60
log(base}
Figure 8. Critical multiplier for a rectangle. Nagumo equation, a = 0.20, L = 12.00.
7. C O N C L U S I O N S We have presented a method for following the evolution of the bistable equation on the stable manifold of the standing wave w. This is an extension of the work presented in [2], which showed that a large class of critical data have constant area but that this common value is not the area of the standing wave. The present results demonstrate that there is no sharp change in the value of the area of a critical data at the wave. This eliminates the possibility of a singularity in the area at the standing wave, which was considered as a possible explanation for the discrepancy. We conclude that the critical surface .Av( is divided into three regions. In one we find data with critical multiplier equal to the threshold value a. Their area can be arbitrarily large. This type of data correspond to the superthreshold result described by Terman [7]. A second region .~v{2 consists of the data obtained in [2] with a constant area. We conjecture that this region occupies a large portion of the critical surface .~4 but the evolution described by (4) and (5) leaves this region exponentially fast. The last region ¢V(3 is a neighborhood of the standing wave. Here, data rapidly evolve into the wave. Moreover we conjecture that in this neighborhood the area has a local maximum at the standing wave. A detailed study of these possibilities is now being conducted. R E F E R E N CES 1. H.P. McKean, Nagumo's equation, Advances in Mathematics 4, 209-233 (1970). 2. V. Moll and S. Rosencraxts, Calculation of the threshold surface for nerve equations, Siam Journal oJ Applied Math. 50, 1419-1441 (1990). 3. P. Fife and B. McLeod, The approach of nonlinear diffusion equations to travelling front solutions, Arch. Rat. Mech. Anal. 65,335-361 (1977). 4. H.P. McKean, Stabilization of solutions of a caricature of the Fitzhngh-Nagumo equation, Comm. Pure Applied Math. 36, 291-324 (1983); 37, 299-301 (1984). 5. H.P. McKean and V. Moll, Stabilization to the standing wave in a simple caricature of the nerve equation, Comm. Pure Applied Math. 39, 485-529 (1986). 6. J. Flores, On a threshold of codimension 1 for the Nagumo equation, Comm. PDE 13, 1235-1263 (1988). 7. D. Terman, Threshold phenomena for a reaction-diffusion system, Jour. Di~. Equations 47,406-443 (1983). 8. J. Smoller, Shock Waves and Reactlon-Di.~uslon Equations, Springer-Verlag, New York, (1982).
0097-4943/93 $5.00 + 0.00 Copyright(~ 1992 Pergamon Press Ltd
P O L Y G O N A L A P P R O X I M A T I O N TO T H E F L O W O N T H E CRITICAL SURFACE FOR THE BISTABLE EQUATION V. MOLL Department of Mathematics, Tulane University New Orleans, LA 70118, U.S.A. (Received November 1991) A b s t r a c t - - W e present a method to follow the evolution of the flow on the stable manifold of the stationary wave for the bistable equation. Previous calculations indicated that a large class of initial data (those whose support above threshold is kept away from 0 and oo) on this surface, are characterized by a constant value of their area. The results presented here show that this is not sufficient to characterize data on threshold. We have observed that the areas converge smoothly under the flow to the area of the standing wave. We conjecture that the critical surface J~ is divided into three distinct regions. In the first the critical multiplier is minimal, in the second data has constant area, and the third is a region containing the standing wave. The procedure is an adaptation of Euler's polygonal method for the solution of ordinary differential equations. 1. I N T R O D U C T I O N T h e s t u d y of s y s t e m s of e q u a t i o n s of t h e f o r m ut = ux~ + f ( u , v),
(1)
v, = ~ g ( u , v)
h a s often b e e n g u i d e d b y r e s u l t s o b t a i n e d for specific choices o f f a n d g. T h e F i t z h u g h - N a g u m o system: ut = u ~ + u(1 - u ) ( u - a) - v,
v, = e(u - ~v)
(2)
t h a t first a p p e a r e d as a simplified version of H o d g k i n - H u x l e y ' s e q u a t i o n m o d e l l i n g n e r v e c o n d u c t i o n a n d its piecewise linear version ut = u~: - u + H ( u -
a)-
v,
vt = ~(u - 7 v )
(3)
i n t r o d u c e d b y M c K e a n in [1] are f r e q u e n t l y used to discover g e n e r a l results for (1). (Here, H is t h e H e a v i s i d e s t e p f u n c t i o n . I t is a s s u m e d t h a t 0 < a < 1, e > 0, 7 > 0.) A c e n t r a l p r o b l e m a s s o c i a t e d w i t h (1) is to classify t r a v e l l i n g wave solutions, a n a l y z e t h e i r s t a b i l i t y p r o p e r t i e s a n d prove t h a t stabilization o f initial data to these w a v e s occurs. Here, we c o n t i n u e t h e s t u d y i n i t i a t e d in [2] for t h e single e q u a t i o n o b t a i n e d f r o m (2) or (3) b y d r o p p i n g t h e v a r i a b l e v: ut = u ~ + u(1 - u ) ( u - a), (4) and u t = ux~ -
u + H(u
-
a).
(~)
Partially supported by the Louisiana Board of Regents LEQSF RD-A-31.
Typeset by JIA/tS-TEX 45
46
V. MOLL
Both equations admit monotone fronts, travelling waves ~o(z), z = x + kt, with ~o(-oo) = 0 and ~o(oo) = 1. In the cubic case, the speed is k = v~(½ - a) and in the piecewise linear case it is given implicitly by
~o°° ~- ~'ooQ(x, s) dz ds = a. Here, -
1
e_(t+x2/4,)
is the kernel of the linear part of equation (5). Fife and McLeod [3] proved that in the cubic case the front is globally stable. This was proven by McKean [4] in the piecewise linear case by studying the free boundary re(t) defined by u(m(t),t) = a. These equations also admit a standing wave w(x) that is symmetric about its peak and vanishes at -t-co. The infinitessimal analysis shows that this wave is unstable with a single unstable mode. Thus, this equilibrium is locally a saddle point. This has been verified in the large in [5] for the piecewise linear case. The global stable manifold A4 for the standing wave w was constructed using a shooting method along rays of the form Cuo, C E R +, where the initial data u0 is normalized by u0(0) = 1. The result is that each ray contains a unique value where convergence to the wave takes place. This is the critical multiplier. This construction was extended by Flores [6] to Nagumo's equation (4). , ~ is a smooth surface of codimension 1 in the space of initial data and it consists of all functions attracted to the standing wave w under the evolution described by (4) and (5). It separates data that collapse to 0 from those that expand to 1. In this sense, the critical manifold .M is a threshold surface. Data converging to 1 are above Iv/and are referred as super-threshold, data collapsing to 0 are below A4 and are called sub-threshold. Calculations for the critical multiplier in the piecewise linear case were described in [2]. The main finding of the numerical experiments was that critical data have the same area, provided the set where the initial data is above the threshold parameter a is kept away from 0 and oo. These results were confirmed by an asymptotic expansion of the free boundary re(t) that controls the evolution in the piecewise linear case. It was quite surprising that this common area is not the area of the standing wave. The critical multiplier attains its minimum value a for the class of data considered by Terman [7]. He has proven the existence of a constant M(a) with tl~,e following property: if uo(x) is above the threshold level a on {x G R :lxl _< M(a)}, then u0 is superthreshold. Therefore, the critical multiplier for this type of data must be a. The area for data in this class can be made arbitrarily large by considering examples that are above threshold in larger and larger domains. The main objective of this paper is to expand and clarify some conjectures stated in [2]. In particular, we will describe an Euler-type method used to approximate the evolution of the partial differential equations (4) and (5) for initial data on A4. These calculations will show that the area is not always constant and moreover that along a polygonal path on A4 one obtains convergence to the area of the wave in a smooth manner. We conjecture that the critical surface M is formed by three distinct regions that are determined by simple parameters. The first one A41 consists of the data described by Terman's result, there the critical multiplier has its minimum value a. The second region fl42 consists of the data of constant area described in [2]. Finally, the third region A~3 is a neighborhood of the standing wave. We propose that the region A42 occupies a large portion of A4, but the evolution of initial data in this region is exponentially attracted to .bt3. 2. T H E C R I T I C A L S U R F A C E The initial value problem
U¢=Uxx-i-f(u),
(6)
is solved for initial data in A'={Uo:U0(Z)=Uo(-X),Uo(Oo)=O,u~(x)
forx>O}.
Polygonal approximation
47
The nonlinear term f(u) is one of the two functions appearing in (4) or (5). The solution u(z,t) starting at u0 remains in X. In particular, there is a unique point z E R + such that u(z, t) = a; this value z = re(t) is called the median. The construction of the critical manifold A4 in [5,6] is based on a shooting method along rays of the form CUl, where ul is a normalized initial data (Ul(0) = 1). The problem (6) is solved with Uo(Z) = C u l ( z ) , where C is the shooting parameter. If C is too small, then u(z, t) converges to 0 exponentially fast. If C is too large, then u(z, t) propagates to 1 uniformly on compact sets with tails that move away like two separating fronts. These two cases comprise all positive multipliers with one exception: this is called the critical multiplier C[ul]. The normalization on ul is dropped by extending C to the ray as a homogeneous function of degree - 1,
C[ u0] =
C[u0].
(8)
Note that C[ul] >__a and equality holds for the data described by Terman. Therefore, the stable manifold of the standing wave w can be described as a level set for the functional C: (9)
A4 = {u0 e X : C[u0] = I} = C-1({I}).
The function C of initial data is smooth with nonvanishing gradient (see [5,6] for details). Therefore, A,[ is a smooth manifold of eodimension 1 in X. The expressions C[Ul] = sup{C E R + : Cul collapses to 0},
C[Ul] = inf{C e R + : Cul expands to 1}
(10)
will be used in the next section to calculate the critical multiplier. The results described in [2] have motivated the study of the set
A4A =
{
uo e X :
F
uo(z) d x = A(a)
oo
}
,
where A(a) is a constant that depends only upon a. Note that A4a is a convex surface. Moreover, as in the case of the critical surface A4 we can define the multiplier
with the property that for any u0 E X we have CA[uo] Uo E A4A. Therefore,
•A
= 6]1(1).
The results of [2] suggest that M N A4A is nontrivial. The analysis of this set will be subject of future work. 3. C A L C U L A T I O N S F O R T H E P I E C E W I S E L I N E A R C A S E Consider the initial value problem us = u~:~ -
u(x, o) = uo(x), u~:(O,t) = O,
u + g(u
- a),
(11)
ux(L,t) = O. We discretize (11) in a standard way (implicitly in the linear part and explicitly in the nonlinear term). The first boundary condition is needed for symmetry and the second one is used to simulate the behavior at co. This is used to follow the evolution and compute the critical multiplier for a fixed value of L using bisection on the scaling of the initial data. The results shown were obtained with a minimal L which does not influence the cMculations. Details are provided in [2]. C,q4~ L~sI-0
48
V. MOLL
T h e value L = 2 is sufficient in the piecewise linear case where we have used A z = 0.01. This relatively small value for L suffices to simulate the half-line problem because of the discontinuity in the nonlinear term: the evolution pushes data away from the level u = a in an exponential way. A much larger value of L is required for smooth nonlinearities. The speed of the calculations in the piecewise linear case is also improved by the use of the middle stationary solution WL to (5) on [0, L] as a comparison function. If u ( z , t) becomes larger than W L then convergence to the constant solution u = 1 is guaranteed. Similarly, if u(z, t) becomes smaller than wL, then collapse occurs. In this case, wL can be computed explicitly:
{
1 +[ ~[ ~ - ]
WL(X)
=
-
½] e ' - ~
m
tle O
,
a
where D =
+ [ -a- e- °B/ " -
m
and m is given by the quadratic equation
e m -4- e 2 L - m
(12)
e2Le - 4 m -- (1 -- 2 a ) ( e ~L -- 1)e - 2 m -- 1 = O.
Figure 1 represents the height of a critical rectangle as a function of the base computed for a = 0.20 and L = 2.00 in log-log scale. Observe that the slope of the linear part is - 1 . This calculation suggests that the area of a critical rectangle is constant provided the support is bounded away from 0 and cx~. This phenomenon was analyzed in [2] for rectangles as well as other classes of initial data. The result persists. In [2], we conjectured that for a large class of initial data the critical area is constant A ( a ) depending only upon a. We do not claim that area is the only factor determining the criticality of u0. In particular, the constant A ( a ) is not the area of the standing wave. One may conjecture that the area of critical data is always A ( a ) with the only exception being the standing wave. This is not the case. Consider the initial data: ttO(ig)
f w(x), t e2,
if I=1 > el, otherwise.
This represents a one parameter perturbation of the standing wave. The value e2 is adjusted to make u0 critical. Observe that in this case we have good apriori bounds on the location of the critical value for c2 : W(el) < e2 < w(0). Figure 2 shows the critical value of ¢2 as a function of el. The critical multiplier for the standing wave for these values of a and L is 0.24. As el ~ 0 we see convergence of the ~2 to this critical multiplier. This example shows that the wave is not a singularity with respect to the area. Height
Critical HcKean c~se. I n i t i a l 1". 0 0
.
.
.
.
epsi
Critical
d~t~t Rectangles
perturb&lions o f
.
Ion2 the wave
0+. 3 0
0.26
0.40
¢ -~ - 0 . 2 0 -C C~
Z
z
O~ 0 . 2 2
g
@
Z
Z Z
g
Z
z
Z ZZZZZZCZZZ*X
.c c~ o
-0.80
m 0.1B
IZZzI lz I
-1.40
-2.002.00.. ,
0.14
-1 I.4@
,
- 0 .i 8 0
,
_0i. 20 ,.
i , 0.40
I .00
log(basel Figure 1. Height of a critical rectangle. Piecewise linear equation, a = 0.20, L = 2.00.
0.100.00
,
01.20 ,
0 i. 4 0
,
e p s i I on
i 0.60
,
0 i. 8 0
,
.00
I
Figure 2. Critical multiplier for a perturbation of the wave.
Polygonal approximation
49
In the next section, we will present a polygonal approximation to the flow on 2~4 that will verify the absence of sharp transitions at the wave. 4. P O L Y G O N A L A P P R O X I M A T I O N
T O T H E F L O W O N A4
The scaling of u0 to the surface .M is now used to implement an Euler type method to approximate the flow on Ad. We start with any data u0 E X. The bisection method explained in Section 2 is used to compute the critical multiplier for u0 : Ul = C[u0] u0. This is the starting point of the flow. The evolution equation (11) is next solved with initial data ul and this evolution is followed for a short time At. The exact calculation would produce a curve completely contained in .M. Instead, errors in the discretization will produce a curve tangent to .£4 at ul. The bisection method can now be used again to find a scalar multiple of u(x, At) that is critical: u2(x) = C[u(x, At)] u(x, At). This forms the second point of the polygonal approximation. This procedure is iterated until convergence is observed. This procedure is the analogue of Euler's polygonal approximation method to solve initial value problems for ordinary differential equations. Note that the vertices of the polygon are always on A/[. Also we have some control on the propagation of the error. At each step the critical multiplier is computed within a prescribed accuracy. Therefore the error is controlled by the size of the gradient of C at points on the critical surface A4. Unfortunately, the expression for this gradient given in [5, Section 9] forbids any reasonable estimation. Observe that if At is sufficiently small then one has good information on the size of the next critical multiplier. This was the point where the code used in [2] spent most of its time: getting an initial approximation for the critical multiplier. This makes the polygon calculations relatively fast. We have produced a sequence of data lying on the critical surface f14. Figure 3 shows a sequence starting with rectangular initial data of base 0.30 and consisting of 20 points on the polygon. In this case one observes convergence to the standing wave w. Figure 4 shows the evolution of the critical multipliers and Figure 5 shows the corresponding areas. Recall that because of the normalization of initial data, the critical multiplier of the wave is not 1 but its value at the origin. Finally Figure 6 shows the relative difference with respect to the wave measured in the L 2 norm. Calculations were performed with other types of initial data and similar results were obtained. Data
along
the
of
r'e¢{~XngJes.
F'vo|u~lon
0.50
,
,
,
,
,
polygon ,
Hultlpliers
BaseaO.30
,
,
0.50
0.40
.~
X0 " 3 0
D
Evolution
along of
~he
rectangles.
polygon
Bame-0.30
0.40
0.30
0.20
-
0. •- 0.20
0.10
-
0.10
0.00
.-.
0.00
0.30
0.60
0.90
1 .20
o 0"00.00
1 .50
i
1.00
I
2.00
I
3.00
,
4.00
5.00
t*10x~(-31 F i g u r e 3. E v o l u t i o n o n t h e c r i t i c a l s u r f a c e s t a r t i n g from a rectangle.
5. C O M P A R I S O N
F i g u r e 4. C r i t i c M m u l t i p l i e r s a l o n g t h e p o l y g o n M a p p r o x i m a t i o n t o t h e flow.
WITH ASYMPTOTICS
IN T H E P I E C E W I S E
LINEAR CASE
The evolution of the area can be described explicitly in terms of the free boundary re(t). Integrating the equation u,
=
u~,: -
u + H(u - a)
(13)
50
V. Areas 0.50
a|ong the Evolution of r e c t a n g l e l .
MOLL
polygon
Errors &long the polygon
Ba.e=0.30
Evolution of r e c t a n g l e l .
0.20
0.40
Ba$e=0,38
0.1b
0.30
0.12 L 0 C
o
t. a) 0.08
L 0.20
CJ _J
0.10
8"000.00
0.04
2100'4100
6100
0.0 00100
8!00 t 10.00
~.~o t
2.00 '
'
3 I00
'
4.0o '
'
s.oo
t~K10~K ( - 3 )
Figure 5. Critical a r e a s a l o n g t h e p o l y g o n a l a p p r o x i m a t i o n t o t h e flow.
Figure 6. Convergenceto the wave measured in the L2 norm.
over R and introducing A ( t ) = f~_o~ u ( x , t ) dx we obtain + 2re(t).
A'(t) = -A(t)
Therefore, A ( t ) = A(O) e - t + 2 f0 t re(s) e - ( t - ' )
ds.
(14)
Let ~ = lim re(t) (this limit exists for critical initial data) and define ml (t) as re(t) = r=n + m l (t) e - '
(15)
and replace it in (14) to get
A(t) -
= [.4(0) - 2
le-' + 2 e - '
f
mx(s) d,.
(16)
This shows convergence of A(t) to the area of the wave ( = 2 ~ ) , provided that the integral term is small. The result of [2] that m l ( t ) = O(t -1) as t ---, eo it is not sufficient to guarantee this, and an improvement of this error is required. Nevertheless, we compare the approximation
.4(0 =
+ [A(0) - 2,
le-t
(17)
with the areas computed by the polygonal method . The results are shown in Figure 7. An improved estimate of the remainder term in (15) will produce a better agreement. 6. E X A M P L E
FOR NAGUMO'S
EQUATION
In this section, we briefly present corresponding results for Nagumo's equation u , = u~,~ + u ( 1 -
u)(u
-
a).
(18)
As before, the problem is discretized on [0, L] with Neumann boundary conditions at the endpoints. The calculations were performed with several types of initial data for L = 12.00 and a = 0.20. Figure 8 shows the critical height for rectangular initial data in log-log scale. Note that the slope of the linear part is - 1. Calculations with other classes of initial data also reveal the constant critical area phenomenon as in the piecewise linear case. As before, this constant area is not the area of the standing wave. Polygonal approximation to the flow on the critical surface does show however, that there is smooth convergence to the area of the wave.
Polygonal approximation
Comparison 0,.60
Critical
asymptotlcs
with
Z , computed ,, , =
x ,
51
Height
Nagumo case. I n i t i a l
asymptotic
data
Rectangle=
1.00
i
x 0.60
0,4U zzzzzz* 0.36
xx=xxxxxxxxXx
z
z
-w 0 . 2 0
zzzzzzz xx~xxXx xxX z = xxxxx
z z
~t~xxxx
,0
= z
0.24
zzzxzz z
o
•
0"000.00
zZ=Zz I
-0.60
0.12
J
0:.24 '
0:.48 '
0h.72 '
0L.96
'
1.20
t,10~*(-2} Figure 7. Comparison with the asymptotic approximation.
-1.00
-1,00
zz
_0=.68 i
-0L.36 '
-0=.04 '
0=.28 '
0.60
log(base}
Figure 8. Critical multiplier for a rectangle. Nagumo equation, a = 0.20, L = 12.00.
7. C O N C L U S I O N S We have presented a method for following the evolution of the bistable equation on the stable manifold of the standing wave w. This is an extension of the work presented in [2], which showed that a large class of critical data have constant area but that this common value is not the area of the standing wave. The present results demonstrate that there is no sharp change in the value of the area of a critical data at the wave. This eliminates the possibility of a singularity in the area at the standing wave, which was considered as a possible explanation for the discrepancy. We conclude that the critical surface .Av( is divided into three regions. In one we find data with critical multiplier equal to the threshold value a. Their area can be arbitrarily large. This type of data correspond to the superthreshold result described by Terman [7]. A second region .~v{2 consists of the data obtained in [2] with a constant area. We conjecture that this region occupies a large portion of the critical surface .~4 but the evolution described by (4) and (5) leaves this region exponentially fast. The last region ¢V(3 is a neighborhood of the standing wave. Here, data rapidly evolve into the wave. Moreover we conjecture that in this neighborhood the area has a local maximum at the standing wave. A detailed study of these possibilities is now being conducted. R E F E R E N CES 1. H.P. McKean, Nagumo's equation, Advances in Mathematics 4, 209-233 (1970). 2. V. Moll and S. Rosencraxts, Calculation of the threshold surface for nerve equations, Siam Journal oJ Applied Math. 50, 1419-1441 (1990). 3. P. Fife and B. McLeod, The approach of nonlinear diffusion equations to travelling front solutions, Arch. Rat. Mech. Anal. 65,335-361 (1977). 4. H.P. McKean, Stabilization of solutions of a caricature of the Fitzhngh-Nagumo equation, Comm. Pure Applied Math. 36, 291-324 (1983); 37, 299-301 (1984). 5. H.P. McKean and V. Moll, Stabilization to the standing wave in a simple caricature of the nerve equation, Comm. Pure Applied Math. 39, 485-529 (1986). 6. J. Flores, On a threshold of codimension 1 for the Nagumo equation, Comm. PDE 13, 1235-1263 (1988). 7. D. Terman, Threshold phenomena for a reaction-diffusion system, Jour. Di~. Equations 47,406-443 (1983). 8. J. Smoller, Shock Waves and Reactlon-Di.~uslon Equations, Springer-Verlag, New York, (1982).