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:
u_{0 }= f_{0}({numeric values only});
u_{1 }= f_{1}(u_{0}, {numeric values}); (1)
. _{ }. . .
u_{n }= f_{n}(u_{0}, u_{1},…, u_{n1}, {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;
begin
. . . . . .
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:
u_{k} = c_{k} , while t = t_{0} , (Initial values)
(2)
u_{k}' = f_{k}(u_{1}, … u_{m}); k = 1,2,…, m (ODEs)
where u_{k}' means derivatives with respect to t. The equations (2) allow to obtain the set of all first order derivatives u_{k}' by computing the right hand parts. According to the classical Taylor method, these equations allow to obtain the set of all norder derivatives u_{k}^{(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 nth 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 nth 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 nth 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 nth derivative of a product of two functions (it requires only O(n2) operations) is this:
n
(uv)(n) =S u(i)v(ni) (3)
i=0
A couple of other similarly efficient formulas are the following:
n1
if R=u/v, then R(n) = ( u(n)  S R(i)u(ni) )/v ; (4)
i=0
n1
if R=u^{a}, then R(n) = ( S(a(1i/n)  i/n)R(i)u(ni) )/u, (a is a constant).
i=0
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 
Explanation 
0 1 2 . . . . . . 
c1 = NumericValue c2 =C2(c1, NumericValues) . . . . . . . . . . . . . . 
Symbolic constants 
MainStart . . . . . MainStart+m1 
u1 = U1(c1, c2,…, NumericValues) . . . . . . . . um = Um(c1, c2,…, u1,…, NumericValues) 
Initial values for the Main variables 
MainStart+m 
v1 = V1(c1,…, um, NumericValues) v2 = V2(c1,…, um, v1, NumericValues) . . . . . . . . 
Auxiliary variables 
DifStart . . . . . . DifStart+m1 
u1’ = f1(u1,…, um, c1, c2,…, v1, v2,…) . . . . . . . . . . um’ = fm(u1,…, um, c1, c2,…, v1, v2,…) 
The source ODEs 
DifStart+m . . . . . 

Internal auxiliary variables 
Table 1.
The Symbolic constants and Auxiliary variables allow to parameterize Initial Value problems and to optimize computation by eliminating recalculations of common subexpressions. 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 subexpressions 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 nonoptimized equations would look like this (just the last one)
vy3' = m2*(y2y3)/ ((x2x3)^{2} + (y2y3)^{2})^{ 3/2} + m1*(y1y3)/ ((x3x1)^{2} + (y3y1)^{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 threeaddress 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: nonrectangular 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 jorder derivative of that variable (when nonconstant). 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 nonrectangular 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}
end;
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 nth 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 t_{0} 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 CauchyHadamard
___
R 1 = lim an1/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 nonzero coefficients, and assigning ‑4951 (the minimal possible extended exponent) to zero coefficients. These are the coarse log_{2} of the values, so that logarithms of the sequence of the required fraction powers for farther processing are these: ea1/1, …, ean/n=s_{n}. To find the greatest point of condensation for this finite sequence (usually 3040 terms), the following heuristic procedure is applied.
The sequence {s_{n}} is sorted in descending order (suppose, it is already done). Assuming as a density segment width w a value like n/6 (usually 45), we start scanning {s_{n}} for the most narrow difference s_{i}  s_{i+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 s_{i} ,…, s_{i+w}, and then (returning from the binary logarithms back to the basic values) the value 1/2^{s }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 multicomponent 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 {u_{i}(t), u_{j}(t)} called a trajectory may be considered for graphing in the coordinates Ou_{i },Ou_{j}_{ }, 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 {u_{i}(t), u_{j}(t), u_{k}(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 nonregular grid t_{0},
t_{1},…, t_{n}. 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(t_{k}), y(t_{k})}
(plotted on the display with a certain scaling factor), would appear
with gaps
in one place and overcondensed 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 {u_{i}'(t), u_{j}'(t)} of the corresponding component. This is achieved by adding a shortlooped 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 http://TaylorCenter.org/gofen/TaylorMethod.htm
With the current version of the product you can:
5. The Conclusion
The Taylor Center is a powerful numeric (and expressionprocessing) 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 norder 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 aftertheoffice hours to grow into a full featured mathintensive 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 nonrectangular and rectangular Dynamic arrays (introduced since Delphi‑4), which proved to be very efficient. Specifically in this application, the zerobase 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. PrenticeHall: 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, 251269 (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.723737, (1992).
5. Gofen A., From Pascal to Delphi to Object Pascal2000. ACM SIGPLAN Notices, Vol. 36, Num. 6, June 2001, pp 3849. Also found at