**The Modern Taylor Method**** Solver of**** ODEs:
The **

*Keywords*: ODE Solver, Taylor Solver, Taylor Method,
Taylor
Series, Numeric
Integration, High Accuracy, Ultimate Accuracy, Finite Step,
Dynamic Graphics, 3D Stereo, 3D Cursor, Trajectories, Interactive
environment,
Windows,
Linux.

The Modern Taylor Method is a
descender of
its classical counterpart. It is an efficient method for numerical
integration
of the Initial Value Problems for Ordinary Differential Equations
(ODEs). What
distinguishes it from all other numerical methods for ODEs is that only
the
Taylor Method can compute the increments of the solution with
principally
unlimited order of approximation so that the integration step __does
not__
approach zero whichever high accuracy is specified. That is possible
because
the method performs the *automatic differentiation* - exact
computing of
the derivatives up to any desired order *N*, allowing to obtain
the Taylor series
of any length for the solution components.

To download the (32-bit) Demo for
Windows
click here, then download and unzip
the file
("*Save*", don't "*Open*" it in your browser).
Unzip and keep it in an empty folder, *TCenter.exe* being the
only executable
to run. Preserve these files and sub-folders structure (in order that
the
program work properly). See a short guide
navigating you through the DEMO. You can also download the full User Manual (in
MS Word
doc format), or the articles (in pdf) published in the
Proceedings CSC 2005
, DMS 2007 , and
in CODEE
Journal, September 2012.

From the algorithmic point of view, this software parses the right hand sides of the ODEs and auxiliary equations, and compiles them into a sequence of pseudo instructions of Automatic Differentiation. Then the programmatic emulator of those instructions runs them performing the evaluation of the derivatives and integration of the Initial Value Problem.

To make inquires, please contact Alexander Gofen at galex@ski.org .

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 as necessary);

- 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
^{4932}implied by the Intel generic 10-byte float point type*extended*with 64-bit mantissa (contrary to the reduced 8-byte 48 bit mantissa in Microsoft C++ for this type); - Graph terms of 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 (in 2D and 3D) 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**defining the same trajectory, but with respect to different independent variables. For example, it is possible to switch the integration in respect to*several states*of ODEs*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**for an 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*array*of Initial Value Problems*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 30) 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 (in 2D, or a triplet in 3D) of variables of the solution –
up to 200 on one screen – either as plane projections, or as 3D
stereo images 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 looking like "hanging in thin air". At that...

- 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 enhanced 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, along with 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 accepting the real 80 byte type extended) 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.
This provides a way to re-use the integration results obtained in the
Taylor Center in other user applications;

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

In particular, the Demo includes a fascinating example of the so called Choreography for the Three Body motion, an eight-shaped orbit, discovered just recently by Chenciner and Montgomery (2000) . Click here to learn more about the Choreographies of N-body problem. You can "feed" the equations found there into the Taylor Center, integrate them, draw the curves and play the motion in the real-time mode all in the same place. Another recent fascinating example of the four body non-planar trajectories inscribed in a cube discovered by Cris Moore & Michael Nauenberg (also here) is incorporated too.

The Taylor Center is available both as a 32-bit and 64-bit
application
for Windows. The executable module is only 1 Mb (32-bit version) and
near 3 Mb
(64-bit). As a 32-bit application, the program can use no more than 4
Gb of
available memory for variables and their expansions –
the limit far exceeding
any practical needs. (The new 64-bit
version does away with this 4 Gb
limitation, but meanwhile the compiler unfortunately maps 10 byte float
point
numbers into
the 8 byte type – due to
the limitation of the currently existing 64-bit Delphi
compiler, mimicking the Win64 blunder - using SSE instead of the x87
FPU).

The memory consumption depends on the number of variables *VarNum*
(a
function of the number of ODEs and their complexity) and on the
specified *Order*
of approximation. If the expansions are not stored, the program takes *2*VarNum*Order*10*
bytes of memory. If the expansions in *P* points are stored, it
additionally requires *P*VarNum*Order*10* bytes.

Here is a benchmark for the 10 body planar problem comprised of
2*(10+10)+1=41
ODEs and 10*9/2*(2+1)=135 auxiliary equations, which are parsed into
811 AD
instructions. For this problem 10000 steps of integration took 32 s (or
3.2 ms
per step) on 2.4 GHz Pentium with polynomial expressions spelled out,
and it
took 29 s (or 2.9 ms per step) with polynomial expressions encoded.

For a benchmark for the 99 body spatial problem comprised of
3*(99+99)+1=595
ODEs and 99*98/2*(3+1)=19404 auxiliary equations, 100 steps took 65 s
with
polynomial expressions spelled out, and took 42.5 s with polynomial
expressions
encoded.

