Recursive Journey to Three Bodies


[Taylor solver, Taylor Method, ODE, Automatic Differentiation]


The “Recursion Excursion” started in Delphi Informant 8/2000, now continues and will bring us as far as to a full featured Taylor Solver – actually the Taylor Center with powerful GUI and graphics. On the way, we will learn more, where the recursion is good and were it isn’t, and finally we will have an opportunity to play with such an illustrative sample of differential equations as the Problem of Three Bodies.


1. Introduction.  


We have considered the examples where recursive procedures delivered elegant and powerful solution: it was parsing of an arithmetic expression into a postfix sequence, and then its evaluation. Starting from this very point, now we can reach much farther and to implement an essentially more complex application – the modern Taylor method, an interesting and powerful numeric method for integrating the Initial Value Problems for the Ordinary Differential Equations (ODEs).     


In the Recursion Excursion we considered the algorithms evaluating the system of recurrent formulas, the so called Linear Elaborated Declarations:


u0 = f0({numeric values only});

u1 = f1(u0, {numeric values});                                            (1)

.  . . .

un = fn(u0, u1,…, un-1, {numeric values});


where each variable first appears in the left hand part before it may be used in the right hand parts of the subsequent equations (just like the principle of declaring constants in Pascal). The right hand parts are arbitrary arithmetic expressions built of the Delphi predefined functions, the four arithmetic operations (+, -, *, /), raising to power (^) and factorial (!).


          First, we applied a recursive parsing procedure to each line to obtain its postfix sequence, whose operands are either numeric values, or references to the previously computed variables. The evaluation process of the postfix sequence consists of locating triplets (Operand1,  Operand2, Operation) in it, evaluating the elementary equation


Result := Operand1 Operand2


and replacing the triplets with the corresponding Result until the whole sequence reduces to just one term – the final result of evaluation.


In the previous article, to evaluate a postfix sequence, we followed an approach based on a recursive function GetResult, whose key statement looked like this:


function GetResult : extended;


. . . . . .

GetResult := MyFunctions(GetResult, GetResult);

. . . . . .


This method was quite efficient, and it preserved the postfix sequence for repetitive usage in massive computations. Note, that the results of all intermediate simple evaluations were temporarily stored in a growing stack, thus it was neither possible nor necessary to access them.


The modern Taylor Method means basically iterative evaluation of recurrent equations similar to (1).


Let us consider the standard Initial Value Problem for ODEs:


uk = ck  ,      while  t = t0 ,                             (Initial values)       


uk' = fk(u1, … um); k = 1,2,…, m                   (ODEs)


where uk' means derivatives with respect to t. The equations (2) allow to obtain the set of all first order derivatives uk'  by computing the right hand parts. According to the classical Taylor method, these equations allow to obtain the set of all n-order derivatives uk(n) also. To do that, we have to recurrently apply the well known rules of differentiation to both parts of the equations and substitute the already obtained numerical values. As such, the classical Taylor method had not been considered as a numerical method for computers at least for two reasons: applying the formulas for n-th derivative to complex expressions (superposition of several functions in the right hand parts) is both challenging and computationally wasteful due to the exponentially growing volume of calculations in general case).


On the contrary, the Modern Taylor Method capitalizes on the fact that whatever complicated the right hand parts are, they are equivalent to a sequence of simple formulas (instructions) involving one of the four arithmetical operations and several well known simple functions over no more than two variables. (Of course, this sequence of instructions is exactly that resulting from the process of the postfix evaluation). For them the rule of computing n-th derivative requires only no more than O(n2) operations and it is easily programmable. Unlike to its classical counterpart, not only the Modern Taylor Method becomes viable and competitive for a computer implementation, it also demonstrates several unique features.


In Modern Taylor Method the concept of n-th  derivative means the so called normalized derivative, which by definition is


u[n] = u(n)/n!


so that the Taylor expansion looks like this:


u(t) = S u[k](t - t0)k 


(farther on, the conventional notation u(k) is used, although we always mean the normalized derivatives). For example, the Leibnitz formula for n-th derivative of a product of two functions (it requires only O(n2) operations) is this:



(uv)(n) =S u(i)v(n-i)                                                                                                                                       (3)



A couple of other similarly efficient formulas are the following:



if R=u/v, then R(n) = ( u(n)   -    S R(i)u(n-i) )/v  ;                                    (4)



