1
(2.38)
where P2 > Po. By (2.36), (2.37) it suffices that P2 2 ( i - 1)p 1 P22 > 2i(2i + 3) + 2i(2i +.1)(2i + 3) '
i >_ 1,
(2.39)
P2 P22>- 2 i ( 2 i + 3 )
i_>l.
(2.40)
and 2 ( i - 1)(1 + p , ) + 2i(2i+1)(2i+3)'
Clearly (2.40) implies (2.39). In (2.40) the strongest inequalities are for i = 1 and i = 2: 1 p2 >-- 7-6'
P2
022 >- ~
+
1 + p, 7------0-
(2.41)
Thus we set
Pz=max
(11 (1 10' 56 +
~
+
)j2)
"
(2.42)
For Pl as in (2.25) P2 < 0.158.
(2.43)
With yp and ypp approximated by finite sums YPN, YPPN, the corresponding remainders EPs and EPP s
D. Michelson / Elementary particles as solutions of the Sivashinsky equation
511
are estimated by IEPNI, IEPPNI
(2.44)
The first- and second-order derivatives of yp and ypp are approximated by the corresponding derivatives of YPN and YPPN" The respective remainders are denoted by EP/~, EP~, EPP~, EPP~. For example, in EPP~ similar to (2.31) we obtain (2i - 1)(2i - 2)lapPil < p ~ - 2 -t- (1 + p l ) p ~ - 3
(2.45)
and hence IEPP~I
~pN-lr2N-I ( 1 + l + p l ) ( P2
1 --p2 rE) -~ < 7 . 6 × 10 -15 .
(2.46)
It is easy to see that the above bound dominates also IEP~I and hence the first derivatives IEPP~I and IEP~I. Altogether we set + 10 -14 as the error term ERROR in the subroutine SERIES to be added to all results YN,YPN,YPPN and their derivatives up to order 2.
3. The asymptotic expansion for large r
Let us write the following ansatz:
y= ~_~r-nYn(r),
ank(P)e ikr, p = l o g r .
Yn(r)=
n=l
(3.1)
k= - n
If we denote by _z:(y) the non-linear operator - f ( y ) = ( D 2 + 2 D + 1)(D + 2 ) y + y 2
(3.2)
then the result of the application of . f to y in (3.1) will be
b.k(O) eikr,
- : ( Y ) = E r-~ n=l
(3.3)
k= - n
where the coefficients
bnk(p)
have the following form
bnk(p ) = ik( - k 2 + 1)ank(p ) -- kE(aD, - 3n + 7)a~_l.k(p) + ( D p - n + 3)an_l.k(p) + ik[aD 2 - ( 6 n - 17)Dp + an 2 - 17n + 22]a,_2.k(p)
+(Dp-n+2)(Dp-n+a)(Dp-n+5)a~-a.k(P)+
~-,
E
anlklan2k2"
(3.4)
nl +n2ffin kl +k2=k
Here Dp is the differentiation with respect to p. It is implicitly assumed in (3.4) that ank --~ 0 for Ikl > 1.
512
D. Michelson /Elementary particles as solutions o f the Sivashinsky equation
For a real function y it is also necessary that an _ k = a . . k .
(3.5)
This in turn causes bn.-k = b,.k"
(3.6)
In o r d e r for y to be a (formal) solution of (1.5) one should equate all coefficients bnk to zero. Note that for n = 1 (and hence k = 0, _+ 1), bnk in (3.4) vanish automatically. For n = 2 and k = 0, + 1 we obtain the system of O D E s all
=a,oa,,,
a{o=
-(a,o +a2o + 2a, _,a,1),
a, _, =al,.
(3.7)
Since a m is real, the a r g u m e n t of a of aj~, all=
Jail I e i'~,
(3.8)
is a free p a r a m e t e r while am, 1all I satisfy the real system (1.10). In general, for k = 0, 1 the coefficients an ~.k are defined by the system of O D E s (D o - n
+ 3)an_L0 + 2aman_l, 0 + 4Re(all,n_1,1)
-- (2Dp - 2n + 4 ) a n_ 1,1 +
=fn0,
2aman-l.l + 2al,an-1,0 = f n l ,
(3.9)
w h e r e fnk, k = 0, 1 satisfy fnk= -ik[ 3D2-
(6n-
17)Dp + 3n 2 - 17n + 22]an
-(Dp-n+2)(Dp-n+3)(Do-n+5)an-3,k-
2.k
E
~-,
anlk,an2k 2
(3.10)
n t+nz=n k l+k2=k
and the sum does not include t e r m s with an_l, k, [kf < 1. T h e coefficients ank for Ik I > 1 are d e t e r m i n e d by algebraic equations in terms of coefficients an~k, and their derivatives with nl < n. First consider system (1.10). It has two critical points: a m = la~l [ - - 0 and a m = - 1 , ]al~ I = 0. T h e first is a parabolic one with a stable manifold la~ll = 0. T h e second critical point is a hyperbolic one with the unstable manifold la~l I = 0 and a stable manifold Mst. Five trajectories of the system and the manifold Mst are plotted by c o m p u t e r in fig. 2. O u r c o m p u t a t i o n s show that Mst intersects the line al0 = 0 at laHI -~ 0.824. Since for a m > 0, lall[' > 0 and a~o < 0, the line Mst in b a c k w a r d direction tends to laHI = 0, a m = +oo. It is evident from the plot that a~o < 0 along Mst also when a m < 0 . This fact could be proved analytically. Indeed, d e n o t e by ~ the d o m a i n 21a~112 < - a m - a2o, w h e r e a[~ >__0. N e a r the critical point a m = - 1 , [a~l = 0 the asymptotics of Mst (e.g. see (3.22)) shows that Mst lies outside 11. In o r d e r to e n t e r f~ in backward direction the trajectory should cross the b o u n d a r y of f~ with a positive slope dam/dla~l 1. But this is impossible since at the b o u n d a r y of 11 this slope is 0. D e n o t e by 2 0 the d o m a i n b o u n d e d by Mst and the axis [al~ I = 0. Clearly each point in -~0 is attracted by the flow to the origin a m = [a~ I = 0 while the points outside of 2 0 tend in forward direction to [a~[ = 0, a m = -oo. It is also clear that the trajectories of the flow n e v e r cross the line al0 = 0 from below since at a m = 0, a~0 < 0.
513
D. Michelson / Elementary particles as solutions of the Sivashinsky equation
Next let us find the asymptotics of Mst near the critical point al0 = U~ ---- laHI, u 2 ---- 1 + al0 the graph of Mst is described by a power series,
--1, lalll = 0 . In variables
] ~ a n U '~.
u2=g(ul)=
(3.11)
n>2
Differentiation of (3.11) with respect to p and substitution of u[, u~ from (1.10) leads to the identity
u 2- u 2 - 2u 2 = ( ~ nanu~ -1)ul(u 2 - 1 ) .
(3.12)
~n>2
Hence u2+ ~nanU~=U2(U2+ n>2
Y'.nanu~)+2u 2
(3.13)
n>2
or in terms of u 1 only
Y'. (n+ 1)anu~=2u2+( ~.anu~) ~ (n+ 1)anu ~. n~2
~n~2
(3.14)
In~2
The coefficients a n are now defined by the recurrent relation n-2
( n + l ) a n = ~., (m+l)arnan-m+262,n.
(3.15)
mffi2
In particular a 2 = 2/3, a 3 = 0, a 4 ffi 4/15. Since all a , _> 0, the partial satisfy the inequality
sums
1 2 SN( UO < 2u 2 + "$SN( Ul).
SN(Ul) =
Y~n=Nl(n + 1)anU,~ (3.16)
The corresponding equality has positive roots ½(3 + gt9 - 24u 2 ) provided u 1 _<3 / 2 ~
= 01 --- 0.612.
(3.17)
Since for small ul
SN(Ul) < ½(3--
gt9 - 24u 2 ),
(3.18)
the same inequality is correct for all u~ as in (3.17) and all N. I n particular the radius of convergence of the series in (3.11) is at least 01. The curve Mst in fig. 2 was actually computed using the above series for u I < 0.3 and then continuing with a difference approximation of system (1.10). We will need in the future (see fig. 3) the bound on Mst for lalll < 0.23. By (3.18) O
3 2 ( n + 1)p]'
(3.19)
D. Michelson / Elementary particles as solutions of the Sivashinsky equation
514
and hence, for u 1 _
~_~anU'~ <
(1-uJp,)
32 UI=--~(I--Ul//Pl)
1 4
1
- U 4.
(3.20)
n>_4
In particular,
anU'~ <_3.42u~
for u 1 _< 0.23
(3.21)
n>_4
and hence along Mst
2 12 + 3.51all r4 - 1 + ~1al1122 < a l 0 _< - 1 + 31all
for lall I _<0.23.
(3.22)
The above estimate will be used once in section 4 to prove that the lowest rectangle in fig. 3 lies in the convergence domain. Now we turn to system (3.9) with n ___3. Here we should make the following important observation.
Lemma 3.1. If ank , Ikl < n are solutions of the equations bnk = 0 and are bounded as p--> ~, then necessarily ank=in-leik'~nk,
tink is real, a i s f i x e d
(3.23)
and ano=O
for e v e n n .
(3.24)
Proof Indeed, substitution of (3.23) into (3.4) leads to bnk = in eik'~ t)nk
(3.25)
where bnk likewise bnk depends on the coefficients ~i but with real coefficients. Hence the equations Dnk = 0, which define rink, are real. Since tilk, Ikl < 1, by definition are real, so will be all rink. It is clear from the definition of bnk and from (3.5) that bn.-k = bnk-
(3.26)
Hence, by (3.25) bn,-k = ( - 1) n/;n,k
(3.27)
bn, 0 = 0
(3.28)
and for odd n.
The same argument applied to fnk in (3.10) shows that fn,0 = 0
for odd n.
(3.29)
515
D. Michekon /Elementary particles as solutions of the Sivashinsky equation
= 0 and hence u,,_~,~ satisfies the homogeneous
By (3.23) for odd n, Re(a,,Z,_,,,)
equation (3.30)
[DP-n+3+2aio(p)]a,_l,o=0. For n = 3 we obtain
[D,+ 24ob)]~2,0=0.
(3.31)
If one compares (3.31) with the first equation in (1.10) it becomes clear that u2,01u,112 is constant and hence u2,0 - lull1 - ’ + m as p + m. Hence the only bounded solution of (3.31) is u2,0 = 0. For n > 3 the fact that u,,(p) + 0 as p + m is enough to deduce (3.24). Notice that there is a gap in our arguments. Namely, one should prove that system (3.9) for n ~3 does not have non-trivial bounded homogeneous solutions. For n > 3 it is obvious, because a,, and ai, tend to 0 as p + m so that system (3.9) for large p becomes an infinitesimal perturbation of the constant coefficient system (D, -n + 3)~1.0
=fno,
(D, -n + 2)%-i,,
(3.32)
=fn1.
For II = 3 we have a system
q&207
a211
%,-l)T=~ob)
(3.33)
*(%0t%17%.-l)T,
where -a10 4m
=
a11
i-
a11
-41
-a11
1+a,, 0
0 I
+a10
(3.34)
. I
Obviously, Re A,(p) = diag( -a,,, 1 + a,,, 1 + alo) is positive for large p, say p > po. Multiplying eq. (3.33) by the vector (uZo, u2, i, a,, _,> and integrating, leads to [ho12
+ 21%,lz],PO=
/“[
PO
-~oIQ~oI~
+ 2(1 +qo)la2J2]
dp > 0.
(3.35)
Hence, the boundedness of u,~, u2, implies the uniform boundedness of the integral at the right-hand side of (3.35) so that u2i E L2(po,m) and u2i + 0 as p 3 m. This in turn implies that Ju,,(p)l 2 ~u,,(p,)1 for p > p. and sufficiently large. However, the integral 01
/
-u,,dp PO
= ,/
$$dp
= ~~~loglo,,(Po)/a,,(p)
is diverging, which contradicts the boundedness zero.
I= m
(3.36)
of the integral in (3.35). Thus u2k, Ikl I 1 are identically
516
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
Now one can efficiently estimate an_Lk, [kl < 1. For n odd
~,n - l , 1 _ ( n _ 2 + a l o ) 6 , _ l . l
= - ~f~l 1"
(3.37)
where f,k is defined likewise /).k in (3.25) and is real. Hence
IL,(p)I sup 16,_L~(p)l < sup 2 [ n - 2 + a , 0 ( p ) ] " P > PO
(3.38)
P >-PO
For n = 3, this estimate is not good enough since the denominator 1 + alo(p) in the domain 5~o could be arbitrarily close to 0. Instead we rewrite eq. (3.37) for the variable c2~ = a2~/[a~l:
c~,la,l [ + c211a,,[alo- (1 + alo)C211a,l I =
(3.39)
--1/31.
By (3.10) f31 = - i ( 3 D ~ - D p - 2 ) a l ) - 2 a 2 2 a , _ l 1"
2
= - i ( 3 a [ l a l o + 3 a l i a [ o - a l l a l o - 2all) + ~lallal,-i 1 = - i a l l ( - 3 a 2 o + 3a10 + 3a~o + 61aur 2 + a l o + 2 + 31aul e)
(3.40)
= iall(2 + 4a,o + ~ l a u 1 2 ) . Thus c21 satisfies the equation
(3.41)
c~1 - czl = (1 + 2alo + ~ laul2), with the explicit solution
(3.42)
c21(Po) = - f~e",,-" (1 + 2a,0 + ~[all[2) dp. Po
Obviously, (3.43)
Ic21(Po)1 < max 11 + 2alo(p) + ~ l a l , ( p ) f 2 [ P ~ Po
and finally (3.44)
la2~(po)[ _< laH(Po)l max I1 + 2ajo(p) + ~lal~(p)121 • P > PO
The last estimate does not require the negativeness of alo(Po) and applies to any Po, provided lall(po)l, alo(Po) lie in -~o. For even n we rewrite (3.9) as
2an-l,o,an-l,l,an-l.-1
= ½(go, - g , , - g , - , )
-A(p)(an-l.o,a.-l,l,an-1
-1) T (3.45)
D. Michelson / Elementaryparticlesas solutions of the Sivashinsky equation
517
where
A ( p ) = A o ( p ) + diag(½(n - 3 ) , n - 3, n - 3)
(3.46)
and Ao(p) is defined in (3.34). Again, R e A(p) is diagonal: R e A ( p ) = diag(½(n - 3) - alo, n - 2 + alo, n - 2 + alo ). Now multiply eq. (3.45) by e-2A(P-Po)(a,_l,o,a,_l.l,a,_l,_l) infinity. W e obtain
(3.47) with A > 0
and integrate from Po to
¼ ( l a . _ l , o ( P o ) 12 + 4 1 a . _ l , l ( p o ) I z) oo
+
[½(n
-
3 - 2alo - A ) l a . _ l , o l 2 + 2 ( n - 2 - alo - A)la._~, 112] e -2Atp-°o) d o
o oo
-< ~lfo ( l a . - 1.ol If~.ol + 21a.-1,11 [f~,l ~J'-,~-2x(p-po> dp.
(3.48)
Define Ao = n - 3 - 2 m a x alo(p ) = n - 3 - 2 max(0, alo(Po) ).
(3.49)
P>Po
Clearly, alo(Po) < 1 / 2 implies Ao > 0 for all n >_ 4. T h e n for 0 < A < Ao estimate (3.48) implies ¼ [ l a " - l ' ° ( P ° ) 1 2 + 4 1 a " - l ' l ( P ° ) 1 2 ] -< 8 ( A - Ao)
(If~,°12 + Lf",'12)e-E"(°-°°)dP
1
< 16A(A-Ao)
sup [ l f . , o ( p ) 1 2 + [ f . , , ( p ) l / ] .
(3.50)
P~Po
With A = Ao/2 we obtain the optimal estimate
la._l,o(Po ) [2 + 4la._l,l(po)12 < 1 sup (If.,o(P)12 + Lf.,l(p)12) •
(3.51)
"*0 P~Po
N o t e that for odd n it is almost as strong as (3.38) since the minimum of n - 2 + alo(p) in --~o is n - 3 . Since alk(p) , Ikl < 1 tend to 0 as p --* 0% it follows from equations b.k = 0 and estimates (3.44), (3.51) that the coefficients a,k(p) and their derivatives of any order t e n d to 0 as p -* oo. Thus, for fixed r 0 > 0, the p a r a m e t e r a ~ [0, 2"n'] and a point lan(Po)l, alo(po) in -~o where Po = log r o define uniquely a formal series oo
Yas = ~-~ r-nyn(r)
(3.52)
rt=l
as in (3.1). W e will now prove that Ya~ is indeed an asymptotic series to some solution y ( r ) of eq. (1.5) which tends to 0 as r ~ oo. For that sake consider a general system of O D E s ~' = A ( r -~) ~ + f ( ~ , r - l ) ,
(3.53)
518
D. Michelson/ Elementaryparticlesas solutions of the Sivashinskyequation
where A(r -1) = A o + r-lA1 + ~ ( r - 2 ) , A 0 and A 1 are constant matrices, the eigenvalues Ai of A 0 have non-negative real part and are distinct and f ( ~ , r -1) = C([~I 2) depends smoothly on ~ and r -l. Let Ya = ~ ' ( r - l ) b e an approximate solution of (3.53), i.e. Y: - A ( r - l )
Ya - f ( Y , , r - ' ) = g ( r -1) = G ( r - " )
(3.54)
with n to be specified later. In terms of the variable V=(~-Pa)r
m,
m
(3.55)
system (3.53) becomes V' = B ( r - ' ) V + r-'n~(IVl 2) + ~(r--n+m),
(3.56)
where n ( r -1) = h ( r -1) + m r - l I + d~f(ya, r - l ) .
(3.57)
The eigenvalues of the matrix B(r-1) are
Ai+t-ti r-1 + m r - l + ~(lYa[) + C ( r - ~ ) •
(3.58)
Since ~'(lya I) = i f ( r - l ) , for sufficiently large m these eigenvalues have positive real parts as r ~ ~. Note that the derivative ~g, as in the case of eq. (1.5), need not belong to the space L v However, by a similarity transformation T = T O + Tlr -1 one can bring B to the form B = A + ~(Ya) at- i f ( r - 2 ) ,
A =diag{Ai+txi r-1 + m r -1}
(3.59)
so that for sufficiently large m and large r Re B > 0.
(3.60)
System (3.56) could be solved by virtue of the following simple Lemma 3.2. Consider a linear system x' = A ( r ) x + f ( r ) ,
f ~ L I ( 0 , ~ ),
fcontinuous
(3.61)
where A ( r ) is a continuous matrix function such that ReA(r)>0
forr>0.
(3.62)
Then system (3.61) has a unique solution x(r) which tends to 0 as r ~ oo and the following estimate holds Ix(r)[-<2f [f(p)ldp Jr
forall r > 0 .
(3.63)
D. Michelson ~Elementary particles as solutions of the Sivashinsky equation
519
Proof If x(r) is a solution of (3.61) which tends to 0 as r ~ ~, then multiplying (3.61) by x and integrating from r to ~ leads to oo
Ix(r) 12 _
oo
2f tf(p) l Ix(p)Idp < 2sup Ix(p)If If(P) ldp, "r
(3.64)
~r
p~r
which implies estimate (3.63) and proves the uniqueness of such a solution.In order to prove existence, consider solutions x~(r) to the initialvalue problems in backward direction,
xA=A(r) x,+f(r),
O
x,(n)=0.
(3.65)
Multiplying by x n and integrating from r to n leads again to the estimate Ix,(r)l<2 f nIf(p)ldp<2 f~l/(P)l
dp, r
(3.66)
For bounded f, Ix~(r)l are bounded too and hence there exists a subsequence x n which converges on each compact interval to a solution x of (3.61). By (3.66) this solution will satisfy estimate (3.63) and hence will tend to 0 as p --* oo. Now we solve (3.56) by iteration,
r:/+,--B(r-') -v,+I +f,., f~ --r-m:(15,1 z) + ~'(r-"+').
(3.67)
If n > m + 1 > 2 and B satisfies (3.60), the map T: vi ~ Vi÷l is a contraction in a ball I~1~ < C, where the norm is taken in the interval [r0, ~) with sufficiently large r 0. The fixed point of T is the solution of (3.56). Clearly, if B and f depend smoothly on some parameter p, e.g. through the function y~, so does the fixed point of T. The above general framework dearly applies to eq. (1.5). The vector ~ is identified with the triple (y, y', y"). The eigenvalues of the corresponding matrix A(r-1) are +i-r -l+~(r-2),
-2r -1+~(r-2).
(3.68)
Since the coefficients a l k ( P ) , Ik[ < 1 tend to 0 as r ~ , in order to make Re B positive, one needs m > 2 and hence n > 3. This means that in order to estimate y -Yas one has to expand y~ up to order 3 and retain in a3k the coefficients with Ikl = 2,3. We summarize the above arguments in
Theorem 3.1. For each r 1 > 0 and each triple/5 = (la111, a10, a) ~ ~ 0 × $1, there exists a unique solution y(r,/5) of (1.5) which is approximated by the asymptotic expansion Yas in (3.52) determined by the condition la11(Pl)1--la111, al0(Pl)=al0, a r g ( a l l ( P l ) ) = a where p l = l o g r v This solution depends smoothly on/5. With y(r; ,5) defined above let
u(r;/5) --ry(r;/5). For lall I > 0 we consider the map
~- ( ~ 1 , ~ 2 , ~ 3 ) : /5 "~ (½[lu'(rl;/5)l 2+
lu"(rl;/5)12]l/2, u(rl;P)
-tan-l[u'(rfi/5)/u"(rl;/5)]
- rl).
+u"(rl;P),
(3.69)
Clearly ~ = id + d~'(rll). Our next aim is to estimate the image ~ ( . ~ x S 1) for some r I (in our case
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
520
0.85~
0.00
~3
"1J IF -0.50
-t.oo
I~"''
'
0.00
I
. . . .
0.05
I
. . . .
0.t0
I
0. t 5
. . . .
I
0.20
'
'
'
I
0.25
Xt
Fig. 3. The rectangles in x~, x 2 coordinates correspond to the interval values of y'(O) as in table 1. As y'(0) increases, the size of the rectangles grows until the interval of y'(0) is halved. The domain Q1Q2Q3Q4Q5 is displayed in the la111, alo coordinates.
r I = 80), where domain _~ is shown in fig. 3. The purpose of this estimate is as follows. As mentioned in section 1, we are shooting the trajectories of (1.5) in interval arithmetic for interval values of y'(0). As a result of these computations we obtain a box of R 3 which contains the vectors ~(rl) -- (u(rl), u'(rl), u"(rl)). Then we compute the vectors X l = ~1-"U ,2 +U"2) 1/2,
X--~- ( X 1 , X 2 ) ,
X 2 = U q-U"
(3.70)
and obtain a rectangle in the plane xl, x 2 to which the corresponding $(r t) belong. These rectangles for different initial intervals of y'(0) are plotted on fig. 3 and printed in table 1. Suppose that for some domain 2 _ c - ~ 0 we are able to show that ~0(.~ × S')___-~ × S', where ~ includes all these rectangles. This would imply that the solutions y(r) which started at r = 0 with initial conditions y ( 0 ) = 0, y'(0) ~ [ - 0.34, 0.33], y"(0) = 0 exist for all r, tend to 0 as r ~ oo and are approximated by the asymptotics in (3.1). Our first step is to estimate the difference y(rfi i f ) - :9as(rl; ~), where N
yas(r;/5) = Y'~ r-"y,(r),
Yn(r)
a s i n (3.1)
(3.71)
n=l
is a finite asymptotic sum Corresponding to some value of/~ ~ 2 0 × S', while y will always stand for the vector (y, y', y"). Represent u
Y = Yas + rN--------S•
(3.72)
Table 1 The bounds on the displacement of the parameters la~ll and al0 as they are mapped into the plane x~, x:. Computed by the program ASYMP.
R-80. R~80. R-80. R-80.
N-6 -I.<-AI0<-.03 N-6 -I.<-AI0<-.03 N-6 -I.<-AI0<-.03 N-6 -.001<-AI0<-.001
0<-AII<-.230 EI-.2142E-01 E2-.2391E-01 E-.4533E-01 0<-All<-.040 EI-.4479E-02 E2-.4017E-02 E-.8496E-02 0<-All<-.005 EI-.2356E-02 E2-.1246E-02 E-.3602E-02 0<-All<-.020 EI-.4782E-03 E2-.1659E-03 E-.6441E-03
521
D. Michelson / Elementary particles as solutions o f the Sivashinsky equation
Note that (oz+23+r
1)(O+2)rNV__l
=
(
(D+i)(D-i)
rN_ 2
D+
D+i-----7---
rN_l
D-i
r
r
v
(3.73)
and hence v satisfies the equation
(D + i - ~N_2)( D - i N_2)(o N3) r r v + 2YasV = --E(ya~)
v2
r N-1 - - - -
IN-1
(3.74)
,
where f ( Y a s ) is defined as in (3.2). It is convenient to pass from v to the characteristic variables [D + i - ( N - 2 ) r - ' ] [ D - i zl = (-i-r-')(i-r [D + i - ( N Z2=
( N - 2)r-1]v -n)
2)r-ll[D2i(i + r - l )
(N- 3)r-llv '
(3.75)
Z3 = Z 2 "
The above scaling is suggested by the form of the Lagrange interpolating polynomials and should give Z1 + Z2 + Z3 = U + ~ ( r - 2 ) . A direct computation shows ZI+Z2+Z3
=U'}-
~"
N-2 N-3 N-3 ) r 2(--~-_'~+ 1) + 2rZ(ir -1 - 1) + 2 r 2 ( - i r -1 - 1) (3.76)
U.
(1+ r2-~)
In order to obtain the equation for z I one first computes the commutator
[(D+i N-2)(O_ir N-2)r,O N-3]r 2 ( N - 2) =
D2+l
4
_ 2(N-3) r2
(ND+
2)(N-
D - 2(N-3)
N3]r
1)
r2
, D - - -
2 ( N - 2) D + 2 ( N - 2 ) ( N -
2(N-2)(N-3)
r3
r3
r2
1)
r3
(3.77) Eq. (3.74) could then be rewritten as D
N -
r
3
_ r 2 + 1. )(1 + r - 2 ) z l + 2Yasr-T-~+ 2 (Zl + z 2 +z3) v2
---- - - - E ( Y a s ) r N - 1 - - - -
r N-I
+
2v, r2
2 ( N - 1) r3
v
(3.78)
522
D. Michelson /Elementary particles as solutions o f the Sivashinsky equation
or finally D
N-
3
) ZI 4" 2Yas(Zl
4"Z2 4"Z3)
(3.79)
= r-zfl
where
fl
= (1 +
[ r-2)-l[-~(Yas)
2(N-
v2 r N + l -- -r N-_ 3 -4- 2U'
r
1)
1 V + 2r-Izl
J
4r 2 + r-5--~+ 2 Yas( 21 + Z2 + 23).
(3.80)
Since the operators D + i - (N - 2)r- 1 commute between themselves we obtain from (3.74) the equation for z 2 D
-
i
N -r2 )
(3.81)
22--Yas(Zl 4-22 4-23) --r-2 fz,
where
1
f2 = 2 ( 1 - ir -1)
_~(Yas) rN+ 1
c2 + u---'-'-~ r
) r 2 ( i r - 2) -2iz 2 +
"z
r-~+-2 Yas( ,
+z2 +23)
•
(3.82)
Now multiply eq. (3.79) by 2'1, eq. (3.80) by 422, add, integrate from r I to oo and take real part. Since z 1 is real and z 3 = 2,2, the cross terms with ya~ cancel and we obtain
½[izl(r,)12 + 4]z2(rl)12] + [~(NJr I
k
3
]z112 + N-
r
2 r
'2z2'2)dr
oo
4- f~(-2Yas)(lZl12-4]Rez212)dr N frrlr-2(lflz,] -I-4]f2z2l)dr.
(3.83)
rl
For positiveness of the sum of the integrals at the left-hand side of (3.83) one needs -(N-2)
<2Uas(r ) < N - 3
f o r r > r 1,
(3.84)
where uas = rYas al0(p) + 2 R e ( a u ( p ) e jr) + & ( r - l ) . In the domain 2 in fig. 3, alo varies between - 1 and 0.03 while lalll varies between 0 and 0.23. Hence, the leading part of 2Uas may vary between -2.92 and 0.98. Thus N should be at least 5. Assume that Ua~ satisfies (3.84). Then =
o¢
½[Izl(r') 12+ 4lz2(rl) 12] -< fr, r-2(Iz'12
+
4122t2)1/2(Iflt2 4. 41f212)1/2 dr
1
_< -- sup (Iz, 12 + 412212) 1/2 sup (If112 + 41f212) 1/2. rl
r>r i
r>_r I
(3.85)
D. Michelson / Elementary particles as solutions o f the Sivashinsky equation
523
Since r 1 in (3.85) may be replaced by any number greater than r 1, we obtain (3.86)
IlzlL - 211flL, where
Ilzll= = sup [ Izl(r ) Iu + 4lZz(r ) 12] ,/2
(3.87)
r>r I
and similarly for f. A stronger estimate could be obtained if we assume that --3<2Uas<2
and
N>6.
(3.88)
Then
½[Iz1(rl)l2+41z2(rl)l2] +frS-N~(Iz,l2+41z212)dr <
frS~-~( Iz'12+4lz212)dr+
4 ( N 1- 5) frSr-3 ([f'12 + 4Lf212) dr
(3.89)
and hence 1
IlzlL ~ zv N - '3' c~r 7r-,- - ~ Ilfllo~.
(3.90)
Although IlflL depends on IlzlL, with estimates (3.86) or (3.90) one can find bounds of IlzlL. Recall that by (3.76) df
Ivl ~ Izl = IZll +21221.
(3.91)
In order to express v' in terms of z we consider the sum N-3+2Re((i+N-2)
Zl-----~ -
+2Re
r
z2
)=v,
+
i+----7--- 2r2(~_--1-_1 ) v
(N-3)(N-2) r3(r+r_2)
+r(r2+l)
v V.
(3.92)
Hence Iv'l-< 21z21 + r - ' ( g - 2 ) ( I z ,
I + 2[z21 + Izlr-2)
21z21 + r - l ( N - 2)Izl(1 + r-2).
(3.93)
For v"we use
zl=v,,+v
2(N-
2)v' + r
(N-2)(N-1)v r2
(3.94)
D. Michelson /Elementary particles as solutions o f the Sivashinsky equation
524
Hence, in view of (3.76), Iv"l <21z21 + 2r-l( N - 2)[v'l + r-2lcl[( N - 2)( N - 1 ) +1].
(3.95)
We will need also the estimate
Iv + v"l ~< Iz, I +
2 ( N - 2) (Nr Iu'l +
2)(Nr2
1)
Ivl.
(3.96)
Now let us estimate Ilfll~ in terms of IIzlPo=.By (3.80) [ft 1= < [ J (
Yas) r N + 11 oo -t" Izl2/r N - 3 + 4 Iz21~ + 2( N - 2) r~ 1]zloo
+ 2(N-
1)ri-11z]~ + 2ri- 11zl [~ + 4r~-lluasl~ Izloo,
(3.97)
where (3.98)
Izl~= Iz~l~+21z21~= sup[Iz~(r) l + 2 1 z 2 ( r ) l ] . r>~rI
Similarly
2]fz[~ < I f ( Y a s ) rN+ll~
+
Iz]2/r N-3 + 21z21~ + 21Uasl~ IZI~
(3.99)
and hence tlfll~-<
g~-l-zP(Yas) rN +l l~ + V~ ]Zl2/ r N - 3 + (201z212 + 81Uasl~ IZ21~ Izl~ + 41Ua~l 2 IZ12) ~/2 + 4r~-1( N - 1 + lUa~l~)Izl~.
(3.100)
Clearly,
21z21~ -< Irzll~,
(3.101)
Izl~ -< v~llzll~.
Hence
Ilflloo --- 7~- [~zP(Yas) rg+llo~ + 2v~llzll2/rN-3 + C~llzll~
(3.102)
where C1
=
(5 + 4V~IUa~I~ + 8lUaslZ) 1/2 + 4 v ~ r [ l ( g -
1 + lUasl~).
(3.103)
For example, if Ua~ is bounded as in (3.88), r 1 = 80 and N_< 6 C 1 < 6.1.
(3.104)
D. Michelson /Elementary particles as solutions o f the Sivashinsky equation
525
With estimates (3.86) and (3.90) rewritten in a uniform way C2
IIzL= ~ -~(llflL,
C2 = 2
or
1 C2 = 2 Nil--S_ 5
(3.105)
we obtain the final bound on Ilzll®: 2C2v/2 "
IlzlL®-<
rI
l-z:(Yas) r N + I ~//C 3,
(3.106)
where 16CzZl.z:(Yas) rN+ll~lrlU-') '/2,
C 3 = C 4 + (C42 -
C 4 =
1-
ClC2//rl
.
(3.107)
Now, for the map ~0 in (3.69) we are able to estimate the difference A~p = q~ - id = ( A ~ , Aq~2, A~3 ) .
(3.108)
The vector R(r; 13) = (Uas , Uas , U;ls) + r2-N(V, v' -- ( N - 2)v/r, v" - 2 ( N - 2)v'/r + ( N - 1 ) ( N - 2) v/r2), (3.109) where N
Uas = E r-n+l ~, ank eikr,
(3.110)
Ikl~n
n=l N
u's=
E r--+' E { [ i k - - ( n - - 1 ) / r l a n k + r - l D p a n k } e
ik',
(3.111)
Ikl<_n
n=l N
u~s= E r -"+1 Y'. { - - k 2 a n k + 2 i k [ ( - n + 1)a~k + Dpa,k]/r n=l
Ikl <_n
+[n(n-
1)a~k-- 2 ( n - - 1)Doa~k-- D~a~k]/r }e
2
2
ikr
.
(3.112)
When Uas, u',~ and u"~ are estimated, all terms at the right-hand side are replaced by their absolute values. In order to obtain Aqh, Aq~2 one has to subtract from Ua~, u~ and u~s the terms of order 0, i.e. AUas= U a s - ( a l o + 2 R e a n e i r ) , Au~
= U a"s
- 2 Re
a l l e ~r.
Au'~=Ua~+2Ima n
e ir '
(3.113)
526
D. Michelson /Elementary panicles as solutions of the Sivashinsky equation
Then 21A~,I ~ (IAu;~l 2 + [AUasl2) 1/2 + (U '2 q+(g-
u,,2)l/2//r N-2
2){[vl 2 + [2lt"[ + ( g - 1 ) I u [ / r ] 2 } l / 2 / r g - 1 = 2 E l ,
(3.114)
while [A~21 < ImUas+ m U a s l - 4 - I c + v " l / r +(g-2)[Ivl
N 2
+ 21v'l + ( g -
1 ) l v l / r ] / r N-L.
(3.115)
Since we fix r = r 1, the norms of Ivl, IL"I and Iv"l are estimated as in (3.91), (3.93) and (3.95) with and Iz[ bounded as in (3.101). In Auas + AUas there are some cancellations,
IAUa, +
N AUasl < Y', r - " + ' E / l k 2 - 11(1 - a n ~ ) a , e
+ 21kl [(n
- 1)lankl
Iz21
+ [Doa,kl]/r
n=l
+ [n(n-
1)la,k [ + 2 ( n -
1)lDpa,k I +
ID2ank[]/r 2}
(3.116)
and also It' + v"l is estimated as in (3.96)with Iz~l -< I[zl[~. As for the difference Aq~3 notice that - t a n - l ( u ' / u " ) is the argument of the complex vector u " - iu' while a + r is the argument of its leading part 21a1~ [ e i(r+a). In order to prove that for fixed lal~l > 0 and al0, ~03 maps S I onto S 1, one has to verify
lu"-
i u ' - 21all I ei(r+")[ < 21all I
(3.117)
for all a ~ S 1. However, the left-hand side of (3.117) is bounded by the right-hand side of (3.114). Thus, one has to check only E1 < lalll.
(3.118)
Now, all the machinery for the computation of A~01,A~2 has been prepared. The computations, however, are very cumbersome to be carried out by hand. Indeed, the order N is at least 5 or 6 so that the number of coefficients a,k and their derivatives needed is about 32N3 = 144. Instead we wrote a computer program to do all the estimates which are described in the next section.
4. The numerical estimate of the asymptotic expansion Let 2 be a domain in the plain lall [, alo bounded by the contour PIP2P3P4P5 or contour P1P2P3P4'P5' as shown in fig. 4. Here P1P2 is a horizontal line a~o = al0max with aloma x > 0, P2P3 is a trajectory of the flow (1.10), P3P4 or P3P~ is a vertical line lall] =allmax, P4P5 a horizontal line alo=alomin, where - 1 < alomin < 0, while P4'Ps' is a segment of the stable manifold Mst of the critical point Ps' = ( 0 , - 1 ) . Clearly the flow in (1.10) enters 2 at the boundaries P1P2,P3P4 and at a part of P4Ps. Observe that P3P4 is tangent to the arcP2P3 at P3. In order for 2 to be forward invariant, the flow should enter 2 at the whole segment P4Ps. For that sake one has to request that on PaPs, a~o > 0, namely al0 min -F a20 mi. + 2a 2, max
< O.
(4.1)
D. Michelson ~Elementary particles as solutions o f the Sivashinsky equation
527
OlotPl
P2 *]P3 t OtIp5 D it°lP4lmox I%111
OlOma:
°'°n"
i/Ms,
_[,L p;
Fig. 4.
In case of the boundary PIP~, the domain -~ is always forward invariant. The lalll coordinate of P2 could be estimated as follows. By (1.10) for al0 > 0, dla111 lalllalo > - lall I . dalo = - 4 1 o+a2o+21a1112
(4.2)
Hence logJalll(P2) - l o g l a ] l l ( P 3 ) > -a]omax,
(4.3)
e--al0maxallmax ~ [alll(P2) ~ allmax.
(4.4)
i.e.,
In the case of the boundary P~P; we set alomi. = - 1. Thus for lall(P])l, a t o ( P l ) ~ ..~ there is a uniform estimate sup 1a11( p l ) j
<
allmax,
P >P]
sup lalo(p) J < max(alomax, - a l 0 m i n ) .
(4.5)
P ~P]
In our computer code ASYMP the variable sup
sup
A(n, k, l) bounds
IDta,,k(p)l
(4.6)
(lajl(ol) I, alo(Pl)) ~ -~ P >Ol
Hence we set A ( 1 , 0 , 0 ) = max(a]oma x , - alomin),
Z ( 1 , 1,0) = a]lmax.
(4.7)
For the derivatives lahJ and Ja;ol by (1.10) we have A(1, 1, 1) =A(1,O,O) A(1, 1,0)
(4.8)
A(1,0, 1) = max (a]o + a2o + 2[a]11 z) = max(Yma ~ + 2a11 ma x , JYminl),
(4.9)
and
D. Michelson / Elementary particles as solutions of the Sicashinsky equation
528
where 2 Ymax = a l 0 rnax q" a l 0 max, Ymin = a l 0 m i n "k- a 2 0 m i n
if al0mi n > - 0 . 5 ,
Ymin =
otherwise
-- 0 . 2 5 .
(4.10) In order to decrease the values of Imq91,21 to an admissible level we were forced to use refined estimates also for the second derivatives a10, ae~, It
a7, = -a,~(a,o + 2[aHI2),
vt
a~'o=a,, + 3a2o + 2a3o + 21a1112(1+a~o),
(4.11)
which leads to the definition A ( 1 , 1 , 2 ) = all A ( 1 , 0 , 2) =
max
max(al0 max +
max
2a~lmax, -
(al0 + 3a20 +
(4.12)
al0 min),
2a30) + 2al2lmax(1 + al0max).
(4.13)
t/10 max ~
For la21(p)l we use estimate (3.44), which implies (4.14)
~ A(2, 1,0) = all max ( 1 -1- 2al0 max + T19a l2l maxl"
All the rest of the coefficients a~k and their derivatives are estimated with the help of equations b~k(p) = 0 (see 3.4) and estimates (3.38), (3.51). Formula (3.4) is used to define the function COEF(n, k, l) which majorates ID~b~k I. Namely,
COEF( n, k, l ) = k l - k 2 +
l l A ( n , k , l ) + 13k 2 - l [ A ( n - 1 , k , l + 1)
+ 1(3n - 7)k 2 - n + 31A(n - 1, k, l) + k[3A(n - 2, k, l + 2) + 16n - 17lA(n - 2, k , l + 1) + 13n 2 - 17n + 22[A(n - 2, k , / ) ]
+ A ( n - 3 , k , l + 3) + l - 3n + l O I A ( n - 3, k , l + 2) + l ( n - 2 ) ( n - 3) + ( n - 2 ) ( n - 5) + ( n + [(n - 2 ) ( n -
+
E
E
3)(n - 5)lA(n-
3, k , l + 1)
3 ) ( n - 5 ) l A ( n - 3, k , l )
E
( tll ) A ( n " k l , l l ) A ( n z , k 2 ' /2)
(4.15)
nj +n2=n kl +k2=k 11+12=l
Suppose that all coefficients A(nl,kl, l 1) with n l < n - 1 , Ikll < n 1 and l l < n + l - - n I were already computed while A(n, k, l) with k 4~ 0, 1 was set to be 0. Then we compute COEF(n, k, l) by formula (4.15) and define
A ( n , k, l) = C O E F ( n , k, l ) / ( k [ -
k 2 + 11).
(4.16)
or 1, the functions fnk in (3.10) could be also estimated by COEF(n, k, 0). Indeed, if A ( n - 1 , k,l) for Ikl < 1 and l < l are 0, then COEF(n, k, 0) majorates the functions fnk. Hence If k = 0
D. Michelson / Elementary panicles as solutions of the Sivashinsky equation
529
A(n - 1,k, 0) for even n and k = 0, 1 could be defined as 2 1/2
A ( n - 1,0,0) -- [ C O E F ( n , 0 , 0 ) 2 + C O E F ( n , 1,0) ]
/ ( n - 3 - 2al0m~x) ,
A ( n - 1, 1,0) = A ( n - 1 , 0 , 0 ) / 2 .
(4.17) (4.18)
For odd n,
A(n-
1,0,0) = 0,
A(n-
1,1,1) = C O E F ( n , l , 0 ) / ( n -
2 +al0min).
(4.19)
If the terms A(n - 1, k, l~) for k = 0, 1, l 1 _
A(n-
1 , k , l + 1) = C O E F ( n , k , l ) / 1 3 k 2 - 11.
(4.20)
For negative k we use identities
A(n,-k,l)=A(n,k,l),
COEF(n,-k,l)=COEF(n,k,l).
(4.21)
Since for odd n, b~o - 0 we define C O E F ( n, 0, /) -- 0
f o r o d d n.
(4.22)
Recall that we wish to estimate Yas as in (3.71) and [.E(yas)rN+tl~. In y,~ we assume that aN, k for Ik[ < 1 are identically zero. The function -f(Ya~) is 2N
"~(Yas)
=
E n=N+l
r-" E b.k eikr,
(4.23)
Ikl
where b~k are defined by (3.4) but with a~k =- 0 for n > N and aN, k = 0 for [k I < 1. Thus for n > N + 1, b~k are estimated by
B(n, k) = COEF(n, k,0)
(4.24)
and I.t(Yas)[® is bounded by 2N
LYAS=
Y~. r? ~ E B ( n , k ) . n=N+l
(4.25)
Ikl
Then [ - ~ ( Y a s ) r N + l [ ~ --< r ~ +t LYAS.
(4.26)
The order of computations is as follows. Initially the array A ( n , k , l ) for - 2 < n < 2 N , Ik[ < 2 N , 0
530
D. Michelson / Elementaryparticles as solutions of the Sicashinsky equation
loop (we mimic some F o r t r a n n o t a t i o n s ) do 20 n = 2 ,
N
if n < 3 go to 10 if n even c o m p u t e A ( n - 1,0,0), A ( n - 1, 1 , 0 ) b y (4.17), (4.18) if n odd c o m p u t e A ( n - 1, 1,0) by (4.19)
A ( n - 1, - 1,0) = A ( n - 1, 1,0) 10do20
l=0, N+3-n
do 20 k = 0, n if k > 2 A ( n , k , l )
= COEF(n,k,l)/k/(k
2 - 1)
A ( n , - k , l ) = A ( n , k, l) else if n = 2 a n d l < 1 go to 20; (A(1, k, 1) were already c o m p u t e d )
A ( n - 1, k, 1 + 1) = C O E F ( N , K, L ) / 1 3 k 2 - 11 if n = 2 a n d k = 0, A ( 1 , 0 , l + 1) = A ( 1 , 0 , l + 1) - 2 ( A ( 1 , 0 , 0 ) - al0max) A(1,0, l)
A ( n - 1, - k , l
+ 1) = A ( n - 1, k , l , + 1)
20 e n d l o o p do30 n=N
+ I,2N
do 30 k = 0, n 30 B ( n , k ) = C O E F ( n , k , 0 ) LYAS = 0 UAS = 0 do 50 n = N + SUM = 0
1,2N
do 40 k = 1, n 40 S U M = S U M + B ( n , k ) S U M = 2 S U M + B ( n , O) 50 L Y A S = L Y A S + S U M / r ~ ' do 70 n = 1, N SUM = 0 do 60 k = 1, n 60 S U M = S U M + A ( n , k,O) SUM = 2SUM +A(n,O) 70 U A S = U A S + S U M / r ~ ' - 1 T h e c o m p l e t e p r o g r a m A S Y M P a p p e a r s in the appendix. H e r e we show the logic of the loop. Notice the
D. Michelson /Elementary particles as solutions of the Sivashinsky equation correction of A(1, 0,1 + 1) in the loop "do Ii')1+1-L,,p Ul0 =
IDtoal0 + 2al0D~al0 t +
20
531
k = 0, n". Indeed
+2 ~t (111)D~,lallID~-t~lal~ I .
l~- 1 ( 11 l ]x..,p]r~t~',,10,..or~t-t",,10 11= 1
(4.27)
li = 0
When we use formula (4.20), the first two terms on the right-hand side of (4.27) are replaced by A(1, 0, l) and 2A(1, 0, 0)A(1, 0, l) but actually their sum is ID~alo(1 + 2alo) I < A ( 1 , 0 , I ) (1 + 2al0max).
(4.28)
Thus we correct A(1, 0,1) by the difference Z ( 1 , 0 , l) [1 + 2 A ( 1 , 0 , 0 ) ] - h ( 1 , 0 , 1 ) ( 1 + 2al0max) = 2 A ( 1 , 0 , I) [ h ( 1 , 0 , 0 ) - al0max].
(4.29)
The domain ~ in fig. 3 bounded by the c o n t o u r Q I Q 2 Q 3 Q 4 Q 5 Q I is covered by a union of domains -~i of the type P1P2P3P~P ~ shown in fig. 4. The point QI = (0.004, 0.03), Q4 has the x 1 coordinate 0.23, QIQ5 is a vertical line and the boundary Q1Q2Q3Q4Q5 has the same meaning as P~P2P3P~P~ in fig. 4. In table 1 the results of computation by the program ASYMP are exhibited. For example, in the first line of the table R = 8 0 means r 1 =80, N = 6 is the order of the expansion )'as" The bounds al0min = --1, a10 max = 0.03, all max = 0.23 characterize uniquely the corresponding domain -~1 of the type PIP2 P3 P4'P~. By E1 and E2 we denote correspondingly the maximal values of IA~,11 and IA~21 defined in (3.114) and (3.115), while E bounds (A~Ol2 + A~o22)1/2. Clearly the domain .~ lies in -~1 and hence IA~0tl and IAq,21 there are bounded by E1 --- 0.0215 and E2 -- 0.0240. Let ff -- (qh, ~02) denote the restriction of the map in (3.69) to its first two components. It is clear from fig. 3 and the above bounds El, E2 that the image of the boundary Q1Q2Q3Q4 under the map ff for all values of the parameter a lies away from the rectangles in fig. 3. Along the boundary Q4Q5 the same is true for [al~l >_ 0.04. Also for (laHI, alo) ~ -~l and laHI > 0.022, IA~otl < E1 < laH[, i.e. (3.118) holds. The domain -~2 with a t l mxa = 0.04 takes care of the remaining portion of the boundary Q4Qs. We should be concerned only with the lowest rectangle whose lower edge has coordinates xl ~ [0.007,0.0114] and x2 = -0.986 (see the last line of table 2). According to (3.22), the ato coordinate of P4'P~ for la111 < 0.04 satisfies a10 < - 1 + ~ • 0.04 2 + 3" 0.043 < -0.9989. With the bound E2 in ~ 2 clearly -0.9989 + 0.0041 < -0.986, which means that ff(P4Ps) stays away from all rectangles. Also for 0.005 < lanl < 0.04 the condition in (3.118) holds. In order to treat the boundary QIQ5 we need the additional domain .~r 3 with a l l mxa = 0.005, where IAqhl < 0.0024, IA~021< 0.0013. The laHI coordinate of QIQ5 is 0.004 and hence q~l(Q1Qs) is away from x I = 0.0064, while the left edge of the lowest rectangle has x 1 coordinate 0.007006. Also, for (laHI, a~0) ~ -~3 and lat~l -- 0.004 condition (3.118) holds. There is only a problem with the rectangle which originates at x~ = 0 and corresponds to ly'(0)l < 0.005. But otherwise tables 1 and 2 show that for -0.34015625 < y'(0) < 0.33 and ly'(0)l > 0.005 the solution of (1.5) falls in the domain of existence of the asymptotic solution and hence tends to 0 as r ---}oo. In order to treat the range of initial conditions ly'(0)l < 0.005 one has to normalize the variable y(r) by the factor e = y'(0). Namely, we change variables
yAr)
y'(O)
=
(4.30)
The new variable y~ satisfies the equation (4.31)
532
D. Michelson /Elementary particles as solutions of the Sit~ashinsky equation Table 2 The coordinates of the rectangles in the xl, x 2 plane for the corresponding interval values of y'(0). Computed by the program RUN. Y'(0)
U+U''
. 5 * ( U ' * * 2 + U ' ' * * 2 ) * * I / 2
-.50000000E-02 .50000000E-02 -.II01E+00 .II03E+00 .I139E+01 .1911E+01 T h e a b o v e v a l u e s of U,U' and U'' are d i v i d e d b y t h e c o r r e s p o n d i n g Y'(0) -.50000000E-02 .50000000E-02 .15000000E-01 .25000000E-01 .35000000E-01 .45000000E-01 .55000000E-01 .65000000E-01 .75000000E-01 .85000000E-01 .95000000E-01 .10500000E+00 .II500000E+00 .12500000E+00 .13500000E+00 .14500000E+00 .15500000E+00 .16500000E+00 .17500000E+00 .18500000E+00 .19500000E+00 .20500000E+00 .21500000E+00 .22500000E+00 .23500000E+00 .24500000E+00 .25500000E+00 .26000000E+00 .26500000E+00 .27000000E+00 .26500000E+00 .27000000E+00 .27500000E+00 .28000000E+00 .28500000E+00 .29000000E+00 .29500000E+00 .30000000E+00 .30250000E+00 .30500000E+00 .30750000E+00 .31000000E+00 .31250000E+00 .31500000E+00 .31750000E+00 .32000000E+00 .32125000E+00 .32250000E+00 .32375000E+00 .32500000E+00 .32625000E+00 .32750000E+00 .32875000E+00 .32937500E+00
50000000E-02 15000000E-01 25000000E-01 35000000E-01 45000000E-01 55000000E-01 65000000E-01 75000000E-01 85000000E-01 95000000E-01 10500000E+00 II500000E+00 12500000E+00 13500000E+00 14500000E+00 15500000E+00 16500000E+00 17500000E+00 18500000E+00 19500000E+00 20500000E+00 21500000E+00 22500000E+00 23500000E+00 24500000E+00 25500000E+00 26000000E+00 26500000E+00 27000000E+00 27500000E+00 27000000E+00 27500000E+00 28000000E+00 28500000E+00 29000000E+00 29500000E+00 30000000E+00 30250000E+00 30500000E+00 30750000E+00 31000000E÷00 31250000E+00 31500000E+00 31750000E+00 32000000E+00 32125000E+00 32250000E+00 .32375000E+00 .32500000E+00 .32625000E+00 .32750000E+00 .32875000E+00 .32937500E+00 .33000000E+00
-.1405E-03 -.I074E-02 -.2974E-02 -.5849E-02 -.9693E-02 -.1448E-01 -.2015E-01 -.2664E-01 -.3387E-01 -.4174E-01 -.5014E-01 -.5898E-01 -.6813E-01 -.7750E-01 -.8701E-01 -.9657E-01 -.I062E+00 -.I158E+00 -.1254E+00 -.1352E+00 -.1452E+00 -.1557E+00 -.1673E+00 -.1809E+00 -.1985E+00 -.2252E+00 -.1941E+00 -.2022E+00 -.2116E+00 -.2227E+00 -.2116E+00 -.2227E+00 -.2360E+00 -.2526E+00 -.2735E+00 -.3014E+00 -.3414E+00 -.3334E+00 -.3533E+00 -.3763E+00 -.4033E+00 -.4354E+00 -.4742E+00 -.5222E+00 -.5839E+00 -.5926E+00 -.6277E+00 -.6674E+00 -.7129E+00 -.7653E+00 -.8265E+00 -.8990E+00 -.9264E+00 -.9693E+00
1405E-03 1583E-03 7457E-03 2608E-02 5439E-02 9225E-02 1393E-01 1950E-01 2585E-01 3289E-01 4051E-01 4859E-01 -.5700E-01 -.6560E-01 -.7428E-01 -.8290E-01 -.9134E-01 -.9948E-01 -.I072E+00 -.I145E+00 -.1212E+00 -.1272E+00 -.1325E+00 -.1364E+00 -.1378E+00 -.1327E+00 -.1829E+00 -.1895E+00 -.1968E+00 -.2050E+00 -.1968E+00 -.2050E+00 -.2142E+00 -.2246E+00 -.2360E+00 -.2480E+00 -.2576E+00 -.3115E+00 -.3275E+00 -.3457E+00 -.3663E+00 -.3897E+00 -.4163E+00 -.4462E+00 -.4788E+00 -.5547E+Q0 5843E+00 - 6173E+00 6542E+00 6956E+00 - 7424E+00 - 7952E+00 - 8759E+00 - 9131E+00 -
.0000E+00 .7412E-02 .2268E-01 .3801E-01 .5328E-01 .6832E-01 .8298E-01 .9710E-01 .II05E+00 .1231E+00 .1348E+00 .1453E+00 .1548E+00 .1629E+00 .1698E+00 .1754E+00 .1796E+00 .1823E+00 .1835E+00 .1827E+00 .1796E+00 .1746E+00 .1677E+00 .1583E+00 .1454E+00 .1265E+00 .1506E+00 .1449E+00 .1387E+00 .1319E+00 .1387E+00 .1319E+00 .1243E+00 .I157E+00 .I059E+00 .9381E-01 .7751E-01 .9068E-01 .8560E-01 .8022E-01 .7447E-01 .6819E-01 .6115E-01 .5288E-01 .4249E-01 .5262E-01 .4941E-01 .4605E-01 .4247E-01 .3859E-01 .3430E-01 .2938E-01 .3385E-01 .3206E-01
.7612E-02 .2287E-01 .3824E-01 .5361E-01 .6877E-01 .8357E-01 .9786E-01 .III5E+00 .1243E+00 1363E+00 1472E+00 1571E+00 1659E+00 1735E+00 1799E+00 1851E+00 1892E+00 1921E+00 1939E+00 1951E+00 1962E+00 1964E+00 1960E+00 1955E+00 1959E+00 2006E+00 1596E+00 1550E+00 1500E+00 1448E+00 1500E+00 1448E+00 1395E+00 1343E+00 1295E+00 1262E+00 1268E+00 9891E-01 9487E-01 9086E-01 8699E-01 8340E-01 8038E-01 7853E-01 7918E-01 5952E-01 5728E-01 5518E-01 5329E-01 5173E-01 5071E-01 5056E-01 4156E-01 4081E-01
T h e f o l l o w i n g v a l u e s of Y'(0) d o not lie in the d o m a i n of c o n v e r g e n c e .33000000E+00 .33062500E+00 -.1016E+01 -.9536E+00 .3019E-01 .4021E-01 .33062500E+00 .33125000E+00 -.I068E+01 -.9977E+00 .2824E-01 .3981E-01
D. Michelson / Elementary particles as solutions of the Sivashinsky equation Table 2
15000000E-01 25000000E-01 35000000E-01 45000000E-01 55000000E-01 65000000E-01 75000000E-01 85000000E-01 95000000E-01 10500000E+00 II500000E+00 12500000E+00 13500000E+00 14500000E+00 15500000E+00 16500000E+00 17500000E+00 18500000E+00 19500000E+00 20500000E+00 21500000E+00 22500000E+00 23500000E+00 24500000E+00 25500000E+00 26500000E+00 27500000E+00 28000000E+00 28500000E+00 29000000E+00 29500000E+00 30000000E+00 30500000E+00 31000000E+00 31250000E+00 31500000E+00 31750000E+00 32000000E+00 32250000E+00 32500000E+00 32625000E+00 32750000E+00 32875000E+00 33000000E+00 33125000E+00 33250000E+00 33375000E+00 33500000E+00 33562500E+00 .33625000E+00 .33687500E+00 .33750000E+00 .33812500E+00 .33843750E+00 .33875000E+00 .33906250E+00 .33937500E+00 .33968750E+00 .34000000E+00 .34015625E+00
Below are the -.50000000E-02 -.15000000E-01 -.25000000E-01 -.35000000E-01 -.45000000E-01 -.55000000E-01 -.65000000E-01 -.75000000E-01 -.85000000E-01 -.95000000E-01 -.I0500000E+00 -.II500000E+00 -.12500000E+00 -.13500000E+00 -.14500000E+00 -.15500000E+00 -.16500000E+00 -.17500000E+00 -.18500000E+00 -.19500000E+00 -.20500000E+00 -.21500000E+00 -.22500000E+00 -.23500000E+00 -.24500000E+00 -.25500000E+00 -.26500000E+00 -.27500000E+00 -.28000000E+00 -.28500000E+00 -.29000000E+00 -.29500000E+00 -.30000000E+00 -.30500000E+00 -.31000000E+00 -.31250000E+00 -.31500000E+00 -.31750000E+00 -.32000000E+00 -.32250000E+00 -.32500000E+00 -.32625000E+00 -.32750000E+00 -.32875000E+00 -.33000000E+00 -.33125000E+00 -.33250000E+00 -.33375000E+00 -.33500000E+00 -.33562500E+00 -.33625000E+00 -.33687500E+00 -.33750000E+00 -.33812500E+00 -.33843750E+00 -.33875000E+00 -.33906250E+00 -.33937500E+00 -.33968750E+00 -.34000000E+00
533
Continued.
negative values of Y' -.I017E-02 .1327E-03 (0!7373E-02 -.2700E-02 -.7431E-03 .2195E-01 .3601E-01 -.5104E-02 -.2419E-02 -.8132E-02 -.4809E-02 .4942E-01 -.I168E-01 -.7816E-02 .6210E-01 -.1565E-01 -.I134E-01 .7396E-01 8494E-01 -.1993E-01 -.1526E-01 -.2442E-01 -.1948E-01 9499E-01 -.2901E-01 -.2390E-01 1041E+00 -.3362E-01 -.2838E-01 i122E+00 I193E+00 -.3816E-01 -.3285E-01 -.4254E-01 -.3722E-01 1254E+00 -.4670E-01 -.4141E-01 1305E+00 -.5058E-01 -.4536E-01 1346E+00 -.5413E-01 -.4900E-01 1377E+00 -.5734E-01 -.5230E-01 1399E+00 -.6019E-01 -.5521E-01 1404E+00 -.6268E-01 -.5773E-01 1396E+00 -.6485E-01 -.5987E-01 1377E+00 -.6677E-01 -.6166E-01 1349E+00 -.6854E-01 -.6312E-01 1310E+00 -.7036E-01 -.6431E-01 1260E+00 -.7241E-01 -.6538E-01 1200E+00 -.7505E-01 -.6645E-01 I125E+00 I034E+00 -.7888E-01 -.6761E-01 -.8488E-01 -.6885E-01 9236E-01 -.9519E-01 -.6942E-01 7806E-01 -.9137E-01 -.8504E-01 8336E-01 -.9734E-01 -.8939E-01 7752E-01 -.I049E+00 -.9473E-01 .7134E-01 -.I146E+00 -.1013E+00 .6474E-01 -.1272E+00 -.I093E+00 .5755E-01 -.1441E+00 -.1187E+00 .4935E-01 -.1682E+00 -.1291E+00 .3899E-01 -.1726E+00 -.1579E+00 .4458E-01 -.1874E+00 -.1698E+00 .4108E-01 -.2051E+00 -.1837E+00 3743E-01 -.2266E+00 -.1998E+00 3352E-01 -.2531E+00 -.2187E+00 2916E-01 -.2869E+00 -.2407E+00 2396E-01 -.2996E+00 -.2789E+00 2744E-01 -.3208E+00 -.2970E+00 2567E-01 -.3450E+00 -.3175E+00 2383E-01 -.3731E+00 -.3408E+00 2189E-01 -.4059E+00 -.3673E+00 1979E-01 -.4448E+00 -.3979E+00 1745E-01 -.4917E+00 -.4334E+00 1471E-01 -.5495E+00 -.4745E+00 I129E-01 -.5741E+00 -.5383E+00 1511E-01 -.6102E+00 -.5697E+00 .1403E-01 -.6509E+00 -.6047E+00 .1285E-01 -.6972E+00 -.6439E+00 .I155E-01 -.7503E+00 -.6881E+00 .1006E-01 -.7770E+00 -.7438E+00 .I096E-01 -.8084E+00 -.7724E+00 .1030E-OI -.8424E+00 -.8032E+00 .9595E-02 -.8793E+00 -.8364E+00 .8826E-02 -.9195E+00 -.8723E+00 .7984E-02 -.9633E+00 -.9112E+00 .7052E-02 -.9860E+00 -.9556E+00 .6998E-02
.2225E-01 .3637E-01 .4986E-01 .6262E-01 .7457E-01 .8565E-01 .9581E-01 .1050E+00 .I133E+00 .1205E+00 .1268E+00 .1321E+00 .1365E+00 .1399E+00 .1423E+00 .1438E+00 .1448E+00 .1454E+00 .1450E+00 .1435E+00 .1411E+00 .1377E+00 .1333E+00 .1281E+00 .1224E+00 .I165E+00 .1119E+00 .9098E-01 .8593E-01 .8082E-01 .7576E-01 .7100E-01 .6699E-01 .6502E-01 .4986E-01 .4704E-01 .4435E-01 .4192E-01 .3995E-01 .3888E-01 .3050E-01 .2910E-01 .2776E-01 .2653E-01 .2546E-01 .2464E-01 .2425E-01 .2467E-01 .1800E-01 .1740E-01 .1686E-01 .1642E-01 .1612E-01 .1369E-01 .1336E-01 .1304E-01 .1274E-01 .1247E-01 .1224E-01 .I132E-01
The following values of Y'(0) do not lie in the domain of convergence .34031250E+00 -.34015625E+00 -.1011E+01 -.9784E+00 .6499E-02 .III8E-01 .34046875E+00 -.34031250E+00 -.I036E+01 -.1002E+01 .5969E-02 .ll04E-01
D. Michelson /Elementary particles as solutions of the Siuashinsky equation
534
with the initial condition y~(0) =y~'(0) = 0,
y~(0) = 1.
(4.32)
We solve this equation numerically in interval arithmetic in the same way as eq. (1.5) but with the interval value of e = [ - 0.005, 0.005]. The resulting interval values of x 2 = u~ + u~n and x~ = ~1,"u ~,2 + U~,2)1/2 at r = r~ = 80, where u~ = rye a p p e a r in the first line of table 2. Clearly
Ix2/xll ~ 0.1103/1.13 < 0.1.
(4.33)
The same ratio is maintained for x~ and x 2, which correspond to the original variable u = eu,.. From the second line of table 1 it is apparent that for u the corresponding x t < 0.007612. Denote by 2 c' the cone in the plane l a u l , al0 - ~ ' = { (la~, I, a,0)l la~01/lau I -< 0.1, 0 < la~ll - 0.008}.
(4.34)
Consider a larger cone . ~ = {( lal~ I , alo ) ] ]al01/1all [ _< 0.2, 0 < ]a H I < 0.009}.
(4.35)
We wish to prove that ,~c' x S 1 is included in the image q~(.~ x $1). In order to estimate IAqh I and Ia~21 in _~ include it in a forward invariant domain .~ of the type P1P2P3P4P5 as shown in fig. 4 with alomax = 0.001,
alorain = - 0.001,
allma x = 0.01.
(4.36)
Clearly the above numbers satisfy (4.1). Also, by (4.4), the laul coordinate of the point P2 is bigger than 0.01 x e -°'°°1 > 0.009 and hence _~c c_ 2 . The bounds on IA~II and IA~21 in .~ a p p e a r in the last line of table 1. They are E1 = 0.4782 x 10 -3 and E2 = 0.1659 x 10 -3. Now consider the conical domains e2~, with 0 < e < 1. We claim that the values of IA~ll and IA~o21 in e2~ are bounded by e . E1 and e . E2 correspondingly. Indeed, in order to estimate Aqh,2 on e . ~ we should include it into domain 2~ similar to .~ but with the bounds al0raax = 0.001e,
al0mi n = - 0 . 0 0 1 e ,
a l l m xa = 0.01e.
(4.37)
Clearly 5~, contains e-~c. Since al0mi n = -al0raax, our algorithm of computing E1 and E2 depends only on the two p a r a m e t e r s ~o = (al0max, all max)" Let us prove that E l ( e ~ 0) _< e • E l ( ~ o ) ,
E 2 ( e ~ 0) _< e • E 2 ( ~ o ) .
(4.38)
It is clear that if f ( $ ) is a polynomial with positive coefficients which depends on a positive vector $, then e-~f(eYc)
D. Michelson ~Elementary particles as solutions of the Sivashinsky equation
535
scaling and have positive coefficients. This finally proves (4.38). Hence, for (latll, a l 0 ) ~ .~c we have
]A~ol(la111, alo, a)I _.<
0.4782 × 10-alaHI < 0.05321a111, 0.009
IA~2(lalll,alo,a) I
0.1659 X lO-3lall[ < 0.01851aHI. 0.009
(4.39)
If la~0l/lanl = 0.1, then (lal0[ - IA~o21)/(la111 + IA~0~l)>0.172 and hence the image of the boundary a.~¢ under the map ~0 stays apart from the boundary of .~¢'. By (4.39), condition (3.118) holds in .~¢. This proves our claim that ~P(-~c × $1) -~ ~¢' × $1 and establishes theorem 1 also for [y'(0)l < 0.005. Of course, the hardest part of theorem 1 is the interval computation of the trajectories for r up to r~ = 80, which will be discussed in section 6. Let us remark that the program ASYMP is not written in interval arithmetic. However, all numbers in our formulas are positive and hence there is no loss of significance. Nevertheless to be sure we multiplied all expressions in the main loops by a factor (1 + K- DEL), where D E L -- 10-~3 bounds the maximal relative error in our computer and K is the number of operations in the expression. By the way, this correction did not change any digit in table 1. Thus the results in table 1 should be considered as mathematically rigorous.
5. Analytical proof of global existence for small ly'(0)l In this section we wish to prove theorem 2. Notice that the function 3 cos r Yo =
~
3 sin r +
(5.1)
r2
is the solution of the linear equation
with the initial condition y(0) = y " ( 0 ) -- 0,
y ' ( 0 ) = 1.
(5.3)
This solution is found by the change of variables y - u / r so that u satisfies ( 0 2 + 1)(O + r - l ) u = O.
(5.4)
Hence (D + r - 1 ) u = a cos r + b sin r, and formula (5.1) is implied by the initial conditions (5.3). Now consider eq. (1.5) with initial condition y(0) = y"(0) = 0, y'(0) = e. Rescale the problem as in (4.30)-(4.32). For bounded r problem (4.31)-(4.32) is an ~ ( e ) perturbation of problem (5.2)-(5.3). Fix r ffi r t and consider the vectors
(5.5)
536
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
at r -- r 1. Here u~ = rye, u 0 = ry 0. Clearly with small e the difference $~ - 20 could be made arbitrarily small. By (5.1) (5.6)
2 0 = 3(111 + C l ( r i - l ) l , ~2(ri-2)) where ] ~ , ( r l ' ) l < r i - ' + 2 r 1 2 + 2 r i -3,
t~2(ri-l)[ < 3 ( 3 r i - 2 + 2 r i - 3 ) .
(5.7)
For r 1 = 80 the slope Ix02/x011 < 0.001 and hence for lel < 0.005 vector (le[xm, EX02) belongs to the cone , ~ ' in (4.34). Thus for lel sufficiently small, the vector (lelx~l, ex~2) also belongs to -~c'. We have shown in section 4, however, that -~c' lies in the domain of existence of the asymptotic solution. This proves theorem 2. Actually, we do not need here the whole machinery of section 4. Since the claim of theorem 2 is local, we may let r 1 be arbitrarily large at the expense of very small [el. The point (lelx~l,ex~2) is made to belong to the cone _~'. It is clear that the bounds E1 and E2 on IA~IL, IA~21 in the domain _~ defined by (4.36) tend to 0 as r 1 ~ ~. Hence for large rx, the inequalities in (4.39) will be surely fulfilled so that ~o(_~c× S') _~ 2 c' × S'.
6. Numerical solution of the differential equation in interval arithmetic
Consider an initial value problem for a system of differential equations Y' = f ( Y ) ,
Y(r0) =Y0,
Y~Ra-
(6.1)
The idea of how to solve this problem with rigorous bounds in interval arithmetic was explained in ref. [1] but for completeness we present it here. Given y(r) we compute y(r + h) by the Taylor formula y(n)(r) h n
y(r+h)=y(r)+y'(r)h+...+
n!
+E~
(6.2)
l<_i
(6.3)
with the remainder
En'i=
y}n+l)(r + Oih ) h n+l (n+l)! , 0
The coefficients y(k)(r), 1 _
m
(6.4)
D. Michelson ~Elementary particles as solutions of the Siuashinsky equation
537
where F A C T is a "safety" factor, e.g. F A C T = 1.5, while m < n, e.g. n = 8, m = 5. The box Yh0 is about to include all values {yi(r + Oih)}/d=l, 0 < 0 i < 1 of the true solution. (3) Call D E R I V with the interval [r, r + h] instead of r and the box Yh0 instead of y. The resulting box values of the derivatives are denoted by .,,~k) rh0 • (4) Compute the box
Yh = ~ Ytk)(r) [0, h]k/k! +y~h~+l)[0, h]~+l/(n + 1)[
(6.5)
k=O
and check whether Yh lies in Yh0" If not, the program fails. If yes, go to step (5). (5) Compute y(r + h) by (6.2) with E n replaced by the box y~h~o+l)hn+l/(n+ 1)!. The test in (4) shows d that the derivative {y/~n+l)(r + Oih)n=l}, 0 < 0 i < 1 of the true solution lies in yh0"t~+l).Hence y(r + h) contains the true solution at r + h. Note that the above algorithms requires only two calls to the most time-consuming subroutine DERIV. Suppose now that the true solution y ( r ) tends to a critical point, say y = 0. Then for large r
y' =Ay +g(y),
g(y) = ~ ( y 2 ) .
(6.6)
The problem with the interval arithmetic is that it results in exponential growth of the solution even when A has no positive eigenvalues. For example, if y is a box with the center at 0, then Ay = abs(A) y, where abs(A)---{[au[}. In case of eq. (1.5), the matrix A has the eigenvalues A -- + i and h = 0. The interval arithmetic will transform A = + i into h -- 1. Clearly, with such exponential growth we cannot advance the solution of (1.5) till r - - 80. Here we suggest a simple remedy to this problem. Assume for simplicity that A is diagonal with pure imaginary eigenvalues A1, h 2 , . . . , Ak. Change variables y = z e "4",
z ' = e-'4" g ( y )
=
e-Arg(ze Ar) =
~(z2).
(6.7)
Actually, instead of solving equation y' =Ay we now compute in interval arithmetic the imaginary exponent e Ar. The latter is a periodic function and hence there is no accumulation of errors. The interval arithmetic causes instability also in the equation z' = ~'(z z) but to a much lesser extent since Izl is small. The computations are now carried out in terms of the z variable, i.e. given z(r) we compute the derivatives z
m
y(r;p)=
~
3,~,vk.r /
,',x
k----!
3~ , m y
.
"
m
~pktr;u)p / K . + ~ p m ( r ; p ) p /m!.
(6.8)
k=O
For eq. (1.5) we used m = 2, while p was the initial slope y'(0). The derivatives O~y/apk satisfy differential equations with respect to r which are solved in the same way as (6.1) above. The equation for t h e ruth derivative of y depends on the interval p = [ - 8 , 8]. However, since Omy/Op'' is multiplied by a small interval p " = [ - 8 m, 8 m] the resulting y(r, p) is not affected too much. This way we could afford 8 = 5 x 10 -3 in eq. (1.5) and still reach r = 80.
D. Michelson /Elementary panicles as solutions of the Sivashinsky equation
538
It turns out that eq. (1.5) in z variables as in (6.7) possesses a non-linear stability. Actually, to th leading order, z is the vector all, al0. As we saw, system (1.10) has a stable critical point (0,0). Fc example, the second equation in (1.10) to the leading order is a~0= - a l 0 . Suppose that al0 is symmetric interval al0 = [ - e , e]. If one computes alo(r + h) = alo(r) + halo(r) = alo(r) - halo(r) th result is [ - e , e ] - h [ - e , e ] = [ - ( l + h ) e , ( l + h ) e ] , i.e. the interval is amplified! Instead comput a l o ( r ) - - h a l o ( r ) = a l o ( r ) ( 1 - h ) = [ - ( 1 - h ) e , ( 1 - h ) e ] . This idea applies also to the first equation i (1.10) if al0 is negative, especially as we solve the equation for the second derivative zpp. Indeed, becaus of the interval value of p, zpp becomes a large almost symmetric interval. Equation all = aHa~o for zp becomes a~lpp = a~xppa~o+ allalopp. Since a~0 turns out to be negative, one should compute aHp p h(alxppalo) in the form allpv(1 + halo). The above ideas are implemented in eq. (1.5) as follows. First switch to the variable u =yr whic satisfies the differential equation
u"+u'= -(u"+u)r-l-f,
f = 2 u ' ( r - l ) ' + u ( r - X ) " +u2r-1.
(6.9
Then define the vector function ~ = (a~, a2, a3) by al=
-u'sinr-u"cosr,
a2=u'cosr-u"sinr
,
a 3 = u + u".
(6.10
For u = ax cos r + a 2 sin r + a 3 with constant a~, a2, a 3 formulas (6.10) recover the coefficients al, az, a: Clearly, vector ~ plays the role of z in (6.7). One can readily check that ~ satisfies the system of ODl
a{ = ( f +a3r-a)cosr,
a~ = ( f +a3r-l)sinr,
a'3= - a 3 r - ' - f
(6.11
with f as in (6.9). The term a3r-1 is not included in f since it enhances the stability of a 3. Namely, whe adding a 3 + ha~ one should combine a 3 - h a 3 r-1 = a3(1 - h r - l ) . The leading term in f with respect t r - l is u2r -1 = ( al cos r + a2 sin r + a3)2r -1.
(6.12
In order to additionally stabilize the first two equations in (6.11) we add a i + ha~ for i = 1,2 as follow,
al + ha~ =al(l + hC°s2r (alcosr + 2xx)) + flhcosr,
(6.13
where
x l = a 2 s i n r + a 3, f l = u ( r
--lxtt
) + 2 u ' ( r - l ) ' + ( a 3 + x 2 ) r -I
(6.14
and similarly for a 2 sin2r
h ( a 2 sin r a 2 + ha~ = a2(1 + h sin2r r
+
2X2) ) +f2 h sin r,
(6.15
where
x2=a2cosr+a3,
f z = u ( r - 1 ) " + 2 u ' ( r - ' ) ' + ( a 3 + x 2 ) r -'.
(6.16
D. Michelson ~Elementary particles as solutions of the Sivashinsky equation
539
As already mentioned
(6.17)
a 3 + ha'3 = a3(1 - h / r ) - hf.
The higher derivatives of K are computed by differentiating eq. (6.11). There are two ways of computing d k f / d r k. One is to express f in terms of ~ and then differentiate. The problem is that f will include many terms, and the differentiation will be costly. The second way is, given ~(r), find ~ ( r ) = (u(r), u'(r), u"(r)) by the formula
u=a]cosr+a2sinr+a3,
u'= -alsinr+a2cosr
, u"= -alcosr-a2sinr
(6.18)
and then compute the higher derivatives of u and f by eq. (6.9). We preferred the second way because it is computationally more efficient. The precise algorithm is as follows. For i = 0 . . . . . n compute successively the derivatives E CijU(J)u (i-j), J
(U2)(i)=
f(') = Y'~ cij [(u2)(J)(rJ
(6.19)
l ) ( i - j ) -1" u(J)(r-')(i-J+2)
+ 2u~J+l)(r-1)(i-j+l) ]
(6.20)
and U (i+3) =
--
..... ECij[U (J" + 2 ) + u{"J(r - 1 ). ( i - j ) - f " ) , J
(6.21)
where ci/= ( i ) are the binomial coefficients, Next, for i - 0 . . . . . n compute ~
/ ~(i)
--
(j).
= ~.~cija 3
J
(r
-1.(i-J)
)
,
(6.22)
a~,+ l) = y" cii( f,i, + ( a3/r)'J))(sin r)'i-i + '), J
(6.23)
a~2i + ') = ~ cij ( f~J) + ( a3/r)(/))(sin r) "-j), )
(6.24)
a~3i+,) = _f(i) _ ( a3/r ),).
(6.25)
In order not to duplicate the term f(/) + (a3/r) (J) in (6.23), (6.24), we store it under the name TERM. Formulas (6.19)-(6.25) are computed in the subroutine DERIV (see the main program R U N in the appendix). The five-step algorithm outlined above is carried out in subroutine TRAJ with vector playing the role of y in (6.1) and ff the one of z in (6.7). In step (1) DERIV computes the derivatives ~
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
540
evaluated by (6.5) and the condition ~h c ~h0 is verified. Note that we do not check whether ~h c_ ~h0Indeed, since z(n+~) is based on ~h0 it is enough to check that fih Crib0. In the final step ~(r+h) is -h0 computed by formula (6.2) with y replaced by ~ and En by ~Z("+l)hn+ ~/(n + 1)! h0 Now we incorporate in the code the dependence on the parameter p = y ' ( 0 ) . For p = [ P l , P 2 ] = [ Po - 6, Po + 8] denote by T 0 - the value of ~ which corresponds to Po, tip0- the derivative O~/Op at Po, the box value of u corresponding to the whole interval p, ~p and ~pp respectively, the derivatives Ou/~p and O2u/Op2 at p = [Pl, P2]- Similar notations are used for all primary and auxiliary variables like x~, fl, etc. For example the variable X1P0 in subroutine T R A J stands for the derivative 8x~/Op at Po. Given a general functional relation y =F(x),
(6.26)
where y and x depend on parameter p, and given x o, x, xp,,, xp and Xpp we compute
Yo=F(xo), yp=F'(x)Xp,
y=F(x),
yvo=F'(xo),
ypp=F'(x) Xpp+F"(X)Xp'Xp,
(6.27)
where F' and F" are the derivatives of F with respect to x. Actually, one can save time by computing y and yp by Y=Yo+Ypo[Pl,Pz]+Ypp[Pl,P2]2/2,
Yp = Yp,,+ Ypp[Pl, P2]
(6.28)
instead of the expressions in (6.27). This approach is used constantly in our program whenever F is a one-step function for which the derivatives F' and F" are easily calculated. Hence we store only the variables of the type Y0, Yp0 and ypp while the variables y and yp are computed by (6.28) only when it is necessary for another step, e.g. the variables U and UP in T R A J which are used later on in FPP. The above procedure is applied to all steps of the computer code. For example, the variable a3/r in (6.22) is denoted in D E R I V by A3R. Instead of (6.22) we compute A3R0(I) = EC(I,J).AD0(3,
J),RD(I-J),
A 3 R P 0 ( I ) = }~C(I, J) * ADP0(3, J ) * RD( l - J ) , A 3 R P P ( I ) = ~_,C(I,J),ADPP(3, J ) , R D ( I - J ) ,
(6.29)
where the latter D in AD and RD stands for derivative and R D ( I - J ) stands for ( r - l ) ~;-j~ The computations in T R A J are carried out for r starting at 2 until r = 80. The stepsize h for r < 0 is 1 / 8 and becomes 1 / 4 for r > 10. The order n of the Taylor series is denoted by N O R D E R and equals 8. The a'-priory box fih0 is computed by formula (6.4) with m = 5. Actually we compute in TRAJ the boxes UH00, UHOP0 and UHOPP for the variables u o, Up,, and Upp. If the test in step (4) fails, T R A J is rerun with an interval p = [Pl, P2] reduced by half. At r = 80 the results A l l S Q R = ~ 2 + a 2 and A ( 3 ) = a 3 are printed in interval form under the names I A l l l and A10, respectively. These results together with the corresponding range of y'(0) appear in table 1 and are plotted in fig. 3. The program ran on CDC Cyber 180/855 in single precision which has a 48 digit binary mantissa. The CPU time for one interval y'(0) ~ [pl, P2] was about 30 s, so that table 1 took about an hour of CPU. Some of the runs resulted in too large boxes in fig. 3, so that the corresponding intervals [p~, P2] were halved and T R A J was rerun
D. Michelson / Elementaryparticles as solutions of the Sivashinsky equation
541
twice. Hence, in order to obtain table 1 in a single run, one has to use end-points of the intervals y'(0) which appear in table 1 as an input to TRAJ. The normalized box for lY'(0)I < 0.005 which appears in the first line of table 1 was computed by a modification of TRAJ for eq. (4.31). Note that r~ --- 80 is close to the minimal value of r 1 for which the program ASYMP could prove convergence. Also the order of the Taylor method may not be reduced below 8 and h below 1/4 unless the error E, becomes too large and the computation explodes. It is remarkable that all the special efforts like the change of coordinates in (6.1), the expansion with respect to y'(0), the addition formulas (6.13)-(6.17) were necessary to do the computation with a reasonable amount of computer time. The computations for the program RUN were done in interval arithmetic. We are aware of the existing compiler developed by U. Kulisch and W. Miranker and co-workers to handle interval arithmetic. Unfortunately this compiler was not available to us in the CDC system. Instead, we wrote a library INTAR (which stands for interval arithmetic) of elementary interval operations SUM, DIF, MUL and DIV which replace the usual operations of summation, difference, multiplication and division. A preprocessor COMPINT (compiling intervals) was designed in order to split all the algebraic Fortran expressions (indicated by letter I) into a sequence of simple CALL statements to the above subroutines. More about the library INTAR, the preprocessor and the final translated program will be disclosed in the next section.
7. The interval arithmetic programming As mentioned in section 6 we wrote a library INTAR of elementary arithmetic operations SUM/j, DIFij, MULij, DIVij, where i and j are integers 1 or 2. For the commutative operations SUM and MUL, always i < j. The calling sequence for SUM is CALL
(7.1)
SUMij(A,B,C),
which results in the interval sum C = A + B. If A is a point value then i = 1, if A is an interval i = 2. Similarly, j relates to the type of B. In any case A, B and C are real and C is always an interval. What we said about the type of variables applies to all operations. The interval variables should be declared as arrays in the calling program. It is also assumed for an interval variable A that A(1)
D(2) < C(2) < D ( 2 ) ( 1 + e),
(7.2)
where D(i)=A(i)+B(i),
i=1,2,
e = 2 - 4 7 + 2 -48.
(7.3)
Namely, we add A ( i ) + B(i) by the computer with symmetric rounding and obtain SUM(i), i = 1, 2. Then C(1) and C(2) are defined by C(1) = SUM(l) - 2 el-48,
C(2) = SUM(2) + 2 e2-48
(7.4)
where e I and e 2 are the exponents of SUM(l) and SUM(2) in the binary floating point representation.
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
542
Number 48 here stands for 48 digit binary mantissa in our computer. The rounding in (7.4) may cause a relative error up to 2 - 4 7 with additional maximal relative error of 2 - 4 8 in the symmetric rounding. The c o n s t a n t s 2 ei-48, i = 1, 2 are formed by some logical operation at the level of the binary words SUM(l) and SUM(2). This procedure of rounding could be avoided would the computer have an option of directional rounding, i.e. either up or down. The other operations are defined in a similar way. For example, MUL22(A, B, C) forms an interval C(1), C(2) which satisfies (7.2) where
D(1)=
min
i 1= 1,2; i2= 1,2
{A(il).B(iz) },
D(2)=
max ,z{A(i,).B(i2) }.
i,= 1,2; i2= 1
(7.5)
Again, because of the symmetric rounding we obtain instead of D(i) machine numbers MUL(i), which are further rounded as in (7.4). The complete subroutine MUL22 is presented in the appendix to serve as an example of the interval arithmetic. Despite the many IF statements, on the average MUL22 takes three multiplications and three IF statements. In addition the library I N T A R contains the subroutines ISIN22(X,Y), ICOS22(X, Y) which for an interval X compute the interval Y_~sinX and Y _ ~ c o s X correspondingly. First we replace the interval X by a shifted interval X1 such that X](1) = X(1) + m'rr/2, Xl(2) - X I ( 1 ) = X ( 2 ) - X(1) and IXl(1)l < w/4. Now, c o s X l or sinX1 are approximated by the Taylor series up to order 16. The computation of the series and also of the interval X1 is performed in interval arithmetic utilizing the above elementary subroutines SUM, etc. Since Fortran does not allow multivalued functions, the elementary operations should be represented by subroutines. This requires each algebraic expression to be split into a sequence of elementary CALL statements. As a result, the program becomes unreadable. Instead we wrote a preprocessor C O M P I N T which does the job of translating these expressions into such sequences. The expressions to be processed are marked in the main program by I in the first column. C O M P I N T keeps track of the type of the variables - interval with type 2 and single with type 1. All variables in the interval expressions should be real since such are the arguments A, B and C in the elementary subroutines. The interval-type variables are declared in the beginning of the program by the statement I N T E R V A L followed by a sequence of names and the first line (if there are continuation lines) is marked by ID in the first two columns. All variables defined by the interval expressions enter automatically into the list of interval variables. In addition there is a working array S W O R K of interval variables which is used by C O M P I N T to store the intermediate results of computations. The most inconvenient in our program is the necessity to add to all interval variables additional first dimension of length 2. This takes care of the two end points of the interval. Because of that, the Fortran expressions become more cumbersome and the D I M E N S I O N statement very lengthy. Of course one can design the preprocessor to add automatically first dimension to interval variables but this was beyond the scope of this paper. We also do not have exponents, squares and sign inversions. Thus an interval expression of a type Y = - X + Z will be replaced by Y = 0 - X + Z and square should be written as a product. Still we believe that the listing of the program RUN with all the assisting comments could serve as a formal proof of correctness of the computations. As an example we also attach in the appendix a single expression of RUN with its translation by COMPINT. Both programs C O M P I N T and I N T A R are available upon request. All our programs were extensively tested. However, the main practical test of correctness of our programs is that eq. (1.5) was also solved by a straightforward Taylor method of the same order and same h as in T R A J but in the single-valued arithmetic. For the point value of y'(0) the results agreed with at least eight significant digits as far as R = 80. In fig. 3 one can see points exactly in the middle of the rectangles. These points correspond to the above single-valued solution with y'(O) in the middle of the corresponding interval [pl, P2]-
D. Michelson / Elementary particles as solutions of the Sivashinsky equation
543
8. Conclusion A combination of analytic and numerical methods was used to prove the existence of stationary elementary particles in the Sivashinsky equation. According to Sivashinsky, these are the classical particles. The quantum particles should oscillate in time in order to derive in a natural way the Schr6dinger equation (see ref. [3]). It is not clear how to attack the time-dependent equation (1.1), but it is possible to reduce it by the Galerkin method to an ODE system for a few time harmonics and then use a shooting method in r. Nevertheless, the suggested method of interval solutions by computer could have a very broad application in mathematical problems. For example, the matching of asymptotic expansions for differential equations could be justified by this method. This paper shows that with all the power of modern computers the solution, if possible, may require delicate stabilization techniques and good programming skills. In order not to duplicate efforts, standard libraries of tested and efficient programs in interval arithmetic should be developed. These programs will be declared theorems or lemmas and will be quoted by other programs in the same way as already proved mathematical results are quoted in other papers. This may change the way papers are published - instead of printed articles, they should appear as an electronic file on a diskette or sent to the readers by electronic mail.
Appendix The Computer Programs MUL22, ASYMP & RUN d C C C C
SUBROUTINE MUL22(A,B,MUL) Given machine intervals A and B this subroutine computes a machine interval M U L w h i c h includes the product A*B. First the interval C-A*B is computed in single precision. Then interval MUL is defined by MUL(1)-C(1)-EPS1,MUL(2)-C(2)+EPS2 where EPSI-2**(EI-48),EPS22.*(E2-48) and EI,E2 are the exponents of C(I),C(2) respectively REAL MUL BOOLEAN EXPI,EXP2,EPSO DIMENSION A(2),B(2),MUL(2) EQUIVALENCE (CIPOS,EXPI),(C2POS,EXP2),(EPS,EPSO) IF (A(1)) 60,10,10 i0 IF (B(1)) 30,20,20 20 CI=A(1)*B(1) C2-A(2)*B(2) GO TO 170 30 C1-A(2)*B(1) IF (B(2)) 50,40,40 40 C2=A(2)*B(2) GO TO 170 50 C2=A(1)*B(2) GO TO 170 60 IF (A(2)) 70.70,120 70 IF (B(1)) 90,80,80 80 CI=A(1)*B(2) C2=A(2)*B(1) GO TO 170 90 C2-A(1)*B(1) IF (B(2)) 110,100,100 i00 CI-A(1)*B(2) GO TO 170 110 C1-A(2)*B(2) GO TO 170 120 IF (B(1)) 140,130,130 130 CI=A(1)*B(2) C2=A(2)*B(2) GO TO 170 140 IF (B(2)) 150,150,160
D. Michelson /Elementary particles as solutions of the Siuashinsky equation
544
150 CI=B(1)*A(2) C2=B(1)*A(1) GO TO 170 160 CI=AMINI(A(2)*B(1),A(1)*B(2)) C2=AMAXI(A(2)*B(2),A(1)*B(1)) 170 IF(CI.EQ.0E0) GO TO 190 CIPOS=ABS(CI) EPSO=(EXPI.AND.O"37776000000000000000")-O"00570000000000000000" MUL(1)=CI-EPS 180 IF(C2.EQ.0E0) GO TO 200 C2POS=ABS(C2) EPSO=(EXP2.AND.O"37776000000000000000")-O"00570000000000000000" MUL(2)=C2+EPS RETURN 190 MUL(1)=CI GO TO 180 200 MUL(2)=C2 RETURN END C Here
is a loop of the original program TRAJ
************************************************************************
I I
DO 52 I-0,2 U(I,I)=U0(I,I)+DAI(1)*UP0(I,I)+DAI(1)*DAI(1)/2.E0*UPP(I,I) UP(I,I)=UP0(I,I)+DAI(1)*UPP(I,I) 52 CONTINUE
***********************************************************************
C Below is the corresponding
loop of the processed program
************************************************************************
DO 52 I=0,2 CALL DIV21(DAI(1),2.E0,SWORK(1)) CALL MUL22(DAI(1),UP0(I,I),SWORK(3)) CALL MUL22(DAI(1),SWORK(1),SWORK(5)) CALL MUL22(SWORK(5),UPP(I,I),SWORK(7)) CALL SUM22(U0(I,I),SWORK(3),SWORK(9)) CALL SUM22(SWORK(9),SWORK(7),U(I,I)) CALL MUL22(DAI(1),UPP(I,I),SWORK(1)) CALL SUM22(UP0(I,I),SWORK(1),UP(I,I)) 52 CONTINUE PROGRAM ASYMP C This program estimates from above the coefficients A(N,K,0) of the C asymptotic solution Y-asymp. and their derivatives A(N,K,L) up to C order NMAX and also estimates from above the truncation error in C the norm SQRT(U' **2+U' '**2)/2.+ABS(U+U' ' ). C The upper bound AIIMAX on All and lower and upper bounds A I O M I N , A I O M A X C on AI0 should be prescribed. The resulting truncation error in the above C norm is denoted by ERROR. The parts of E R R O R w h i c h correspond respectively C to the first and the second term in the norm are called ERROR1 and ERROR2. REAL LYAS COMMON/AS/A(-2 :18,-18 :18,0 :20) ,CB(0:18,0: 18) DIMENSION B(18,0:18,0:18),Y(18),YR(18),X(4),Yl(4) R~80. NMAX=6 DEL-I. E-13 C DEL is the maximal relative round-off error (our computer has 48 digit C binary mantissa). The factor (I.+K*DEL) enters all expressions where C the round off error could be accumilated. Since all terms are positive C the computed results are always greater than the true ones. C Computation of the binomial coefficients CB(0,0)=I. DO 3 I=I,2*NMAX
CB(I,0)=I. CB(I, I ) - l . DO 3 J - l , I - 1
D. Michelson / Elementary particles as solutions of the Sivashinsky equatio.
C C C C C
C
C
C C
3 CB(I,J)-CB(I-1,J-1)+CB(I-1,J) DO 5 N--2,2*NMAX DO 5 L-0,2*NMAX+2 DO 5 K--2*NMAX,2*NMAX 5 A(N,K,L)-0. AIIMAX-0.23 AIOMAX-0.03 AIOMIN--I.0E0 A(I,0,0)-AMAXI(ABS(AIOMAX),ABS(AIOMIN)) A(1,1,0)-AIlMAX A(1,-1,0)-A(1,1,0) A(2,0,0)-0. A(2,1,0)-A(1,1,0)*(1.+2.*A10MAX+19./6.*A(1,1,0)**2) A(2,-I,0)-A(2,1,0) YMAX-AIOMAX+AIOMAX**2 YMIN-AIOMIN+AIOMIN-*2 IF(A10MIN.LT.-.5) YMIN--.25 A(1,0,1)-AMAXI(YMAX+2.*AIlMAX**2,ABS(YMIN)) A(I,0,1) is the maximum of ABS(A10+A10**2+2*All**2) in the dumaln AIOMIN<-AI0<-AIOMAX, 0<-AII<-AIIMAX. It is assumed AIOMIN(-0<-AIOMAX Note that ABS(A10+AI0**2)-<.5**2 A(1,1,1)-A(1,0,0)*A(1,1,0) A(1,-1,1)-A(1,1,1) Find max. of AI0+3*AI0**2+2*AI0**3 for AIOMIN<-AI0<-AIOMAX. This expression enters AI0'' X(1)--(3.+SORT(3.))/6. X(2)-(-3.+SQRT(3.))/6. X(3)-AIOMIN X(4)-A10MAX DO 6 I-1,4 6 YI(I)-X(I)+3.*X(I)**2+2.*X(I)**3 YMAX-AMAXI(ABS(YI(3)),ABS(YI(4))) IF(X(1).GT.X(3).AND.X(1).LE.X(4)) YMAX-AMAXI(YMAX,ABS(YI(1))) IF(X(2).GT.X(3).AND.X(2).LE.X(4)) YMAX-AMAXI(YMAX,ABS(YI(2))) A(1,0,2)-(YMAX+2.*AIlMAX**2*(1.+A10MAX))*(1.+10.*DEL) A(1,1,2)-AIIMAX*AMAXI(A10MAX+2.*AIlMAX**2,-A10MIN)*(1.+5.*DEL) A(I,-I,2)-A(I,I,2) DO 30 N-2,2*NMAX IF(N.LT.4.0R.N.GT.NMAX) GO TO i0 SUM-SQRT((COEF(N,0,0)**2+COEF(N,1,0)**2)*(1.+2.*DEL)) N0-MOD(N,2) IF(N0.EQ.I) GO TO 7 A(N-I,0,0)-SUM/(N-3-2.*AIOMAX)*(I.+DEL) A(N-I,I,0)-SUM/(N-3-2.*A10MAX)/2.*(1.+2.*DEL) GO TO 8 7 A(N-I,I,0)-COEF(N,I,0)/(N-2+AIOMIN)/2~*(I.+2.*DEL) These estimates come from the differential equations. For even N, A(N,0,0)-0 8 A(N-1,-1,0)-A(N-I, i, 0) I0 CONTINUE DO 30 L-0,2*NMAX-N IF(N.GT.NMAX.AND.L.GT.0) GO TO 30 No derivatives are needed for the truncation term L(Y-asympt) DO 30 K-0,N IF(K.NE.0.OR.N0.NE.I) GO TO 15 B(N,K,L)-0. GO TO 30 For odd N, B(N,0,L)-0 and A(N-1,0,L)-0 15 B(N,K,L)-COEF(N,K,L) IF(N.GT.NMAX) GO TO 30 No computation of coefficients IF(K.LE.I) GO TO 20 A(N,K,L)-(B(N,K,L)/K/(K*K-I))*(I.+2.*DEL) A(N,-K,L)-A(N,K,L) GO TO 3O
545
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
546
20 IF(N.EQ.2.AND.L.LE.I) GO TO 30 C A(I,K,I),A(I,K,2) were already computed A(N-I,K,L+I)=(B(N,K,L)/IABS(3*K*K-I))*(I.+DEL) A(N-I,-K,L+I)=A(N-I,K,L+I) IF(N.NE.2.0R.K.NE.0) GO TO 30 A(I,0,L+I)-(A(I,0,L+I)-2.*(A(I,0,0)-AIOMAX)*A(I,0,L) I*(I.-2.*DEL))*(I.+DEL) C There is a cancellation in the (L+l)-th derivative of AI0+AI0**2 C which diminishes the previously computed A(I,0,L+I) 30 CONTINUE LYAS-0. C LYAS=L(Y-asymptotic) DUAS=0. C DUAS=DeI-U-asympt. DO 42 N=NMAX+I,2*NMAX SUM=0. DO 40 K=I,N 40 SUM-(SUM+B(N,K,0))*(I.+DEL) SUM=(2.*SUM+B(N,0,0))*(I.+2.*DEL) 42 LYAS-(LYAS+SUM/(R*(I.-DEL))**N)*(I.~DEL) DO 46 N=2,NMAX SUM-0. DO 44 K-I,N 44 SUM-(SUM+A(N,K,0))*(I.+DEL) SUM=(2.*SUM+A(N,0,0))*(I.+2.*DEL) 46 DUAS=(DUAS+SUM/(R*(I.-DEL))**(N-I))*(I.+DEL) UASMAX=DUAS+2.*A(I,I,0)+AIOMAX UASMIN=-DUAS-2.*A(I,I,0)+AIOMIN UAS-AMAXI(UASMAX,ABS(UASMIN)) IF(UASMAX.GT.(NMAX-3)/2.0R.UASMIN.LT.-(NMAX-2)/2.) STOP C Condition (3.84) has been violated
CI=SQRT(5.+4.*SQRT(2.)*UAS+8.*UAS**2)+4.*SQRT(2.)*(NMAX-I+UAS)/R C See formula (3.103) C2=2. IF(NMAX.LT.6) GO TO 48 IF(UASMAX.GT.I.0.OR.UASMIN.LT.-3./2.) GO TO 48 C Check condition (3.88)
C2-I./2./SQRT(NMAX-5.)
C
C
C
C
C
48 C4-1.-CI*C2/R C3-C4+SQRT(C4**2-16.*C2**2*LYAS*R**2) ZL2=2.*SQRT(2.)*C2/C3*LYAS*R**NMAX see formulas (3.105)-(3.107) and (4.26) ZI-ZL2 Z2=ZL2/2. Z-SQRT(2.)*ZL2 Z-ABS(ZI)+2.*ABS(Z2), ZL2-SQRT(ZI**2+4*Z2**2) V-Z VP-2.*Z2+(NMAX-2)*Z*(1.+I./R**2)/R VPP-2.*Z2+2.*(NMAX-2)*VP/R+V*((NMAX-2)*(NMAX-I)+I)/R**2 VVPP-ZI+2.*(NMAX-2)*VP/R+(NMAX-2)*(NMAX-1)*V/R**2 VVPP-V+V'' ; see formulas (3.91),(3.93),(3.95),(3.96) ERII-(SQRT(VP**2+VPP**2)+(NMAX-2)*SQRT(V**2+ I(2.*VP+(NMAX-I)*V/R)**2)/R)/R**(~I4AX-2)/2. ERI2-(VVPP+(NMAX-2)*(V+2.*VP+(NMAX-I)*V/R)/R)/R**(NMAX-2) ERI-ERII+ERI2 see formulas (3.114),(3.115) UP-0. UPP=0. UUPPm0. UP-U',UPP-U'',UUPP-U+U'' where U stands for del-U-asymp, in (3.113) DO 52 N=I,NMAX SUMP-0. SUMPP-0. SUMUP-0. Ii=l IF(N.EQ.I) Ii=0 DO 54 K=-N,N
1). Michelson / Elementary particles as solutions of the Sivashinsky equation C Ii-0 removes the zero order terms which come from N-I SUMP~(SUMP+(II*ABS(K)+(N-1)/R)*A(N,K,0)+A(N,K,1)/R)*(1.+5.*DEL) SUMPP-SUMPP+II*K*K*A(N,K,0)+2.*ABS(K)*((N-1)*A(N,K,0)+ 1A(N,K,1))/R+(A(N,K,2)+2.*(N-1)*A(N,K,1)+A(N,K,0)*(N-1)*N)/R**2 SUMPP=SUMPP*(I.+10.*DEL) SUMUP-SUMUP+II*ABS(K*K-1)*A(N,K,0)+2.*ABS(K)*((N-1)*A(N,K,0) I+A(N,K,1))/R+(A(N,K,2)+2.*(N-1)*A(N,K,1)+A(N,K,0)*(N-1)*N)/R**2 SUMUP-SUMUP*(1.+I0*DEL) 54 CONTINUE UP-(UP+SUMP/(R*(1.-DEL))**(N-1))*(1.+DEL) UPP=(UPP+SUMPP/(R*(I.-DEL))**(N-1))*(1.+DEL) UUPP-(UUPP+SUMUP/(R*(I.-DEL))**(N-!))*(I.+DEL ) 52 CONTINUE ER21-SQRT(UP**2+UPP**2)/2. ER22=UUPP ER2-ER21+ER22 C ER2 comes from the asymptotic expansion Y-asymp. ERRORIEERII+ER21 ERROR2-ERI2+ER22 ERROR=ERI+ER2 PRINT 55, NMAX,R,AIOMIN,AIOMAX,AIIMAX 55 FORMAT(' NMAX~',I2,' R=',F6.1, ' AIOMIN=',F6.2,'AIOMAX-', IF6.2,' AIIMAX-',F6.4) PRINT 56, ERRORI,ERROR2,ERROR 56 FORMAT('EI-',EI2.4,' E2-',EI2.4,' E-',EI2.4) WRITE(I,57) R,NMAX,AIOMIN,AIOMAX,AIIMAX,ERRORI,ERROR2,ERROR 57 FORMAT('R-',F3.0,' N-',II,IX,F5.2,'(-AI0<-',F3.2,' 0<-All<-', IF5.4,' EI-',E9.4,' E2-',E9.4,' E-',E9.4) STOP END FUNCTION COEF(N,K,L) COMMON/AS/A(-2:I8,-18:I8,0:20),CB(0:I8,0:I8) DEL-I.E-13 IF(K.NE.0) GO TO 5 M-N/2 N0-N-M*2 IF(N0.EQ.0) GO TO 5 COEF-0. C See (3.24) and (3.29) RETURN 5 SUM-0. DO I0 NI-I,N-I DO i0 K1--N1,N1 DO i0 L1-0,L N2-N-NI K2-K-KI L2-L-LI IF(IABS(K2).GT.N2) GO TO i0 SUM-(SUM+CB(L,L1)*A(N1,K1,L1)*A(N2,K2,L2)*(1.+2.*DEL))*(1.+DEL) i0 CONTINUE COEF-K*IABS(K*K-I)*A(N,K,L)+IABS(3*K*K-I)*A(N-I,K,L+I)+ IIABS(K*K*(3*N-7)-N+3)*A(N-1,K,L)+K*(3.*A(N-2,K,L+2)+ 2IABS(6*N-17)*A(N-2,K,L+l)+IABS(3*N*N-17*N+22)*A(N-2,K,L)) 3+A(N-3,K,L+3)+IABS(-3*N+I0)*A(N-3,K,L+2)+IABS((N-2)*(N-3)+ 4(N-2)*(N-5)+(N-3)*(N-5))*A(N-3,K,L+I)+IABS((N-2)*(N-3)*(N-5))* 5A(N-3,K,L)+SUM COEF-(1.+I2.*DEL)*COEF RETURN END
547
548
D. Michelson /Elementary particles as solutions of the Sil~ashinsky equation P R O G R A M RUN (INPUT, OUTPUT, TAPE1 ) D I M E N S I O N A1 (2 ) ,DAI (2 ) ,SWORK (i00 ) ID INTERVAL SWORK AI0 = .5E-2 D E L A I - 1. E-2 AIMAX-. 34E0 D A I M I N - 1. E-4 i0 A I I - A I 0 + D E L A I A1 ( 1 )- (AI0+AII)/2. E0 DAI ( 1 )= (AII-AI0)/2. E0 DAI (1 )=-DAI ( 2 ) IF(AII.GT.AIMAX) GO TO 30 C A L L TRAJ(AI,DAI, IER) IF(IER.EQ.0) G O TO 20 DELAI~D~I/2 .E0 I F ( D E L A I . L T . D A I M I N ) GO TO 50 G O TO i0 20 A I 0 - A I I GO TO i0 30 WRITE(I, 40) 4O F O R M A T ( ' T H E RUN S U C C E S S F U L L Y COMPLETED') STOP 50 W R I T E ( I , 6 0 ) A I 0 , D E Z ~ I 60 FORMAT( ' RUN S T O P P E D AT AI0-' ,E12.4, ' BECAUSE DELAI-' ,E12.4, i' WAS TOO SMALL') STOP END S U B R O U T I N E TRAJ (A1, DAI, IER ) ID I N T E R V A L A1, DAI, SINR, COSR, SINRI, COSRI, SINRH, COSRH, R, RDH, H, H0, SWORK ID I N T E R V A L Y0, YP0, YPP, U0, U, UP0, UP, UPP,A0, A, AP0, AP, APP, UH0, UHP0, IUHPP, F0, FP0, FPP, ADO, ADP0, ADPP, ADER0, ADERP0, ADERPP, UD0, UDP0, UDPP, 2UDER0, UDERP0, UDERPP, UEPS0, UEPSP0, UEPSPP DIMENSION IB(0:20,0:20),AI(2),DAI(2),H(2),H0(2),SINR(2),COSR(2), ISINRI (2), COSRI ( 2), SINRH(2), COSRH(2), 3R(2), RI(2), RD(2,0 :20), RDI(2,0 :20), RDH(2,0 :20), HJ(2), 4AI ISQ ( 2 ), A1 ISQI ( 2 ), AIISQ2 (2 ), A I I S Q R (2 ), SWORK (I00 ) DIMENSION Y0(2,0:2),YP0(2,0:2),YPP(2,0:2),U0(2,0=3),UP0(2,0:3), IUPP(2,0:3) ,A0(2,3) ,A(2,3) ,AP0(2,3) ,AP(2,3) ,APP(2,3) ,UH0(2,0 :2), 2UHP0(2,0: 2) ,UHPP(2,0: 2) ,F0(2) ,FP0(2) ,FPP(2) ,AD0(2,3,0: 20), 3ADPO(2,3,0:20),ADPP(2,3,0:20),ADERO(2,3,0:20),ADERPO(2,3,0:20), 4ADERPP( 2,3,0:20), UD0(2,0 :20), UDP0(2,0 :20) ,UDPP(2,0 :20), 5UDERO(2,0:20),UDERPO(2,0:20),UDERPP(2,0:20),AHO(2,3),AHPO(2,3), 6AHPP(2,3), UH00(2,0 :2) ,UHOP0(2,0 :2) ,UHOPP(2,0 :2) ,U(2,0: 3) ,UP(2,0 :3) DIMENSION XI0(2),XI(2),XIP0(2),XIP(2),XIPP(2),FI0(2),FIP0(2) I,FIPP(2) ,X20(2) ,X2(2) ,X2P0(2) ,X2P(2) ,X2PP(2) ,F20(2), 2F2P0 (2), F2PP(2) ,U02(2) ,UEPS0 (2), UEPSP0 (2) ,UEPSPP(2) ,HRC(2) C O M M O N / T / C B (0 :20,0 :20 ) U P P M A X = 2. E5 C U P P M A X is the m a x i m a l a l l o w ed value of the d e r i v a t i v e UPP. If exceeded, C the d i a g n o s t i c s IER-I is produced, the p r o g r a m returns to RUN and the C increment of A1 is halved. IER=0 NORDER = 8 C This is the order of the Taylor m e t h o d C C o m p u t e the b i n o m i a l c o e f f i c i e n t s and store them in the array CB C The integer array IB is used to avoid round-off errors IB(0,0)=I CB(0,0)=I.E0 DO I0 I = I , N O R D E R + I IB(I,0)=I
IB(I,I)-I CB(I,0)=I.E0 CB(I,I)=I.E0 DO i0 J - l , I - i IB(I,J)=IB(I-I,J-I)+IB(I-I,J) i0 C B ( I , J ) = I B ( I , J ) R0=2.E0
D. Michelson / Elementary particles as solutions of the Sivashinaky equation
C
C C C C C I I I I I I I I I C I I C C C C I I I I I I I I I C C C I I
H(1)-0.125EO H(2)-H(1) RMAX-80.0E0 R0 is the radius at which the power expansion for Y is used R(1)-RO R(2)-R0 CALL SERIES(R,AI,DAI,Y0,YP0,YPP) Set the initial condition for the variables U0,UP0 and UPP. The letter P here and below stands for a "prime",i.e. one derivative with respect to the parameter At. The digit 0 i n t ~ e name means that A1 is fixed at the central point of the increment while the absence of 0 means that A1 is the whole increment U0(I,0)-R0*Y0(1,0) U0(1,1)=Y0(1,0)+R0*Y0(I,1) U0(I,2)-2.E0*Y0(I,I)+R0*Y0(1,2) UP0(1,0)-R0*YP0(1,0) UP0(1,1)-YP0(1,0)+R0*YP0(1,1) UP0(1,2)-2.E0*YP0(I,I)+R0*YP0(I,2) UPP(I,0)-R0*YPP(1,0) UPP(I,1)-YPP(1,0)+R0*YPP(I,1) UPP(I,2)-2.E0*YPP(I,I)+R0*YPP(I,2) Below is the computation of the derivatives of 1/R RD(I,0)-I.E0/R(1) DO 20 I=I,NORDER+2 EXP=-I RD(I,I)=EXP*RD(I,I-I)/R(1) 20 CONTINUE ISIN22,ICOS22 are user written interval subroutines which compute SIN and COS of an interval in interval arithmetic using a Taylor series expansion up to a suitable order. The maximal possible truncation error of the series is added to the results SINR and COSR CALL ISIN22(R,SINR) CALL ICOS22(R,COSR) A0(I,I)--U0(I,I)*SINR(1)-U0(I,2)*COSR(1) A0(I,2)-U0(I,I)*COSR(1)-U0(1,2)*SINR(1) A0(1,3)=U0(I,0)+U0(I,2) AP0(I,I)=-UP0(I,I)*SINR(1)-UP0(I,2)*COSR(1) AP0(I,2)=UP0(I,I)*COSR(1)-UP0(I,2)*SINR(1) AP0(I,3)-UP0(I,O)+UP0(I,2) APP(I,I)=-UPP(I,I)*SINR(1)-UPP(I,2)*COSR(1) APP(I,2)=UPP(I,I)*COSR(1)-UPP(I,2)*SINR(1) APP(I,3)-UPP(I,0)+UPP(I,2) 30 CONTINUE Compute the point values of the derivatives UD and AD at R up to the order NORDER-I. The values of A0,AP0 AND APP are stored by DERIV in the arrays AD0(.,.,0),ADP0( .... 0) AND ADPP(.,.,0) respectively CALL DERIV(RD,DAI,SINR,COSR,U0,UP0,UPP,A0,AP0,APP, IUD0,UDP0,UDPP,AD0,ADP0,ADPP,NORDER-I,0) RI(1)~R(1)+H(1) RDI(I,0)=I.E0/RI(1) DO 40 I=I,NORDER+2 EXP--I
I 40 C The
50 C The
C are I I 52
RDI(I, I)-EXP*RDI(I, I-I)/RI(1) CONTINUE values of derivatives of 1/R in the interval (R,R+H) are computed DO 50 I-0,NORDER+2 RDH(1,I)-AMINI(RD(1,I),RDI(1,I)) RDH(2,I)-AMAXI(RD(2,I),RDI(2,I)) CONTINUE interval values of U and UP for A1 in the interval (A1-DAI,AI+DAI) computed by the corresponding Taylor forlaulas DO 52 I-0,2 U(1,I)-U0(1,I)+DAI(1)*UP0(1,I)+DAI(1)*DAI(1)/2.E0*UPP(1,I) UP(I,I)-UP0(I,I)+DAI(1)*UPP(I,I) CONTINUE H0(1)-0.E0 H0(2)-H(2)
549
550
D. Michelson /Elementary particles as solutions of the Sivashinsky equation C Compute the initial box UH0 FACT=I.5D0 DO 53 K-I,2 UEPS0(K)-0.E0 UEPSP0(K)-0.E0 UEPSPP(K)-0.E0 53 CONTINUE DO 60 I=0,2 DO 55 J-5,1,-i C The 5-th order Taylor series expansion based on the derivatives UD0 etc. C is used to estimate the size of the boxes UH00 etc. In order to assure C that U0 at the interval (R,R+H) lies in UH00 w e multiply the size of C UH00 by FACT-I.5 RJ-J I UEPS0(1)-(UEPS0(1)+UD0(I,J+I))*H0(1)/RJ I UEPSP0(1)-(UEPSP0(1)+UDP0(I,J+I))*H0(1)/RJ I UEPSPP(1)-(UEPSPP(1)+UDPP(I,J+I))*H0(1)/RJ 55 CONTINUE I UH00(I,I)-U0(I,I)+UEPS0(1)*FACT I UHOP0(I,I)-UP0(I,I)+UEPSP0(1)*FACT I UHOPP(I,I)-UPP(I,I)+UEPSPP(1)*FACT 60 CONTINUE C The above computation does not need the control of the interval arithmetic C since the boxes UH00 etc. could be computed roughly. C Compute sin and cos of R1 and of the interval (R,R+H) CALL ISIN22(RI,SINRI) CALL ICOS22(Rl,COSR1) SINRH(1)-AMINI(SINR(1),SINRI(1)) SINRH(2)-AMAXI(SINR(2),SINRI(2)) COSRH(1)-AMINI(COSR(1),COSRI(1)) COSRH(2)-AMAXI(COSR(2),COSRI(2)) C Compute the derivatives UD up to order NORDER in the interval (R,R+H) C This will provide the truncation error in the Taylor series 70 CALL DERIV(RDH,DAI,SINRH,COSRH,UH00,UH0P0,UHOPP,AH0,AHP0,AHPP, IUDER0,UDERP0,UDERPP,ADER0,ADERP0,ADERPP,NORDER,I) C New compute the interval values of U in the interval (R,R+H) DO 80 K=I,2 DO 80 I=0,2 UH0(K,I)=UDER0(K,NORDER+I+I) UHP0(K,I)=UDERP0(K,NORDER+I+I) UHPP(K,I)=UDERPP(K,NORDER+I+I) 80 CONTINUE DO 90 J=NORDER,0,-I DO 90 I-0,2 RJ-J+I I UH0(I,I)-UH0(I,I)*H0(1)/RJ+UD0(I,J+I) I UHP0(I,I)-UHP0(I,I)*H0(1)/RJ+UDP0(I,J+I) I UHPP(I,I)-UHPP(I,I)*H0(1)/RJ+UDPP(I,J+I) 90 CONTINUE C Check whether UH0 belongs to the interval box UH00 etc. DO i00 I-0,2 IF(UH0(I,I).LT.UH00(I,I)) GO TO ii0 IF(UH0(2,I).GT.UE00(2,I)) GO TO ii0 IF(UHP0(I,I).LT.UHOP0(I,I)) GO TO Ii0 IF(UHP0(2,I).GT.UHOP0(2,I)) GO TO ii0 IF(UHPP(I,I).LT.UHOPP(I,I)) GO TO ii0 IF(UHPP(2,I).GT.UHOPP(2,I)) GO TO ii0 i00 CONTINUE C The result is positive GO TO 130 Ii0 IER-I
WRITE(I,115) AI(1),DAI(1),DAI(2),R(1) 115 FORMAT(' THE BOX WAS EXCEEDED',/,'AI-',EI2.4,' i' R-',EI2.4) RETURN 130 CONTINUE C Compute the values of A at R1 DO 140 K-I,2
DAI-',2EI2.4,
D. Michelson / Elementary particles as solutions of the Sivashinsky equation
I I I I C C I I C C C C I I I C C I C I I I I I I I
I I I I
C C C I C I I I I I I I I I
DO 140 I'1,3 A0(K,I)-ADER0(K,7,NORDER+I) AP0(K,I)-ADERP0(K,I,NORDER+I) APP(K,I)-ADERPP(K,I,NORDER+I) 140 CONTINUE DO 150 J-NORDER, 2,-I RJ-J+I HJ(1)-H(1)/RJ DO 150 I-1,3 A0(I,I)-A0(1,I)*HJ(1)+AD0(I,I,J) AP0(I,I)-AP0(I,I)*HJ(1)+ADP0(I,I,J) APP(I,I)-APP(I,I)*HJ(1)+ADPP(I,I,J) 150 CONTINUE In the above loop terms in the Taylor series expansion for A0(R+H) etc. are added up starting frc~ the NORDER term till the second order term DO 152 I=1,3 A(I,I)-AD0(I,I,O)+DAI(1)*ADP0(I,I,0)+DAI(1)*DAI(1)/2.E0*ADPP(I,I,0) AP(I,I)-ADP0(I,I,0)+DAI(1)*ADPP(I,I,0) 152 CONTINUE The values of A(.,I)and AP(.,I) at R are obtained from the values of A0,AP0 and APP at R. Recall that AD0(.,0),ADP0(.,.,0) and ADPP(.,0) coincide with A0,AP0,APP at the point R. (See the comment at the first call to DERIV) DO 155 I-i,3 A0(I,I)-A0(I,I)*H(1)/2.E0 AP0(I,I)-AP0(I,I)*H(1)/2.E0 APP(I,I)-APP(I,I)*H(1)/2.E0 155 CONTINUE The above A0,AP0 and APP are multiplied below once more by H in order to obtain the second order remainder of the Taylor series X10(1)-AD0(1,2,0)*SINR(1)+AD0(I,3,0) For most trajectories Xl,Xl0 are negative XI(1)~A(I,2)*SINR(1)+A(I,3) XlP0(1)=ADP0(I,2,0)*SINR(1)+ADP0(1,3,0) XIP(1)=AP(I,R)*SINR(1)+AP(I,3) XIPP(1)=ADPP(I,2,0)*SINR(1)+ADPP(I,3,0) F10(1)-RD(1,2)*U0(1,0)+R.E0*RD(1,1)*U0(1,1)+ I(AD0(I,3,0)+XI0(1)*X10(1))/R(1) FIP0(1)-RD(I,2)*UP0(I,0)+2.E0*RD(I,1)*UP0(1,1)+ I(ADP0(I,3,0)+2.E0*XlP0(1)*Xl0(1))/R(1) FIPP(1)-RD(1,2)*UPP(1,0)+R.E0*RD(1,1)*UPP(1,1)+ I(ADPP(1,3,0)+2.E0*(XlPP(1)*Xl(1)+XlP(1)*XlP(1)))/R(1) HRC(1)-H(1)/R(1)*COSR(1)*COSR(1) A0(1,1)-AD0(1,1,0)*(1.E0+HRC(1)*(AD0(1,1,0)*COSR(1)+ 12.E0*X10(1)))+(F10(1)*COSR(1)+A0(1,1))*H(1) AP0(1,1)-ADP0(1,1,0)*(1.E0+HRC(1)*R.E0*(AD0(1,1,0)*COSR(1)+ lXl0(1)))+2.E0*HRC(1)*AD0(1,1,0)*XIP0(1)+(FIP0(1)*COSR(1) 2+AP0(1,1))*H(1) APP(1,1)-ADPP(1,1,0)*(1.EO+HRC(1)*R.EO*(A(1,1)*COSR(1)+XI(1)))
I+2.E0*HRC(1)*(AP(1,1)*(AP(1,1)*COSR(1)+2.E0*XIP(1)) 2+A(1,1)*XIPP(1))+(FIPP(1)*COSR(1)+APP(1,1))*H(1) For most trajectories A(.,1)*COSR+XI is negative. The above splitting of ADPP(.,.,I) prevents the exponential growth of the interval APP when it includes the 0 point X20(1)-AD0(I,I,0)*COSR(1)+AD0(I,3,0) For most trajectories X2,X20 are negative X2(1)-A(1,1)*COSR(1)+A(1,3) X2P0(1)-ADP0(1,1,0)*COSR(1)+ADP0(1,3,0) X2P(1)-AP(1,1)*COSR(1)+AP(1,3) X2PP(1)-ADPP(1,1,0)*COSR(1)+ADPP(1,3,0) FR0(1)-RD(1,2)*U0(1,0)+2.E0*RD(1,1)*U0(I,1)+ l(AD0(1,3,0)+X20(1)*X20(1))/R(1) F2P0(1)-RD(1,2)*UP0(1,0)+2.E0*RD(1,1)*UP0(1,1)+ I(ADP0(I,3,0)+2.E0*XRP0(1)*X20(1))/R(1) F2PP(1)-RD(1,2)*UPP(1,0)+2.E0*RD(1,1)*UPP(I,I)+ I(ADPP(1,3,0)+2.E0*(X2PP(1)*X2(1)+X2P(1)*X2P(1)))/R(1) HRC(1)-H(1)/R(1)*SINR(1)*SINR(1) A0(1,2)-AD0(I,2,0)*(1.E0+HRC(1)*(AD0(1,R,0)*SINR(1)+
551
D. Michelson ~Elementary particles as solutions of the Sivashinsky equation
552
12.E0*X20(1)))+(F20(1)*SINR(1)+A0(I,2))*H(1)
I I C I C
APO(I,2)=ADPO(I,2,0)*(I.EO+HRC(1)*2.EO*(ADO(I,2,0)*SINR(1)+ IX20(1)))+2.E0*HRC(1)*AD0(I,2,0)*X2P0(1)+(F2P0(1)*SINR(1) 2+AP0(I,2))*H(1) APP(I,2)=ADPP(I,2,0)*(I.E0+HRC(1)*2.E0*(A(I,2)*SINR(1)+X2(1))) I+2.E0*HRC(1)*(AP(I,2)*(AP(I,2)*SINR(1)+2.E0*X2P(1)) 2+A(I,2)*X2PP(1))+(F2PP(1)*SINR(1)+APP(I,2))*H(1) See the comment at the definition of APP(.,I) U02(1)=U0(I,0)*U0(I,0) U02(1)-AMAXI(0.E0,U02(1)) The above modification assures that U0**2 is positive also when U0 includes 0
I
FO(1)=(2.EO*(UO(I,I)-UO(I,O)/R(1))/R(1)-UO2(1))/R(1)
I I
FP0(1)-(2.E0*(UP0(I,I)-UP0(I,0)/R(1))/R(1)-2.E0*U0(I,0)*UP0(I,0))/R(1) FPP(1)=(2.E0*(UPP(I,I)-UPP(I,0)/R(1))/R(1) -
I
12.EO*(UPP(I,O)*U(I,O)+UP(I,O)*UP(I,O)))/R(1) AO(I,3)=ADO(I,3,0)*(I.EO-H(1)/R(1))+(FO(1)+AO(I,3))*H(1)
I AP0(I,3)-ADP0(I,3,0)*(I.E0-H(1)/R(1))+(FP0(1)+AP0(I,3))*H(1) I APP(I,3)-ADPP(I,3,0)*(I.E0-H(1)/R(1))+(FPP(1)+APP(I,3))*H(1) C The last modification causes A0(3) and APP(3) to diminish also when they C include 0 C Compute the values of U at R1
I
UO(I,O)-AO(I,I)*COSRI(1)+AO(I,2)*SINRI(1)+AO(I,3)
I
U0(I,I)=A0(I,2)*COSRI(1)-A0(I,I)*SINRI(1)
I I
UO(I,2)=-(AO(I,I)*COSRI(1)+AO(I,2)*SINRI(1)) UPO(I,O)=APO(I,I)*COSRI(1)+APO(I,2)*SINRI(1)+APO(I,3)
I I I I I
I I I I I
UP0(I,I)=AP0(I,2)*COSRI(1)-AP0(I,I)*SINRI(1) UP0(I,2)J-(AP0(I,I)*COSRI(1)+AP0(I,2)*SINRI(1)) UPP(I,0)=APP(I,I)*COSRI(1)+APP(I,2)*SINRI(1)+APP(I,3) UPP(I,I)-APP(I,2)*COSRI(1)-APP(I,I)*SINRI(1) UPP(I,2)--(APP(I,I)*COSRI(1)+APP(I,2)*SINRI(1)) R(1)-RI(1) R(2)-RI(2) UPPL2=0.E0 DO 156 I-0,2 DO 156 K-I,2 156 UPPL2-UPPL2+UPP(K,I)**2 UPPL2-SQRT(UPPL2/6.E0) IF(UPPL2.LT.UPPMAX) GO TO 158 IER-I WRITE(I,157) AI(1),DAI(1),DAI(2),R(1) 157 FORMAT(' UPP IS TOO LARGE',/,'AI-',EI2.4,' DAI-',2EI2.4, i' R=',EI2.4) RETURN 158 CONTINUE IF(R(2).GE.RMAX) GO TO 170 IF(R(2).LE.10.E0) GO TO 159 H(1)-.25E0 H(2)-H(1) 159 CONTINUE DO 160 K=I,2 DO 160 I-0,NORDER+2 160 RD(K,I)-RDI(K,I) DO 165 K-I,2 SINR(K)-SINRI(K) 165 COSR(K)=COSRI(K) GO TO 30 170 CONTINUE DO 175 I-1,3 A(I,I)-A0(I,I)+DAI(1)*AP0(I,I)+DAI(1)*DAI(1)/2.E0*APP(I,I) AP(I,I)-AP0(I,I)+DAI(1)*APP(I,I) 175 CONTINUE AIISQI(1)-A(I,I)*A(I,I) AIISQI(1)-AMAXI(0.E0,AIISQI(1)) AIISQ2(1)-A(I,2)*A(I,2) AIISQ2(1)-AMAXI(0.E0,AIISQ2(1)) AIISQ(1)=(AIISQI(1)+AIISQ2(1))/4.EO AIISQR(1)=SQRT(AIISQ(1))*(I.E0-I.E-13) AIISQR(2)-SQRT(AIISQ(2))*(I.EO+I.E-13)
D. Michelson ~Elementary particles as solutions of the Sivashinsky equation WRITE(I,180) AI(1)+DAI(1),AI(2)+DAI(2),R 180 FORMAT(6HY'(0)-,2EI5.8,' R-',2EI5.8) WRITE(I,190) A(I,3),A(2,3),AIISQR(1),AIISQR(2) 190 FORMAT('AI0ffi',2EI5.8,' ABS(AI1)-',2EI5.8) END SUBROUTINE SERIES(R,AI,DAI,Y0,YP0,YPP) C This subroutine computes the power series expansion of Y(R) up to C order Nffi80 at R DIMENSION AI(2),DAI(2),A(2,100),A0(2,100),AP0(2,100),AP(2,100), IAPP(2,100),Y0(2,0:2),YP0(2,0:2),YPP(2,0:2),SUM0(2), 2SUMP0(2),SUMPP(2),R(2),R2(2),ERROR(2),SWORK(100) ID INTERVAL AI,DAI,A,A0,AP,AP0,APP,SUM0,SUblP0,SUMPP,y0,YPO,ypp, ISWORK,R,ERROR N-80 ERROR(2)=I.E-14 ERROR(1)=-I.E-14 I R2(1)©R(1)*R(1) A0(I,I)=AI(1) A0(2,1)=AI(2) DO i0 Kffii,2 AP0(K,I)-I.E0 i0 APP(K,I)-0.E0 DO 60 I=I,N DO 30 Kffil,2 SUM0 (K) -0. E0 SUMP0 (K)-0. E0 SUMPP (K) -0. E0 30 CONTINUE I A(I,I)-A0(I,I)+DAI(1)*AP0(I,I)+DAI(1)*DAI(1)/2.E0*APP(I,I) I AP(I,I)-AP0(1,I)+DAI(1)*APP(1,I) DO 40 J-l,I-i I SUM0(1)-SUM0(1)+A0(1,J)*A0(1,I-J) I SUMP0(1)-SUMP0(1)+AP0(I,J)*A0(I,I-J)+A0(I,J)*AP0(I,I-J) I SUMPP(1)-SUMPP(1)+APP(I,J)*A(I,I-J)+2.E0*AP(I,J)*AP(I,I-J) I+A(I,J)*APP(1,I-J) 40 CONTINUE C Our interval arithmetic subroutines apply only to real numbers. Hence C integers should be first converted to reals RI21ffi2*I+l RI23-2,I*(2-I+3) I A0(I,I+I)--(A0(I,I)÷SUM0(1)/RI21)/RI23 I AP0(I,I+l)--(AP0(1,I)+SUMP0(1)/RI21)/RI23 I APP(I,I+I)m-(APP(I,I)+SUMPP(1)/RI21)/RI23 60 CONTINUE DO 70 I-0,2 DO 70 K-I,2 Y0(K,I)-0.E0 YP0(K,I)-0.E0 70 YPP(K,I)ffi0.E0 DO 80 I-N,1,-1 RII-2*I-I RI2-(2"I-1)*(2,I-2) I Y0(1,0)-Y0(1,0)*R2(1)+A0(1,I) I YP0(1,0)-YP0(I,0)*R2(1)+AP0(1,I) I YPP(I,0)-YPP(I,0)*R2(1)+APP(I,I) I Y0(I,I)-Y0(I,I)*R2(1)+A0(I,I)*RII I YP0(I,1)-YP0(1,1)*R2(1)+AP0(1,I)*RI1 I YPP(1,1)-YPP(1,1)*R2(1)+APP(1,I)*RI1 I Y0(I,2)ffiY0(I,2)*R2(1)+A0(I,I)*RI2 I YP0(I,2)-YP0(1,2)*R2(1)+AP0(I,I)*RI2 I YPP(I,2)-YPP(I,2)*R2(1)+APP(I,I)*RI2 80 CONTINUE I Y0(I,0)-Y0(I,0)*R(1) I YP0(I,0)-YP0(I,0)*R(1) I YPP(I,0)-YPP(I,0)*R(1) I Y0(I,2)-Y0(I,2)/R(1)
553
D. Michelson
554
/Elementaryparticles as solutions of the Sivashinsky equation
YP0(I,2)-YP0(I,2)/R(1) I YPP(I,2)=YPP(I,2)/R(1) I C The addition of ERROR takes care of the remainder of the series DO 90 1=0,2 Y0(I,I)=Y0(I,I)+ERROR(1) I YP0(I,I)-YP0(I,I)+ERROR(1) I I YPP(I,I)=YPP(I,I)+ERROR(1) 90 CONTINUE RETURN END SUBROUTINE DERIV(RD,DAI,SINR,COSR,UO,UPO,UPP,AO,APO,APP, IUD0,UDP0,UDPP,AD0,ADP0,ADPP,NORDER, IFLAG)
COMMON/T/CB(O:20,O:20) DIMENSION RD(2,0:20),SINR(2),COSR(2),SINDR(2,0:20),UO(2,0:3), IUPO(2,0:3),UPP(2,0:3),AO(2,3),APO(2,3),APP(2,3),UD(2,0:20), 2UD0(2,0:20),UDP(2,0:20),UDP0(2,0:20),UDPP(2,0:20),
ID
C I I
C I I I I I I I I I I
3U2D(2,0:20),U2DO(2,0:20),U2DP(2,0:20),U2DPO(2,0:20),U2DPP(2,0:20), 4FD(2,0:20),FDO(2,0:20),FDP(2,0:20),FDPO(2,0:20),FDPP(2,0:20), 5ADO(2,3,0:20),ADPO(2,3,0:20),ADPP(2,3,0:20),A3RO(2,0:20), 6A3RPO(2,0:20),A3RPP(2,0:20),DAI(2),TERM(2),SWORK(IO0) INTERVAL RD,SINR,COSR,UO,UPO,UPP,AO,APO,APP,UD,UDO,UDP,UDPO, IUDPP,U2D,U2DO,U2DP,U2DPO,U2DPP,FD,FDO,FDP,FDPO,FDPP,ADO,ADPO, 2ADPP,DAI,SWORK DO i0 Iz0,2 DO i0 J-l,2 UD0(J,I)=U0(J,I) UDP0(J,I)=UP0(J,I) UDPP(J,I)=UPP(J,I) i0 CONTINUE Compute the box values of UP AND UPD DO 15 I=0,2 UD(I,I)=UD0(I,I)+DAI(1)*UDP0(I,I)+DAI(1)*DAI(1)/2.E0*UDPP(I,I) UDP(I,I)-UDP0(I,I)+DAI(1)*UDPP(I,I) 15 CONTINUE DO 60 I-0,NORDER DO 20 J-l,2 UD0(J,I+3)-0.E0 UDP0(J,I+3)=0.E0 UDPP(J,I+3)-0.E0 U2D0(J,I)=0.E0 U2DP0(J,I)=0.E0 U2DPP(J,I)=0.E0 FD0(J,I)=0.E0 FDP0(J,I)=0.E0 FDPP(J,I)=0.E0 20 CONTINUE 11-1/2 12-I-2,Ii U2D is the derivative of U**2 DO 30 J=0,11 U2D0(I,I)=U2D0(I,I)+CB(I,J)*UD0(I,J)*UD0(I,I-J) U2DP0(1,1)=U2DP0(I,I)+CB(I,J)*(UDP0(I,J)*UD0(1,I-J)+ IUD0(I,J)*UDP0(1,I-J)) U2DPP(1,1)=U2DPP(I,I)+CB(I,J)*(UDPP(I,J)*UD(I,I-J)+ 12.E0*UDP(I,J)*UDP(1,I-J)+UD(I,J)*UDPP(1,I-J)) 30 CONTINUE U2D0(I,I)-2.E0*U2D0(I,I) U2DP0(I,I)=2.E0*U2DP0(I,I) U2DPP(1,1)=2.E0*U2DPP(I,I) IF(I2.EQ.1) GO TO 40 U2D0(I,I)=U2D0(I,I)-CB(I,II)*UD0(I,II)*UD0(I,II) U2DP0(I,I)=U2DP0(I,I)-CB(I,II)*(UDP0(I,II)*UD0(I,II)+ IUD0(I,II)*[~P0(I,II)) U2DPP(I,I)=U2DPP(I,I)-CB(I,II)*(UDPP(I,II)*UD(I,II)+ 12.E0*UDP(I,II)*UDP(I,II)+UD(I,II)*UDPP(I,II)) 40 DO 50 J=0,1 FD0(I,I)-FD0(I,I)+CB(I,J)*(U2D0(I,J)*RD(I,I-J)+UD0(I,J)*RD(I,I J+2) I+2.E0*UD0(I,J+I)*RD(I,I-J+I))
19. Michelson / Elementary panicles as solutiona of the Sivashinsky equation I I I I I I I I I I C
C
I I I I I I I I I I I I I I I
FDP0(I,I)=FDP0(I,I)+CB(I,J)*(U2DP0(I,J)*RD(I,I-J)+UDP0(I,J)~ IRD(I,I-J+2)+2.E0*UDP0(I,J+I)*RD(I,I-J+I)) FDPP(I,I)=FDPP(I,I)+CB(I,J)*(U2DPP(I,J)*RD(I,I-J)+UDPP(I,J)* IRD(I,I-J+2)+2.E0*UDPP(I,J+I)*RD(I,I-J+I)) UD0(I,I+3)=UD0(I,I+3)-CB(I,J)*(UD0(I,J+2)+UD0(I,J))*RD(I,I-J) UDP0(I,I+3)=UDP0(I,I+3)-CB(I,J)*(UDP0(I,J+2)+UDP0(I,J))*RD(I,I-J) UDPP(I,I+3)-UDPP(I,I+3)-CB(I,J)*(UDPP(I,J+2)+UDPP(I,J))*RD(I,I-J) 50 CONTINUE UD0(1,I+3)-UD0(1,I+3)-FD0(1,I)-UD0(1,I+I) UDP0(I,I+3)-UDP0(1,I+3)-FDP0(1,I)-UDP0(1,I+I) UDPP(I,I+3)-UDPP(I,I+3)-FDPP(I,I)-UDPP(1,I+I) UD(1,1+3)-UD0(1,I+3)+DAI(1)*UDP0(1,I+3)+DAI(1)*DAI(1)/2.E0*UDPP(1,I+3) UDP(I,I+3)-UDP0(1,I+3)+DAI(1)*UDPP(1,I+3) 60 CONTINUE SINDR(.,I) is the I-th derivative of SIN(R) DO 70 K-I,2 SINDR(K,0)-SINR(K) 70 SINDR(K,1)-COSR(K) SINDR(I,2)--SINR(2) SINDR(2,2)--SINR(1) SINDR(I,3)--COSR(2) SINDR(2,3)=-COSR(1) DO 80 K-I,2 DO 80 I=4,NORDER+2 80 SINDR(K,I)mSINDR(K,I-4) IF(IFLAG.EQ.I) GO TO 160 DO I00 K-I,2 DO 90 J-1,3 AD0(K,J,0)-A0(K,J) ADP0(K,J,0)-AP0(K,J) ADPP(K,J,0)-APP(K,J) 90 CONTINUE I00 CONTINUE DO 150 I=0,NORDER DO 120 K-1,2 DO 110 J-1,2 AD0(K,J,I+I)-0.E0 ADP0(K,J,I+I)-0.E0 ADPP(K,J,I+I)-0.E0 ii0 CONTINUE A3R(I) is the I-th derivative of A3/R A3R0(K,I)=0.E0 A3RP0(K,I)-0.E0 A3RPP(K,I)-0.E0 120 CONTINUE DO 130 J-0,I A3R0(I,I)-A3R0(1,I)÷CB(I,J)*AD0(1,3,J)*RD(1,I-J) A3RP0(I,I)-A3RP0(I,I)+CB(I,J)*ADP0(I,3,J)*RD(I,I-J) A3RPP(I,I)=A3RPP(I,I)+CB(I,J)*ADPP(I,3,J)*RD(I,I-J) 130 CONTINUE DO 140 J-0,I TERM(1)-CB(I,J)*(A3R0(I,J)+FD0(I,J)) AD0(I,I,I+I)-AD0(I,I,I+I)+TERM(1)*SINDR(I,I-J+I) AD0(I,2,1+I)-AD0(I,2,I+I)+TERM(1)*SINDR(I,I-J) TERM(1)=CB(I,J)*(A3RP0(I,J)+FDP0(I,J)) ADP0(I,I,I+I)-ADP0(I,I,I+I)+TERM(1)*SINDR(I,I-J+I) ADP0(1,2,I+I)-ADP0(1,2,I+I)+TERM(1)*SINDR(1,I-J) TERM(1)=CB(I,J)*(A3RPP(I,J)+FDPP(I,J)) ADPP(1,I,I+I)=ADPP(I,I,I+I)+TERM(1)*SINDR(1,I-J+I) ADPP(1,2,I+I)-ADPP(I,2,I+I)+TERM(1)*SINDR(1,I-J) 140 CONTINUE AD0(I,3,I+I)=-FD0(1,I)-A3R0(1,I) ADP0(I,3,I+I)=-FDP0(I,I)-A3RP0(I,I) ADPP(I,3,I+I)--FDPP(I,I)-A3RPP(I,I) 150 CONTINUE RETURN
555
D. Michelson /Elementary particles as solutions of the Sivashinsky equation
556
C C C C
I I I I I I I I I C
160 CONTINUE NI-NORDER+I The computation of the remainder of the Taylor series. Here the expression of A in terms of U is used unlike the previous computation of ADO based on.the differential equation for A. This way the computer time is saved although the stability is sacrificed DO 170 K-I,2 DO 170 J=l,3 AD0(K,J,NI)=0.E0 ADP0(K,J,NI)=0.E0 ADPP(K,J,NI)=0.E0 170 CONTINUE DO 180 J=0,NI NJ=NI-J AD0(I,I,NI)=AD0(I,I,NI)+CB(NI,J)*(UD0(1,J+I)*SINDR(1,NJ) I+UD0(I,J+2)*SINDR(1,NJ+I)) AD0(I,2,NI)=AD0(I,2,NI)+CB(NI,J)*(UD0(I,J+I)*SINDR(I,NJ+I) I-UD0(1,J+2)*SINDR(I,NJ)) ADP0(I,I,NI)=ADP0(I,I,NI)+CB(NI,J)*(UDP0(I,J+I)*SINDR(I,NJ) I+UDP0(I,J+2)*SINDR(I,NJ+I)) ADP0(I,2,NI)=ADP0(I,2,NI)+CB(NI,J)*(UDP0(I,J+I)*SINDR(I,NJ+I) I-UDP0(I,J+2)*SINDR(I,NJ)) ADPP(I,I,NI)=ADPP(I,I,NI)+CB(NI,J)*(UDPP(I,J+I)*SINDR(I,NJ) I+UDPP(I,J+2)*SINDR(I,NJ+I)) ADPP(I,2,NI)=ADPP(I,2,NI)+CB(NI,J)*(UDPP(I,J+I)*SINDR(I,NJ+I) I-UDPP(I,J+2)*SINDR(I,NJ)) 180 CONTINUE AD0(I,3,NI)=UD0(I,NI)+UD0(I,NI+2) ADP0(I,3,NI)=UDP0(I,NI)+UDP0(I,NI+2) ADPP(I,3,NI)=UDPP(I,NI)+UDPP(I,NI+2) We have no interval subroutines to performe X--X TERM0=AD0(I,I,NI) AD0(I,I,NI)--AD0(2,I,NI) AD0(2,I,NI)=-TERM0 TERM0=ADP0(I,I,NI) ADP0(I,I,NI)=-ADP0(2,I,NI) ADP0(2,I,NI)=-TERM0 TERM0=ADPP(I,I,NI) ADPP(I,I,NI)=-ADPP(2,I,NI) ADPP(2,I,NI)=-TERM0 RETURN END
References [1] D. Michelson, Bunsen flames as steady solutions of the Kuramoto-Sivashinsky equation, SIAM J. of Math. Anal., submitted for publication (1989). [2] G. Sivashinsky, Turbulence in the motion of a free particle and de Broglie waves, Lett. Nuovo Cimento 27 (1980) 504. [3] G. Sivashinsky, Self-fluctuations in the motion of a free particle and the Schr6dinger equation, Lett. Nuovo Cimento 29 (1980) 135. [4] G. Sivashinsky, Elementary particle as a localized excitation of the action function, Nuovo Cimento 77A (1983) 21.