Generally for each system of ODEs there exists such a small value of
the
accuracy requirement, that at this high accuracy the Taylor methods
beats
any
fixed order method due to the unlimited order of approximation in the
Taylor
method.

The future version will include the following:

- It will be supplied not only as the Taylor Center GUI executable, but also as the separate Delphi component (to include them directly in Delphi projects);

- It will implement the Merge procedure and a library of ODEs – definitions of a large variety of commonly used elementary functions. (Presently, the functions which are not in the allowed list, may be used also – providing that the user declares the ODEs defining them and properly links them with the source ODEs (more about that in Help for Merge). Also, it will include a larger variety of the "calculator" functions for specifying even more complex relationships between constants, parameters and the initial values.

- The application will be ported to windowed Linux;

- The set of the internal differentiation instructions will be translated into the machine code – to reach the highest possible speed for massive computations. (Meanwhile it is an emulator written in Delphi which runs these instructions). Also, the set of the inernal AD instructions may be translated into instructions in Pascal, C or Fortran to be further compiled and linked with other applications;

Here is how the front panel of the Taylor Center looks like:

And here are the trajectories of a slightly disturbed (Lagrange) case of the Three Body Problem (Demo/Three bodies/Disturbed/2D). However there is a dramatic difference between viewing a still image like in this page vs. dynamically evolving real time motion displayed in the running Taylor Center. In the Demo you will be fascinated to watch all possible pairs of the three bodies coupling randomly in turn with acceleration or deceleration.

Here is a Double spiral (*Demo/Spirals/Double
spiral*) – a solution of the system

*x' = x² - y²+
2xy - x -3.5y +
1, x(0)=0.1
y' = -x² + y² + 2xy - y +3.5x - 1,
y(0)=0 *

corresponding to one complex ODE *z' = (1 - i)z² + (-1 +
3.5i)z + 1 -
i, z = x + iy. *Similarly to the previous
example, it is
one thing to see this Double Spiral as* *a static image

It is quite different thing to watch the *real time* motion along
this
spiral observing that every lap of the spiral takes exactly the same
time to run
(with a huge acceleration in big laps), as it follows from the complex
analysis.

And here is an array of 5 Initial Value Problems integrated
simultaneously for
the same system of ODEs with various initial values:

The unique feature of this software is that 3D viewing and real time animation is implemented also in stereo mode. To view stereo images, you need anaglyphic Red/Blue glasses (Left - Red, Right - Blue). Such glasses may be ordered here: look for a Red/Blue type (not Red/Cyan).

If you have the anaglyphic glasses, click to view 3D stereo images. You may view non-planar curves also in conventional axonometry projection (where their spatial properties are lost). However, in order to improve axonometric view of such curves, this software can draw them as though tubes of a required thickness implementing the proper skew resolution at points of illusory intersections like in this family of curves representing the Möbius surface outline:

You may also construct and display
an enhanced field of
directions or a phase portrait
for selected pairs of variables built not with small vectors, but
rather with* analytic elements*
i.e. curvy strokes representing finite pieces of trajectories obtained
in
multi-step integration at points of a grid or at arbitrary points of
your
choice. Due to the high order of the Taylor method, these strokes are typically long
enough to
represent the near continuous family of the solutions or the phase
portrait
like that below for the system:

*x' = -y + *2*x ^{2} + xy - y^{2}*

Here is another beautiful phase portrait made by this software for a
system (credit to John H. Hubbard and Beverly H. West)

x' = 1

y' = y^{2} - x

The image below depicts the two curves of this phase portrait near
the separatrix. The initial values of these curves differ only in the
last binary digit of their 64 bit
mantissa:

Normally, in order to populate and
reach the maximal accuracy in 64 bit mantissa, 19 or 20
significant decimal digits suffice. The routine converting real decimal
numbers into this float binary format may take also extra decimal
digits, yet the binary result will be rounded to 64 significant binary
digits anyway. In this example the two decimal numbers are represented
with 22 decimal digits, however the difference between their binary
form is exactly one last bit of the mantissa (the goal set in this
experiment).

This smallest difference in the
initial values however was enough to finally direct the black curve
upward yet
the red curve downward, demonstrating the ultimate accuracy of
integration for this 64-bit mantissa real type. The relative error
tolerance in this case was set to 10^{-25} < 10^{-19} meaning that at every of 80 finite steps,
all 64 binary digits of the result were true, so that the highest
accuracy at every step (in a format with 64 bit mantissa) was achieved.

For pedagogical aspects of the Taylor Center look here: