An algorithm for autonomously plotting solution sets in the presence of turning points

An algorithm for autonomously plotting solution sets in the presence of turning p oin ts Stev en P ollack a, ∗ , Daniel Badali b, ∗∗ , Jonathan P ollac k c a Dep artment of Mathematics, McGil l University, Montr e al, Queb e c, Canada, H3A 2T5 b Dep artment of Chemic al and Physic al Scienc e, University of T or onto at Mississauga, Mississauga, Ontario, Canada, L5L 1C6 c Dep artment of Me chanic al Engine ering, McGil l University, Montr e al, Queb e c, Canada, H3A 2T5 Abstract Plotting solution sets for particular equations may b e complicated by the exis- tence of turning p oin ts. Here w e describ e an algorithm which not only o v ercomes suc h problematic p oints, but do es so in the most general of settings. Applica- tions of the algorithm are highlighted through tw o examples: the first provides v erification, while the second demonstrates a non-trivial application. The latter is follow ed b y a thorough run-time analysis. While b oth examples deal with biv ariate equations, it is discussed ho w the algorithm may b e generalized for space curv es in R 3 . Keywor ds: turning point, implicit function, cusp, bifurcation curve 1. Introduction In this pap er we consider curv es determined b y equations of the form f ( x, y ) = 0 (1) where f : I → R , I ⊂ R × R is a product of op en in terv als. Equations such as (1) are often referred to as implicit equations since in general there do es not exist an explicit, unique function g such that y = g ( x ). A canonical example is that of the unit circle, with f ( x, y ) = x 2 + y 2 − 1 = 0, which cannot b e rearranged to isolate y as a function of x . The following discussion also applies to equations of the form f ( x ; α ) = 0 where α is a bifurcation parameter. The purpose of this pap er is to address the problem of plotting solution curv es to (1) with turning p oints (a precise definition will be presented shortly). ∗ Corresponding author ∗∗ Principal corresp onding author Email addr esses: steven.pollack@mail.mcgill.ca (Stev en P ollack ), daniel.badali@utoronto.ca (Daniel Badali), jon.pollack@mail.mcgill.ca (Jonathan Pollac k) Pr eprint submitte d to Elsevier Novemb er 20, 2018 While a metho d for dealing with this problem has already b een developed by Keller [7, 8] (via a pseudo-arc-length parametrization), it relies on the turning p oin ts to be anything but cusps, and th us cannot be used on curv es without some understanding of their profile. Con versely , the algorithm that we will presen t not only requires no prior knowledge of the solution curv e, but also approac hes the problem in what w e consider to b e a more intuitiv e manner. It should b e said, ho wev er, that although we ackno wledge that there are other methods to deal with this problem, we will make no attempt to compare them. Essen tially , this pap er is to b e self-contained, and therefore our only concern is the explanation of our prop osed algorithm, and its o wn b enefits, not its relativ e benefits. 2. Background Before a discussion of the algorithm may b egin, we must first define the term turning p oint . Definition 1. Let f : R 2 → R , and S =  ( x, y ) ∈ R 2 : f ( x, y ) = 0  . A p oint ( x ∗ , y ∗ ) ∈ S is a turning p oint of f ( x, y ) = 0 if ( x ∗ , y ∗ ) ∈ S , and there exists δ > 0 such that at least one of the follo wing statemen ts is true: T yp e 1:  ( x, y ) ∈ R 2 : 0 < x − x ∗ < δ, 0 < | y − y ∗ | < δ  ∩ S = ∅ . T yp e 2:  ( x, y ) ∈ R 2 : 0 < | x − x ∗ | < δ , 0 < y − y ∗ < δ  ∩ S = ∅ . T yp e 3:  ( x, y ) ∈ R 2 : 0 < x ∗ − x < δ, 0 < | y − y ∗ | < δ  ∩ S = ∅ . T yp e 4:  ( x, y ) ∈ R 2 : 0 < | x − x ∗ | < δ , 0 < y ∗ − y < δ  ∩ S = ∅ . Figure 1: Illustration of the different types of turning p oints Definition 1 may be envisioned in the follo wing wa y: if one were to restrict themselv es to a small enough neighborho o d ab out ( x ∗ , y ∗ ), one of the op en 2 ∂R + R S k step 1. Type 1 turning point reference point minimal point step 2. v y v x v step 3. Figure 2: Demonstration of algorithm on a Type 1 turning p oint. “half-balls” centered at ( x ∗ , y ∗ ) contained in this neighborho o d would hav e no in tersection with S (see Figure 1). The remainder of this paper will deal with the discrete set: S k = n  ( x j , y j )  N j =1 ⊂ R 2 : f k ( x j , y j ) = 0 , 1 ≤ j ≤ N o (2) where S k is the appro ximation of S , and f k is the discrete appro ximation f . F or a detailed discussion of approximating turning points in discrete spaces see the review b y Cliffe, Sp ence, and T av ener[9]. 3. Algorithm Supp osing that ( x j , y j ) ∈ S k has b een iden tified as a turning point of any t yp e, let R be the closed disk of radius r centered at ( x j , y j ). The algorithm is then performed in the following steps: Step 1. Uniformly scan the boundary of R , ∂ R , for solutions to f ( x, y ) = 0. Denote the set of solutions as R . Step 2. (a) If R = ∅ , we ma y stop here: no other points in R 2 satisfy f ( x, y ) = 0. (b) If R 6 = ∅ , select ( x i , y i ) ∈ S k , i < j , to b e a “reference p oint” and select the minimal p oint ( s 0 , t 0 ) ∈ R suc h that φ ( s 0 , t 0 ) ≤ φ ( s, t ) for all ( s, t ) ∈ R , where φ ( s, t ) = 1 p ( x i − s ) 2 + ( y i − t ) 2 (3) Step 3. (a) Create the v ector v = ( s 0 − x j , t 0 − y j ), 3 (b) Determine the largest axial comp onent of v , d = max  | s 0 − x j | , | t 0 − y j |  . (c) Use d to determine the new direction of iteration. T o clear up any p ossible ambiguit y , w e will elaborate on how one migh t implemen t the abov e strategy on a T yp e 1 turning p oint. Step 1. Since ( x j , y j ) is a T yp e 1 turning p oint, we kno w that there exists an op en ball of radius δ cen tered at ( x j , y j ) whose eastern hemisphere is disjoint from S k . Therefore, if any solution to the equation f ( x, y ) = 0 exists, it must exist in the w estern hemisphere. Hence, we describ e ∂ R + , our reduced search path, b y ∂ R + =  ( x, y ) ∈ R 2 : ( x − x j ) 2 + ( y − y j ) 2 = δ 2 , x < x j  whic h is easily seen to be parametrized b y g ( θ ) = ( δ cos( θ ) + x j , δ sin( θ ) + y j ) , θ ∈ [ π / 2 , 3 π / 2] Th us, if n p oints are used to uniformly sample ∂ R + , the mesh would lo ok lik e g i =  δ cos  π 2 n ( n + 2 i )  + x j , δ sin  π 2 n ( n + 2 i )  + y j  for i = 0 , 1 , . . . , n − 1. The p oints g i are then chec ked to see if f ( g i ) = 0. If g i is solution to f ( x, y ) = 0, it should b e recorded in the set R . It should be noted that the metho d emplo y ed to c hec k whether g i is a solution to f ( x, y ) = 0 is irrelev ant. Step 2. W e hav e found that it is b est to pick an ( x i , y i ) that is contained in R . In this wa y , when we minimize the cost function φ , see (3), we are searching for the the point in R which is furthest from our the reference p oin t. If w e assume ( x j , y j ) is not the final p oin t of S k , then the furthest p oint from the reference p oin t is not exp ected to b e an y of the previous p oints in S k . Picking a reference p oin t that is “to o far” from ( x j , y j ) can p otentially lead to choosing ( s 0 , t 0 ) ∈ R suc h that ( s 0 , t 0 ) is a previous solution. When the algorithm is implemen ted inside an iterated plotting program (one w ould do this for autonomous plotting), this has the consequence of leading the plotting program to retrace its old steps. F rom here it is only a matter of ev aluating φ ( s, t ) for all ( s, t ) ∈ R , and finding the smallest v alue. The co ordinate which minimizes φ shall b e called ( s 0 , t 0 ). While there is no garen tee that there exists a single co ordinate, suc h that φ is minimized uniquely; in practice, we seldom arriv e at t wo, or more p oin ts, that equally minimize φ . Step 3. Mak e the direction v ector, v , from ( x j , y j ) to ( s 0 , t 0 ). Then identify whic h comp onen t of v is larger in magnitude. If the y -comp onent is larger, b egin iterating along the y -axis (the direction is determined by the sign of the y -comp onent). Iteration should thence begin at ( s 0 , t 0 ). This is to comp ensate for the fact that it ma y happ en that the turning p oint was approac hed via iteration along 4 Figure 3: Astroid plot generated using the turning point algorithm. Inset: parametrized curv e (solid), approximated curve (circles) the p ositive x -direction, yet the algorithm has concluded that the new direction is along the negative x -direction. If w e do not b egin at ( s 0 , t 0 ), the program will just bac k-trac k along already found p oints of S k . It should b e noted that while the algorithm has been presented in the t wo- dimensional case, it is easily generalized to apply to equations on R n . F or instance, if w e were to consider the case f ( x, y , z ) = 0, then a simple mo dification of Definition 1 including half-spheres, instead of half-balls, w ould allo w for this algorithm to work with space curv es. 4. V erification T o demonstrate the effectiveness of the turning p oint algorithm, we chose to implemen t it on the astroid, a w ell kno wn implicit equation giv en b y f ( x, y ) = x 2 / 3 + y 2 / 3 − 1 = 0. The resulting curve can b e found in Figure 3. The accuracy of the curve generated with our algorithm w as measured in p ercen t error (4), and dep ended on several differen t parameters: r , the radius of ∂ R + , k how many “steps” back the reference p oint was from the turning p oint, and n the grid size of ∂ R + . Each parameter was v aried until notable v alues w ere found. P ercentError( x j , y j ) = q x 2 j + y 2 j − p x 2 + y 2 p x 2 + y 2 × 100% (4) 5 It was found that a radius r ∈ [0 . 0001 δ, 100 δ ] (where δ is the magnitude of the step-size ∆ x or ∆ y ), 1 ≤ k ≤ 10 and 4 ≤ n ≤ 10 yielded a maxim um percent error  0 . 1%. Ho wev er, for n < 4 the algorithm failed to correctly navigate the turning point, and the plot ceased to contin ue. 5. Example: Lubrication Mo del Consider the problem of determining the thic kness of a thin film of lubri- can t inside a rotating cylinder. Using the lubrication approximation to the Na vier-Stokes equations dev elop ed in [10, 11], one ma y describ e the steady- state thic kness, h ( θ ), by  3 ( h 0 + h 000 ) − 1 3 cos( θ ) = Q h 3 − 1 h 2 (5) where Q is the non-dimensional flux and  is dep endent on the surface tension and rate of rotation. Coupling (5) with M = Z 2 π 0 h ( θ ) dθ (6) it is not unreasonable to ask ho w the curv e of points satisfying (5) and (6) b eha ves for a particular  . Because this curve ma y be c haracterized b y a solution set to some equation f  ( Q, M ) = 0, w e will let f  denote the curve. Since our in tention is to demonstrate the robustness of the algorithm, and not the c hallenge b ehind solving f  ( Q, M ) = 0 when either Q or M is fixed, we refer the curious reader to the works of Ashmore et al. [14] and Benilov et al. [12] for detailed exp osition of that task. Briefly , we discretized h on a uniform grid, and use F ourier sp ectral dif- feren tiation matrices whic h account for the 2 π -p erio dicity of h . A MA TLAB program w as written to accept very small v alues of  and p erform an iterated Newton-Raphson procedure to generate f  giv en a particular parameter and direction. Once a turning point was identified, the algorithm w as p erformed. It w as found that for  ≤ 10 − 3 complicated loop b ehavior similar to that exhibited in Figure 4 app eared. With the algorithm deplo yed to automate plot- ting, eac h successiv e turning p oint w as o v er come with no difficulty . In terested readers may also refer to [12, 13, 14] for recent literature on this topic, as w ell as [15, 16, 17]. 6. Run-Time Analysis Due to the v arious settings in which the algorithm ma y b e applied, a general run-time analysis is impossible. Instead, we will analyze the case when f is not explicitly given, but inferred from differential or in tegral equations. In particu- lar, we will perform a run-time analysis on the Lubrication Mo del example. As demonstrated in Section 5, it is clear that these settings are very useful in the study of certain ordinary and partial differential equations. 6 Figure 4: A p ortion of the ( Q, M )-bifurcation diagram generated for  = 10 − 4 . Mo dified from [18]. In this particular instance, w e had tw o particular equations (5) and (6) whic h allo wed us to solve for Q and M given a function h : [0 , 2 π ] → R + . By appro xi- mating h with h m (an m × 1 column vector), and taking the appropriate T aylor expansion of elements of (5), w e w ere able to p erform parameter con tinuation along the M -axis with a m × m matrix appro ximation of the op erator f . With similar linearization tec hniques w e were able to p erform parameter con tinuation along the Q -axis with a ( m + 1) × ( m + 1) matrix approximation of f . (F or more details ab out these tec hniques see Benilov [12].) In this setting, p erforming the Newton-Raphson metho d prescrib ed by Ash- more [14] and/or Benilov [12], and thus completing Step 1 in the algorithm, is b ottle-nec ked by the sp eed of matrix multiplication. No w, an algorithm for matrix m ultiplication of m × m matrices has been found with run-time O ( m 2 . 376 ) (i.e. the Copp ersmith-Winograd algorithm[19]). Ho wev er, this algorithm is particularly hard to implemen t, and is usually used for theoretic bounds abov e all else. Th us, in the spirit of creating a practical run-time argumen t, w e assume that Strassen’s algorithm is used, and therefore the calculations in Step 1 can b e carried out in O ( m 2 . 807 ) [20]. Also, because Newton’s metho d is usually p erformed in a for-lo op for a fixed num b er of itera- tions, ` say , Step 1 ma y b e preformed in O ( n ` m 2 . 807 ) = O ( n m 2 . 807 ), where n is the num b er of sample p oints on ∂ R + . No w, b ecause calculating φ is reasonably assumed to take constant time, it is clear that the most cum b ersome part of Step 2 is determining the point whic h minimizes φ (which is a sorting problem). Thus, if w e assume R has n p oints in it (the maximal amount), then Step 2 ma y be performed in Θ( n log n ). This follo ws from the fact that [20] sho ws that that an y comparison sort algorithm 7 requires Ω( n log n ) in the worst case. Hence, taking any v ersion of heap or merge-sort to p erform Step 2 will require Θ( n log n ). While there are a m yriad of wa ys one may p erform Step 3, any comp etent implemen tation should run in constan t time since this step has only manipulates t wo p oints in R 2 . Hence, the algorithm can b e found to run in O ( nm 2 . 807 + n log n ). How ever, as we sa w in the astroid example, n need not exceed 10, or so, p oints. Th us, if we approximate h with any reasonable amount of resolution (i.e. m = 2 7 or m = 2 8 ), then n  m , and we ma y consider it something of a constan t. Hence, the algorithm can b e expected to run in O ( m 2 . 807 ). 7. Conclusion It is easy to see that turning points make automated plotting a challenge, if not imp ossible. T o o vercome this issue w e ha v e outlined an algorithm that determines the direction iteration should contin ue in after reaching a turning p oin t. As demonstrated with the astroid and thin-film differential equation examples, the algorithm ov ercomes turning p oin ts accurately , and effectively , while asking little in terms of computational p ow er. It should be said that implicit functions and dynamical systems far from exhaust the areas where this algorithm may be helpful. With the amount of research being done to day throughout the v arious sciences (pure and applied alike), it is clear that this algorithm can b e applied an ywhere data visualization is p ossible. Ac kno wledgments : This w ork was partially carried out during the Fields- MIT ACS Undergraduate Summer Researc h Program in 2010. The authors thank M. Chuguno v a for helpful discussions during this program. References [1] G. Moore and A. Spence, The calculation of turning points of nonlinear equations, SIAM J. Numer. A nal. 17 (1980), pp. 567-576 [2] G. Mo ore and A. Sp ence, The conv ergence of op erator approximations at turning points, IMA J. Numer. Anal. 1 (1981), pp. 23-38 [3] A. Sp ence and B. W erner, Non-simple turning p oints and cusps, IMA J. Numer. Anal. 2 (1982), pp. 413-427 [4] R. Seydel, Pr actic al Bifur c ation and Stability Analysis: F r om Equilibrium to Chaos , Springer-V erlag, New Y ork, 1994 [5] E.J. Do edel, H.B. Keller, and J.P . Kernevez, Numerical analysis and control of bifurcation problems (I): Bifurcation in finite dimensions, Internat. J. Bifur. Chaos Appl. Sci. Engr g. 1 (1991), pp. 493-520 [6] D.F. Da videnk o, On a new metho d of numerical solution of systems of nonlinear equations, MR 14 (1953), pp. 906 8 [7] H.B. Keller, Numerical solution of bifurcation and nonlinear eigenv alue problems, in Applic ations of Bifur c ation The ory , Academic Press, New Y ork, 1977 [8] H.B. Keller, L e ctur es on Numeric al Metho ds in Bifur c ation Pr oblems , Springer-V erlag, New Y ork, 1987 [9] K.A. Cliffe, A. Spence, and S.J. T av ener, The n umerical analysis of bifur- cation problems with application to fluid mechanics, A cta Numer. (2008), pp. 39-131 [10] H.K. Moffatt, Beha vior of a viscous film on the outer surface of a rotating cylinder, J. de M´ ec. , 16 (1977), pp. 651-673 [11] R.E. Johnson, Coating flow stabilit y in rotating molding, Engine ering Scienc e, Fluid Dynamics: A Symp osium to Honor T.Y. Wu (e d. G. T. Y ates) , W orld Scientific, 1990 [12] E.S. Benilo v, M.S. Benilov, and N. Koptev a, Steady rimming flo ws with surface tension, J. Fluid Me ch. 597 (2008), pp. 91-118 [13] E.S. Benilo v, M.S. Benilo v, and S.B.G. O’Brian, Existence and stability of regularized shock solutions, with applications to rimming flo ws, J. Eng. Math. 63 (2009), pp. 197-213 [14] J. Ashmore, A.E. Hosoi, and H.A. Stone, The effect of surface tension on rimming flo ws in a partially filled rotating cylinder, J. Fluid Me ch. 479 (2003), pp. 65-98 [15] R.E. Johnson, Steady-state coating flo ws inside a rotating horizontal cylin- der, J. Fluid Me ch. 190 (1988), pp. 321-342 [16] S.T. Thoro ddsen and L. Mahadev an, Exp erimental study of coating flows in a partially-filled horizontally rotating cylinder, Exp. Fluids 23 (1997), pp. 1-13 [17] P .L. Ev ans, L.W. Sch wartz, and R.V. Roy , Three-dimensional solutions for coating flow on a rotating horizontal cylinder: Theory and exp erimen t, Phys. Fluids 17 (2005), pp. 172102:1-20 [18] D. Badali, M. Ch uguno v a, D.E. Pelino vsky , and S. P ollack, Asymptotic b eha vior of regularized shock solutions in coating flows, submitted to Phys. Fluids (Jan uary 2011) [19] D. Copp ersmith and S. Winograd, Matrix m ultiplication via arithmetic progressions, J. Symb olic Comput. 9 (1990), pp. 251-280 [20] T.H. Cormen, C.E. Leiserson, R.L. Rivest, and C. Stein, Intr o duction to A lgorithms, Se c ond Edition , MIT Press and McGraw-Hill, Cambridge, 2001 9