if R=ua, then R(n) =  (   S(a(1-i/n) - i/n)R(i)u(n-i) )/u,  (a is a constant).



The other similar formulas exist for the natural exponent, logarithm, trigonometric and other elementary functions. (For the linear instructions like a*u, u+v the corresponding formulas are trivial and require O(n) operations).


2. How It Is Implemented


To achieve more flexibility in specifying an Initial Value problem (2), the source data is organized in four sections in the given order: Symbolic constants, Initial values, Auxiliary variables, ODEs (Table 1).



Line number

Variables and equations





. . . . . .

c1 = NumericValue

c2 =C2(c1, NumericValues)

. . . . . . .

. . . . . . .

Symbolic constants


. . . . .


u1 = U1(c1, c2,…, NumericValues)

. . . . . . . .

um = Um(c1, c2,…, u1,…, NumericValues)

Initial values for

the Main variables


v1 = V1(c1,…, um, NumericValues)

v2 = V2(c1,…, um, v1, NumericValues)

. . . . . . . .

Auxiliary variables


. . . . . .


u1’ = f1(u1,…, um, c1, c2,…, v1, v2,…)

. . . . . . . . . .

um’ = fm(u1,…, um, c1, c2,…, v1, v2,…)

The source ODEs


. . . . .



auxiliary variables


Table 1.


The Symbolic constants and Auxiliary variables allow to parameterize Initial Value problems and to optimize computation by eliminating re-calculations of common sub-expressions. For example, in the Three Body Problem the section of Constants contains the equations: 


i15 = -1.5 {exponent index used in subsequent equations}

m1 = 1             {masses} 

m2 = 1

m3 = 1     {assumed > 0 to be used in the formulas below}

x1c = 1    {initial positions of the three bodies}

y1c = 0

x2c = cos(120)

y2c = sin(120)

x3c = cos(240)

y3c = sin(240)

vx1c = 0       {initial velocities of the three bodies}

vy1c = 1

vx2c = cos(210)

vy2c = sin(210)

vx3c = -(m2*vx2c + m1*vx1c)/m3       {these conditions make the sum of moments zero}

vy3c = -(m2*vy2c + m1*vy1c)/m3       {so that the center of masses remains motionless}


The equations in the section Auxiliary variables are these:


dx12 = x1 - x2                 {projections of the distances between the bodies}

dy12 = y1 - y2

dx23 = x2 - x3

dy23 = y2 - y3

dx31 = x3 - x1

dy31 = y3 - y1

r12 = (dx12^2 + dy12^2)^i15      {common sub-expressions of the ODEs}

r23 = (dx23^2 + dy23^2)^i15

r31 = (dx31^2 + dy31^2)^i15


The ODEs themselves are:


t' = 1

x1' = vx1

y1' = vy1

x2' = vx2

y2' = vy2

x3' = vx3

y3' = vy3

vx1' = m3*dx31*r31 - m2*dx12*r12

vy1' = m3*dy31*r31 - m2*dy12*r12

vx2' = m1*dx12*r12 - m3*dx23*r23

vy2' = m1*dy12*r12 - m3*dy23*r23

vx3' = m2*dx23*r23 - m1*dx31*r31

vy3' = m2*dy23*r23 - m1*dy31*r31


If not the auxiliary variables, the more familiar, but non-optimized equations would look like this (just the last one)


vy3' = m2*(y2-y3)/ ((x2-x3)2 + (y2-y3)2) -3/2 + m1*(y1-y3)/ ((x3-x1)2 + (y3-y1)2)-3/2 .



The information of any Initial Value Problem as presented in Table 1 (Strings of TStringList) is to be processed in the following three stages: 


1. Like the list of linear elaborated declarations, it should be parsed and evaluated line by line (units Pars and Postfix), preparing the data for the subsequent differentiation process. 


2. The differentiation process up to the specified order is performed. The set of coefficients of the Taylor expansion (also called the Analytical elements) for each variable are obtained (unit Emulator). 


3. With the obtained Taylor coefficients, the convergence radius R should be estimated, an integration step h<R be determined and the Taylor expansion should be applied to increment the all main variables (units Emulator, AuxMath). Then, the stages 2 and 3 may be repeated as many times as necessary (unless a singularity is reached).


