Dissections, transgressions, and perilous paths

Mathl. Comput. Modelling Vol. 28, No. 3, pp. 91-101, 1998 @ 1998 Eleevier Science Ltd. All righte reserved Printed in Great Britain 08957177/98 $19.00 + 0.00 PII: s0895-7177(98)00101-0

Pergamon

Dissections, Transgressions, and Perilous Paths R. ARE Department

of Chemical Engineering and Materials Science University of Minnesota Minneapolis, MN 55455, U.S.A. Dedicated

to Roe1 Westerterp*

(Received July 1997; accepted August 1997)

Abstract-A

model of the mace transfer from a bubble ricing through a liquid has been shown to have equations of the form s = A(z, y)/B(z, y), 6 = I’(%,y)/A(z, p). Thii ie one of a number of dime&one of the equation 2 = BP/AA, all of which have the earne solution patha. These equations will clearly have singular behavior on the curvea B = 0 and A = 0, but the boundary between physically real&tic and physically unreallltic values of z and y must be transgressed to reach this interesting behavior. It is shown that come of the dissections give families of singular trajectories which owe their perilous “atability” to the opposition of two infinite derivativce. @ 1998 Elaevier Science Ltd. All rights reserved.

Keywords--Phase

plane, Dlecontinuitiee, Singular traneformation, Ordinary differential equa-

tions.

TFtANSGFU3SSION [l]. A small gas bubble rises through a liquid to which the gas is transferred and in which it dissolves. The content of the bubble decreases as it ascends, but its size is governed by two opposing influences, The following

problem

arose in the discussion

of msss

transfer

in a chemical

reactor

the loss of content, which diminishes it, and the decreasing hydrostatic pressure, which allows it to expand. The total pressure in the bubble is the sum of the hydrostatic and the capillary pressures and the perfect gas law obtains. If z is the dimensionless volume, y the dimensionless content, and z the dimensionless depth of liquid that would exert the same hydrostatic pressure (i.e., if P is the pressure at the surface of the liquid, it (the surface) is at a putative depth of P/pg), the equations for the model sre [l]:

dx

z=

