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₀ = f₀({numeric values only});

u₁ = f₁(u₀, {numeric values});                                            (1)

. . . .

u_n = f_n(u₀, u₁,…, u_n-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;

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₀ ,                             (Initial values)       

                                                                                      (2)

u_k' = f_k(u₁, … 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 n-order derivatives u_k⁽ⁿ⁾ 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:

  n

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

                i=0

A couple of other similarly efficient formulas are the following:

n-1

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

                                                                i=0

n-1

if R=u^a, then R(n) = (   S(a(1-i/n) - i/n)R(i)u(n-i) )/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+m-1	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+m-1	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 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)² + (y2-y3)²) ^-3/2 + m1*(y1-y3)/ ((x3-x1)² + (y3-y1)²)^-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}

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⁽ⁿ⁾ := (Operand1 ° Operand2)⁽ⁿ⁾, 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 t₀ 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 log₂ 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 30-40 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 4-5), 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 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 {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 non-regular grid t₀, t₁,…, 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 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 {u_i'(t), u_j'(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  http://TaylorCenter.org/gofen/TaylorMethod.htm

With the current version of the product you can:

• Specify and study the Initial Value Problems for virtually any system of ODEs in the standard format, meaning a system of explicit first order ODEs, derivatives in the left hand sides and arithmetic expressions in the right hand. The standard elementary functions, numeric and symbolic constants and parameters may be used;

• Enter arithmetic expression in the standard Pascal syntax either through the editor windows, or via the Polynomial Designer for cumbersome polynomial expressions;

• Perform numerical integration of Initial Value Problems with an arbitrary high accuracy along a path without singularities, while the step of integration remains finite and does not approach zero (presuming the order of approximation or the number of terms could increase to infinity with the length of mantissa unlimited);

• Apply an arbitrary high order of approximation (by default 30), and get the solution in the form of the set of analytical elements - Taylor expansions covering the required domain;

• Study Taylor expansions and the radius of convergence for the solution at all points of interest up to any high order. An upper limit for the terms in the series is as high as 10⁴⁹³² implied by the Intel generic 10-byte float point type extended with 63-bit mantissa (contrary to the reduced 8-byte 48 bit mantissa in Microsoft C++ for this type);
• Graph Taylor expansions as bar diagrams and vary the step h observing its effect on the bell shape bulge;
  

• Perform integration either "blindly" (observing only the numerical changes), or graphically visualized. The visualization comprises graphs of the solution, real time motion along the trajectories, and the field of directions, or rather a phase portrait comprised of particular families of trajectories;

• Terminate the integration either after a given number of steps, or when an independent variable reaches a given terminal value, or when a dependent variable reaches a given terminal value (as explained in the next item);

• Automatically generate the ODEs and switch integration between several states of ODEs defining the same trajectory, but with respect to different independent variables. For example, it is possible to switch the integration in respect to t to that by x, or by y in order to reach the terminal value (or zeros) of a former dependent variable (x, or y). In particular, if the initial value is a nonzero (guess) value and the terminal value is set to zero, the root (the zero) of the solution may be obtained directly without iterations;

• Automatically generate and integrate an array of Initial Value Problems for the array of initial vectors. The solutions of these IVPs are displayed in one plot resembling a phase portrait. An array of IVPs considered as an IVP with an indefinite parameter helps to estimate the solution of certain boundary value problems;

• Integrate piecewise-analytical ODEs;
• Integrate IVPs in complex variables along an arbitrary path in a complex plane –  by automatic (since version 29) transformation of the source ODEs into a system over the real and imaginary parts of the respective variables;

• Specify different methods to control the accuracy and the step size;

• Specify accuracy for individual components either as an absolute or relative error tolerance, or both;

• Graph color curves (trajectories) for any pair of variables of the solution – up to 99 on one screen – either as plane projections, or as 3D stereo images (for triplets of variables) to be viewed through anaglyphic (Red/Blue) glasses. The 3D cursor (controlled by a conventional mouse) with audio feedback enables "tactile" exploration of the curves virtually "hanging in thin air";
• Graph non-planar curves as though tubes of a required thickness implementing the proper skew resolution at points of illusory intersections;

• Play dynamically the near-real time motion along the computed trajectories either as 2D or 3D stereo animation of moving bullets;

• Graph the field of directions – actually the field of curvy strokes, whose length is proportional to the radius of convergence.

• Explore several meaningful examples supplied with the package such as the problem of Three and Four Bodies. Symbolic constants and expressions allow parameterization of the equations and initial values, and trying different initial configurations of special interest.

• Automatically generate ODEs for the classical Newtonian n-body problem for n up to 99, and then integrate and explore the motion. In the case of n=99 there are 595 ODEs, 19404 auxiliary equations, compiled into over 132000 variables and over 130000 AD processor's instructions: a "heavy duty" integration!
  
• Use a DLL (for any programming languages under Windows) which provides functions computing the vector-function of a solution of ODEs (first obtained with the Taylor Center, but later – without using the Taylor Center). Instead, this DLL implements the optimized (Horner algorithm) for computation of the desired vector function using the polynomial expansions obtained from the Taylor integration and saved in a binary file.

• Integrate a few special instances of singular ODEs having regular solutions at the points of the so called "regular singularities".

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