The syntax checking and parsing for the stage 1 is done just like in the case of an arbitrary arithmetic expression (with the same units – DI 8/2000). But since that point on, the goal becomes different. This time we do need the all intermediate results while computing arithmetic expression: we need all those results together with their derivatives up to the required order. Thus, the recursive evaluation GetResult of the postfix sequences is not applicable here. Instead, we will “compile” the postfix sequences into the three-address instructions of an abstract Processor, and create the software implementation of that Processor. 


3. The Evaluating and Differentiating Processor


The memory for this Processor will be the special dynamic array of auxiliary variables AuxVars:


TAnalElem = array of extended;  {to store u, u', u",   }

. . . . . .

AuxVars : array of TAnalElem;  {auxiliary vars and constants: non-rectangular array}


An element AuxVars[i] corresponds to the variable or symbolic constant in line i of the Table 1, while AuxVars[i, j] means the j-order derivative of that variable (when non-constant). What follows after the section ODEs in Table 1 may be either constants, or variables obtained (in the given order) as the intermediate results of parsing and evaluation. Therefore, some of AuxVars[i] are constants and set by SetLength(AuxVars[i], 0), while the others are variables and set by SetLength(AuxVars[i], MaxOrder). Thus, AuxVars is essentially non-rectangular array.


          The format of the instruction of the Processor is this:


type TInstruction = record

       ty : (cc, cv, vc, vv); {flags indicating the nature of the operands: (C)onstant  or (V)ariable}

        op : char;                {operation to perform}

       i1, i2, i3 : word      {addresses, i.e. indexes in AuxVars}



For each postfix triplet (Operand1, Operand2, Operation) there is an instruction of the type TInstruction, corresponding such a way, that Operand1 is in AuxVars[i1], Operand2 in AuxVars[i2], and the result of performing of the Operation is assigned to AuxVars[i3]. To compile a “program” evaluating a postfix sequence, we should scan it from left to right until a first triad (Operand1, Operand2, Operation) is encountered. The corresponding new instruction then should be added to the list with i3 being the index of the operation result in AuxVars. Then, the postfix triplet is substituted by this new result as an operand, thus the sequence shortens for 2 elements. This iterative (not recursive) process is repeated until the postfix sequence reduces to just one operand – the result of the evaluation. So, the original postfix sequence is destroyed and no more available for repetitive computations. Instead, it was encoded in a set of instructions, preserved for the subsequent massive computations of the corresponding arithmetic expression.


Actually, this Processor is designed to do not only the simple evaluation like


Result := Operand1 Operand2,


but also to perform the iterative differentiation up to any given order 


Result(n) := (Operand1 Operand2)(n), n = 0, 1, 2,…


given the corresponding formulas for computing n-th derivative of certain functions – for example the Leibnitz formula (3) and the others (4) mentioned previously (see Code Sample). 


So, according to the Stage 2, with this Processor we can obtain derivatives of any required order in the initial point t0 for all variables (Table 1) representing the source initial value problem. 


To estimate the convergence radius, we use the algorithm based on the formula of Cauchy-Hadamard


R -1 = lim |an|1/n


where Taylor coefficients an are obtained for each unknown function in the source ODEs. Each component of the solution  may have its own radius, thus the procedure selects the minimal among them, which should be considered as a convergence radius for the whole system. Although the formula is exact, the result is heuristic because the formula is applied to only a finite number of terms, with simplified assumptions for obtaining the upper limit. The fraction powers are obtained by a simplified, coarse, but very fast method using just the exponents in the float point representation. 


So, given the Taylor coefficients a0, a1, …, an, first of all we obtain their float point exponents ea0, ea1, …, ean using the standard procedure FrExp(x, m, result) for all non-zero coefficients, and assigning ‑4951 (the minimal possible extended exponent) to zero coefficients. These are the coarse log2 of the values, so that logarithms of the sequence of the required fraction powers for farther processing are these: ea1/1, …, ean/n=sn. To find the greatest point of condensation for this finite sequence (usually 30-40 terms), the following heuristic procedure is applied.


          The sequence {sn} is sorted in descending order (suppose, it is already done). Assuming as a density segment width w a value like n/6 (usually 4-5), we start scanning {sn} for the most narrow difference si - si+w , i =1, 2,…, meaning that the greatest “condensation” point must be probably there. Thus, the i,  delivering that minimum is farther considered for obtaining the average s of si ,…, si+w, and then (returning from the binary logarithms back to the basic values) the value 1/2s is accepted as a heuristic radius of convergence R.


To meet the given requirements to the accuracy of the integration, there exist different strategies of selecting the integration step h<R and the order of differentiation n. In this Taylor Center prototype the default order is 30 (it may be changed), and the only implemented strategy is such that h=R/2.


4. Graphics


One of the important goals of this application (along with numeric integration and Taylor expansion) is graphing the solution. Given the multi-component solution, we may wish to draw different projection of the solution vector onto different coordinate planes (or cubes). Although this prototype implements only 2D graphics, the next version is scheduled to realize the 3D curves imaging on conventional display, creating red/blue stereo pairs to be watched through anaglyph (blue/red) glasses.


Formally, any pair of the solution components {ui(t), uj(t)} called a trajectory may be considered for graphing in the coordinates Oui ,Ouj , although just few of the all possible pairs may be meaningful. For example, the trajectory of motion {x1(t), y1(t)} of a certain body is meaningful, but, say {x1(t), y2(t)} is not. (In the case of 3D imaging, we would like to consider different triplets {ui(t), uj(t), uk(t)} of the components, in particular the trajectories like {x1(t), y2(t), z(t)}). The application displays the special grid with the all available variables, giving the user an opportunity to specify the required trajectories.


The important particularity of the Taylor integration is that we obtain the solution not as a table of values with a small regular step, but rather as a set of Taylor expansions covering the required domain with a finite non-regular grid t0, t1,…, tn. In certain regions the step may be really big (measuring by the scale of the whole domain), while in others it diminishes (depending on the properties of the solution). Correspondingly, the node points of a trajectory {x(tk), y(tk)} (plotted on the display with a certain scaling factor), would appear with gaps in one place and over-condensed in the others. Thus, to draw the continues curves on the pixels grid we need to design an efficient algorithm of filling described here.


The curve drawing in this application may be done either just for the visualizing the current integration step (when the user clicks the button), or in the so called Play mode, when the curves evolve with the speed proportional to the real velocities {ui'(t), uj'(t)} of the corresponding component. This is achieved by adding a short-looped delay proportional to Dt, controlled by the precise hardware counter QueryPerformanceCounter (or just by simple and coarse Sleep procedure, if the former doesn't work.


5. The Features of the Taylor Center


This DEMO application is intended for numerical integration, exploring and graphing the solutions of ODEs (in particular, it allows to obtain the Taylor series expansions of any given order for these solutions). The DEMO may be downloaded at  



With the current version of the product you can:


5. The Conclusion


The Taylor Center is a powerful numeric (and expression-processing) application with extensive visual features. It demonstrates the usage of both recurrent parsing and iterative procedures such as evaluation and encoding the postfix sequences into instructions of the abstract Processor. The processor is capable to compute not only expressions, but also the n-order derivatives of the expressions. 


It is probably the first such system for PC under Windows, and there is an easy opportunity of porting it to Linux platform as well.


Using only the parsing units developed for the Advanced Calculator, this project took mere three months of after-the-office hours to grow into a full featured math-intensive and graphical application. If not the powerful features and RAD environment of Delphi, it would be hardly possible to do it like that. This application essentially capitalizes on the non-rectangular and rectangular Dynamic arrays (introduced since Delphi‑4), which proved to be very efficient. Specifically in this application, the zero-base indexing of the Dynamic arrays fortunately did not conflict with the mathematical meaning and the requirements of the algorithm, but that not necessarily happens for other applications. That is why the necessity of the general Pascal indexing also for the Dynamic arrays has been brought to the attention of the Delphi team and published in ACM SIGPLAN Notices [5].


The further reading on the Modern Taylor Method:


          1. Moore R. E., Interval Analysis. Prentice-Hall: Englewood Cliffs, NJ (1966)


          2. Chang Y. F., The ATOMFT User Manual. Version 2.50, 976 W. Foothill Blvd . #298, Claremont, CA 91711 (1988)


          3. Chang Y. F., Solving Stiff Systems by Taylor Series. Appl. Math. Comp., 31, Special Issue, 251-269 (1989)


          4. Gofen A. M., Taylor Method of Integrating ODEs: the Problem of Steps and Singularities. Kosmicheskie Issledovaniya (Cosmic Research), Vol. 30, No. 6, pp.723-737, (1992).


5. Gofen A., From Pascal to Delphi to Object Pascal-2000. ACM SIGPLAN Notices, Vol. 36, Num. 6, June 2001, pp 38-49. Also found at