(Y-z?

dY

Y

,=,=;v

(y-q

where ty = 2 and p = 2/3. It can be shown that, even if there is a vacuum at the surface, 9 cannot be leas than 3zp, which means that the line y = zfi, along which 9 is infinite, is not in *Professor of Chemical Engineering at the University of Twente, The Netherlands. A master of the chemical reactor in ite many forma, he hae built both experimental and mathematical modele and has avoided the tranegreeeione here confeeeed.

91

92

R.AFLIS

the physically realiitic part of the z,Y-plane. This is fortunate, as calculations with (1) become impossible in the neighborhood of the curve y = L@. But the pair of equations that constitutes the model have a Me of their own, and we propose to transgressthe physical boundary and examine these equations in the whole quadrant. For the moment, we will consider system (1) in the positive quadrant with cu>l>P>O

(2)

and sketch some generalizationslater.

DISSECTION

OF EQUATIONS

Since the equations we are considering are autonomous, we can divide one by the other and consider the solutions of

The solution paths (a term we shall use to denote a curve that satisfies the equation without putting any arrowson it) of (3) are the same as those of(l), though it is not possible to put arrows indicating the movement of (z(z), y(z)) with increasingt. A solution of (1) can be represented by a curve invested with an arrow and this will be called a trajectory. Calculations with (3) are as doomed to failure as those with (l), for the singularitieswreak havoc with any numerical method in the neighborhoods of z = 0 and (y - z”) = 0. To get the solution paths, we must return to a pair of autonomous equations by dissecting (3). We do thii so that the derivativesof 2 and Y have no singularities,but recognizing that we have distorted the variable z. Thus, we write dx - =z(y-ZO),

$y(y-SD).

dt

Actually, to get any ease of computation, we need to take arc length as the independent variable and then g = X(Y -x9 E da

$=Y(Y-Z@) ds E

’

E2 = {z (y - z~)}’ + {y (y - z~)}~,

’

(4

and z can be recovered by quadrature. It is useful to state these equations more generally by writing W,Y)

= (Y -

B(x, Y) =

~9,

(Y

-

xp),

WAY) = Y,

A@, y) = z.

(5)

Then (1) is dz A -=-9 dz B and (4) is

da: AA -&=E,

dy

-&=E’

BI?

dy I’ -=-9 dz A da BA x=E’

(6) E2 = A2A2 + B212,

(7)

and (3) is J&= BI dx

a’

Equation (6) is said to be a dissection of (8), and (7) its “path-length dissection”.

(8)

P8thS

93

Table 1. Possible dissections of $$ = BI’/AA.

There are 16 differentways of dissectingequation (8). There are 32, if we allow a change of sign in each of the right-hand sides of (6), but thii merely reversesthe arrows on all the trajectories. They may be found in Table 1. The reasons for thii particular arrangementwill become clear later, but it is obvious that the singular curves, or loci along which a denominatorvanishes,will be of particular interest. In the special case r = y, A = x, for example, Type Bd has both axes and the curve B =I0 as singular curves. It will be seen that a,b,c,d, correspond, respectively,to neither the x-axis nor the ~-axis, since both axes are singularand A,B,C,D, connote, respectively,the singularityof neither A = 0 nor B = 0, A = 0, B = 0, both A = 0 and B = 0. The solution paths for all the variants are the same, but we shall see that the character of the phase plane ls quite different for differenttypes. There are two paths from the origin to the point (l,l), i.e., from the point where I’ = A = 0 to that at which A = B = 0, They enclose a region, S, filled with closed paths that all come into the origin from above as y = z@ and are infinitely flat along the sx~& i.e., all their derivativesgo to zero as the origin is approached (see Figure 1). The two arms of the bounclmg paths continue to infinity and divide the plane into four regions-N, S, E, W-within any one of which the character of the paths is the same. The general aspect of the phase plane is that of a saddle point (which, indeed, (1,l) sometimes is) with the southern and westerly separatricestied together at the origin (see Figure 1). By linearizationabout the point (l,l), it can be shown that the slopes of the bounding curves as they pass through (1,l) are ~1and V, where

(9)

0 so that we see that the potentially singular NowcrandparetheslopesofA=OandB= curves lie in the quadrantsN and S. We divide N and S into three parts: SW is the region between B = 0 and the arm of the separatrix that runs between S and W; similarly, SE is the region

R. ARIS

94

Figure 1. The first quadrant of the phase plane and its regions.

between the S/E branch of the separatrix; the central part of S will be S’; i.e., S=S’+SE=SW. Similarly,N=N’+NE+NW.

BEHAVIOR

NEAR

THE ORIGIN

Consider a path coming into the origin tangentially to the y-axis. If we assume that the asymptotic behavior of such a path is of the form y = Cxr, we can substitute this in (3) and keep the dominant terms, giving Cyxr-’ = Cxr-1 { (1 - zP/Cxr) / (1 - x”/Cz’)}. If +y= 0, we shall be able to match the constant on the lowest power of z and have remainderterms which go to zero like x+fl/C. In this way, we find that the common asymptotic behavior of all solutions that come to the origin tangential to the y-axis is Xfl

y=(1-

(10)

It is not possible to extend thii expansion in a general formula without knowing the relationship of o and /3. In particular, if the greatest integer in a/P is n, the first n - 1 terms will be powers of xp. If, on the other hand, x and y both go to infinity, keeping y > x4 (as in region W), we can only match terms if +y= cr and obtain y=$$).

(11)

When we consider the paths that come in tangentiallyto the x-axis, we recognize that y < x0, and we must write the equation in the form $ = (y/~P-~+l) { (1 - y/xfl (1 - y/x”)}, and then the trial function that satisfiesthe lowest order terms is

with C arbitrary. Since p < o, this shows that the asymptotic form and all its derivatives go to zero near the origin. But the same form applies as x goes to infinity and y remains bounded,

Path8

95

for, under these conditions also, y/z” and y/zP tend to zero. The paths in region E can be parametrized by the constant C; the larger C is the more closely will the path follow the eastern border of the region S, which is one arm of the separatrix of the saddle point, and switch to the eastern border of N. The middle part of the path must, of course, be calculated from the differential equation; it is only the two ends that behave in this Arrhenius fashion of (12).

TRAJECTORIES Trajectories are paths with direction given by the grouping of the factors in the dissection under consideration. Fist, let us note that the axes are trajectories; the z-c-axis,I = y = 0, is a path for all Ma and all MC. (“All Ma” means Types Au, Ba, Ca, Da and similarly, Am is the set Aa, Ab, AC, Ad.) Since A and B are both negative when y = 0, the sign of 6

will depend

on how many of the two factors A and B appear in the right-hand side of the equation; thus Types Am and Dm have trajectories that go to the left along the x-axis and those of Types Bm and Cm go to the right. Because all factors are positive, the y-axis is always traversed upwards. But the path can be horizontal, not only by the vanishing of 2 = 0, but also by 2 becoming infinite since then it does not matter what $ is. Thii happens in the Types Mb and Md. Furthermore, the trajectory y = 0 is traversed with infinite speed in Types Mb and Md. Likewise, the yaxis is a streak, as we might call a trajectory of infinite speed. Actually, the v-axis is in the physical region of the phase plane, and, recognizing that the equations for the rising bubble do have negative signs, a starting point on the y-axis instantaneously collapses into the origin. This makes good sense, for if a bubble has no volume, it can have no content. Similar considerations apply to the loci A = 0 and B = 0, for the trajectories run (respectively) vertically and horizontally across them. But, when A and B appear in the denominators, their vanishing makes the trajectories run horizontally and vertically by overpowering the other derivative. Thus, Types Am will have identical open phase planes (z, y > 0), since only I’ and A appear in denominators. They take the form of the usual trajectories about a saddle point knotted at the origin, which is clearly a more complex critical point than is usually encountered. The family of homoclinic trajectories that fill the region S are traversed clockwise, all sprouting from y=xp/(l -p) and ret urning infinitely flatly, so to speak, to the origin.

FINITE

TIME

TRAVERSES

In contrast to the common exponential approach to a singular point, which takes an infinite time, some of these trajectories escape from the origin and/or return to it in finite time. Consider the asymptotic form of the trajectory as it approaches the origin tangentially to the y-axis,

and let Z,(Y) be the “time” to traverse the small interval from (X, Y), a point on this asymptotic curve, and the origin. Keeping only the dominant terms in each formula, we have on thii pathway: A=y-x”=

[A] {1+[(4-~n+-1,!;1-8)]Xu-~},

B=y-zP=

[&]

r=y=

[&)I

{I+ {1+

[H]z@},

[H]z@},

A = x. The inversion of the function I’ gives first terms B/P = I’ = y = zo/(l - p), A = x. If the dissection is of type Ad, g = x’(x) = A/I’ = 1, so that Z,(X) = X. We could alternatively

R. AFLIS

96

Write Y’(Y)= B/A = By/z = pza-‘/(l - p) = p(1 - p)-1/B@-1/8, and the integral of l/y’ with respect to y is [(l - p)yl’Ip = X. Of course, this is only asymptotically valid and it would be a mistake to try to find a trajectory by integrating away from the origin, for the random rounding error would cast one onto a random solution. Rather one should integrate from a point on y = 28, say, where 9 is at a maximum, simultaneouslyadding on the approximation Z,(X) to give an estimate of the total time. Letting X tend to zero will show convergence to the total time from the origin to B = 0. A similar treatment of each dissection is possible, for each type has an asymptotic value of the form x’ = QxQ, and this will give a time between the origin and a point (X, Y) on an asymptotic path of 2, = [l/Q(l -q)]X l-q. Clearly,q must be less than 1 if this is to be finite. Table 2 shows the 16 diitions and gives the formulae or indicates the cases where the time is ix&rite; Y is used as an abbreviation for Xp/(l - p). Table 2. Asymptotic timea to the origin.

zu

33 Q

P+l

00

1

co

.1-P

P

(13) al

Bb

IBcI

I I

Bd

Ca

1

I 0

I I

0

03)

I

1

0

00

-l(a

-PI

(13) rP-a

a-P+1

I

(a - PI (13) 0

Cd

a-P

Da

1-P

P

-P

1

-U+P)

I

-p

I

-O+P)

I

PX

I

(13)

1 DC 1 Dd

I

X0 -iT

.I+8

(1+ PI (13)

To find out whether the flat trajectories converge on the origin in finite time, we observe that 3 <<
Paths

97

i.e., has a factor of l/I’,

z&.(X)= ; 0

w - 1) 7p-1r(p,u), P'L, 7

u/-7 7'

where IQ, V), the integral from U to inf%ty of ~p-~e-”

7= b--P)?

(13)

with respect to U, is the incomplete

gamma function. The singular line B = 0 is also reached from a nearby point in fmite time. To see this, let (X, Xs) be a point on the singular curve B = 0, and let a displaced point of the same ordinate be (X + B,Xp). We choose to put it on the same level Y = Xfl because we know that the path goes horizontally through the singular curve. Now, in seeking the asymptotic form, we may give A, r,and A their values at (X, Y), namely (X” - Xp), X0, and X, respectively, and need only remark that, if only the dominant term is kept, B can be approximated by pXa-lt, where < = z - X. Now this singularity belongs to the dissections of types Cm and Dm and the time will be finite since the vanishing factor appears in the numerator of the integral for z. Let ZB(B) be the “time” to the curve B = 0 when it is singular, then in Type Ca, for instance, 2$(B) = ,8Xp-2 B/Xfl - Xa). A similar treatment can be applied to the approach to A = 0 when it is singular, as it is for Types Bm and Dm. Expressions for ZB(B) and 2,4(H) are given in Table 3. We note that the asymptotic times are proportional to squares of the displacements, and that the vertical and horizontal distances between the two singular curves play a role in Types

Bm and Cm. Table 3. Asymptotic times of approach to singular curves from (X,X (Y + H, Y).

23B (2) -2 =.

2ZA(ff)

Type

+ 2) and

- H2 y-1 (Y - Y@b)

Ba Bb

(Y -:el.) y-1+1/0

Bc

(Y - Yfl/O) yl/a

Bd

(Y - Yfl/(l) pxo--2 (X0 - xq pXw-2

Ca Cb

(x= - xq pX+1 (X” - xq /3X284

cc Cd

(X4 - xq

I

Da

I

y-1

1

Db

1

1

1

DC

1

Y-l+l/‘=

I

I

Dd

I

pX+2

1

pX2+2

1

pxB-1

I

yVa

THE

/3X28-1

FLOWS

We can now put arrows on the paths to make them into trajectories, and so describe the flows corresponding to each possible dissection. We have seen how the axes can be either streaks or

R. ARIS

98

trajectories. Thus, the y-axis is always an upward going trajectory and is a streak for Types MC and Md. The z-axis is a streak for Types Mb and Md, and points to increasing values in x for Types Bm and Cm, but to the origin for Types Am and Dm. There are four distinct flows strictly within the quadrant, namely those of Types Am, Bm, Cm, and Dm, which have none, one, one, and two of the singular curves, respectively. It is easy to see that the trajectories point to increasing x and y above the curves A = 0 and B = 0, i.e., in the regions SW, W, and NW. The paths have positive slope in the complementary regions SE, E, and NE, but go to the right for Types Bm and Cm and to the left for Types Am and Dm. In Figure 2, we show the four flows with the singular curves as broken lines. Since both of the factors A and B are in numerators in the dissected equations, there are no discontinuities on either of the curves A = 0 or B = 0 and these are merely the isoclines of vertical and horizontal flow. Thus, the region S is filled with closed solution paths traversed clockwise, and, in the case of Type Ad, in a finite time.

1.5

1

0.5

Figure 2. ‘Ikajectories in the

SOLUTIONS

quadrant and

on the axea for the 16 dissections.

OF PERILOUS

STABILITY

For Types Bm, the derivatives change sign discontinuously across the curve y = xQ, close to which they are arbitrarily large. Thus, in the region S, a trajectory originating from the origin and starting vertically upwards goes through a maximum, and just where it is moving vertically downwards, it meets another coming in the opposite direction, vertically upwards. Thii takes

Paths

99

place on the curve A = 0, as we can easily see from the equations for type Ba,

Now l/A is not defined on A = 0, but we will let it be defined there as

0; Then 2

-A, (X, XQ) A (X, xa)

(X,XU) =

A, (X, Xa) B (x,X”)

I? (x, xa) ’

(14

will equal -A, (X,X”) /A, (X,X% and the solution will turn a sharp corner and, is positive, will proceed up the curve A = 0. Notice that it goes straight through the

because 2

saddle point (l,l),

and indeed, accelerates, since g = z.

Thii path, at least until it reaches the saddle at (l,l),

has a certain stability, which might be

called perilous. It is as though one were to walk the ridge of a mountain, with steep slopes of scree going down into two valleys on either side. This is not a happy situation, for one false step to one side will land you in one or the other of the valleys. “Whan eek Zephyrus with hi sweete breeth,” helped, perhaps, by an inverted Janus, comes to your rescue and makes the ridge the stagnation plane of two violent winds coming from the two sides. Thii situation is still far from happy, even if the winds are perfectly steady, for now, if you step off to one side, you will be picked up by the wind and flung violently back and forth from one side to the other, like responses in the Puritan conception of the Episcopalian liturgy. The situation is even worse when you get beyond the saddle, for Janus is really Janus, the winds are reversed, and an exit of gadarene violence is ensured under the double jeopardy of storm and steep place. We observe that thii perilous solution has been purchased by special definitions of A and B along the curves where they are zero, namely, A(X, Xcr) =

B (X,X13) =

-A,

(X, Xa) B (X,X”) r (x, x”) ’ A, (X, Xe) A (X, Xu)

-B,

(X, X”) A (X,X0)

(15a)

A (x,x13)

B,(X,Xfl)r(x,xfl)

*

(15b)

These definitions do not depend on the particular dissection that is under discussion, though the existence of perilous solutions does. The experience of being thrown backwards and forwards across these ridges can be had computationally by trying to solve any dissection except those of Types Am, and of these, only Aa is really suitable. The necessity of taking a finite step in any computation makes it impossible to follow a perilous path.

A SINGULAR

TRANSFORMATION

The behavior of thii system is illuminated by a singular transformation

u=--,

X0

v=-

Y

UT V

0

(16)

Y’

with inverses x=

XP

)

Y=

up y

(> ‘vo

9

(17)

where y=&.

(18)

100

R. ARIS T

R

Figure 3. Comparison of phase planes under the singular transformation.

The transformed solutions are clearly adapted to the problem. The line v = 1 - p, for instance, is a map of the asymptotic path y = zfl/(l - 0) to the origin in the z,y-plane, while approaching u = (CY- 1)/a, v = 0 is to go to infhrity along y = (o - l)sQ/a. The differential equation for u and v can be obtained most easily by logarithmic differentiation of (16) and substitution in (3), and is dv x

{vpu-v+l--1) = (u[ou

-v-(a-l)]}’

(12)

Again, the path-length dissection gives a pair of equations that can be integrated without numerical difhculty, and the two parts of Figure 3 show the two planes with corresponding points and asymptotes marked with the same letter. Notice that the asymptote of all solutions that come in to the origin from above has become a stable node, 0. On the other hand, the very flat approach to, or departure from the origin (UMO in the z,y-plane) that we noticed in equation 12 is represented by the asymptote SRO’ in the u,v-plane, since y goes to zero more rapidly than any

Paths

101

power of x. The closed paths in the x,y-plane, such as OPNQO, are mapped into curves such as O’PNQO in the u,v-plane. Trajectories like TVO, TWO, and TX0 become heteroclinic paths between the unstable node T (which corresponds to the asymptotic behavior when x and y are large) and the stable node 0. There is a unique path TWO in the u,v-plane, which arrives at 0 tangentially to the eigenvectorof the lessereigenvalueof the Jacobian at 0. Thii separatespaths in the u,v-plane which come in from above, from those which come in from below. Near 0 can be representedby the line v = 1 - p - 021,u = p(1 - p)/(cy - p). This shows that the asymptotic form of TWO near 0 is y = xp/(l - ,@. Similarly,the eigenvector at u = (o - 1)/a, v = 0 gives the asymptotic behavior at the T-infinity in the x,y-plane, namely y = CYX~/(CX - 1) + pxfl where p = (o - l)/ {((Y - 1)2 - (1 - /?)}. The first terms of these asymptotic formulae agree with equations (lO),(ll). The very flat approach along the x-axis of trajectories in the x,y-plane is mirrored in the asymptotics of the paths going to the O/-infinity in the u,v-plane. On the path OR in the x,yplane, y is very much less than xQ, which, in its turn, is much less than xp; hence, u and v are both large and v >> u. Thus, we look for the solution of equation (19), which gives the path SRO’. But (cr - 1) and (1 - 0) can be neglected as small in comparison with u and v and then equation (19) can be rearrangedas

(k)(g)-($)+;-(;)(g)=o. Thii can be integrated immediatelyto give

;+1n

($ >=

constant,

(20)

which is transformed back into equation (12) by the use of equations (18),(19).

DISCUSSION There does not seem to be much literature on solutions of discontinuous derivative, though there is a vast amount on the two-dimensionalphase plane. Lefschetz defines the structure in the region S of Figure 1 as an elliptical sector [2, Figure 12, p. 2091,while Asnov and Arnold [3] classify it as the phase portrait of the quotient system of a nondegenerateequation in class Wit’. Lefschetz [2] treats the asymptotic behaviour by means of the stereographic projection onto a plane. Some years ago, Hardy [4] established the ultimate monotonicity of solutions of such equations as ours and went on to obtain rate of growth estimates. Since it is surface tension that gives rise to the denominator B = (y - x0) in the mass transfer problem that generated thii transgression,I had hoped at first that a discontinuous solution might have been the implosion of the bubble as the radius became so small that the term a/r took over and dwarfed all others. However, this was not to be, and it is easy to see that the perilous solution cannot be physically realistic.

REFERENCES 1. Ft. Arii, A note on masa transfer froma smallascending bubble,Chem. Enging. Sci. 52,443Q4446, (1997). 2. S. Lefschetz, Diffferential Equations: Geometric Theory, Interscience, New York, (1957). 3. V.I. Arnold and Yu.S. Il’yeshenko, Ordinary differential equations, In Art. I of Encyclopedia of MathematicaI Scienceu, (Edited by D.V. Asnov and V.I. Arnold), Springer-Verlag, Berlin, (1985). 4. G.H. Hardy, Some results concerning the behaviour at infinity of a real and continuous solution of an algebraic differential equation of the first order, In Proc. London Math. Sot. Ser. 3, pp. 451-468, (1910).