The Modern Taylor Method Solver of ODEs:
The Taylor Center


Keywords: ODE Solver, Taylor Solver, Taylor Method, Taylor Series, Power Series Method, Differential Transform Method, Adjoint equations, 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 (IVPs) for holomorphic Ordinary Differential Equations (ODEs) whose Right Hand Sides are general elementary functions. What distinguishes it from all other numerical methods for ODEs is that …

·        The Taylor method computes the increments of the solution with principally unlimited order of approximation so that the integration step must 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.

·        The Taylor method does not use any finite difference formulas such as f'(t) ≈ ( f(t+h) – f(t) ) / h  prone to (catastrophic) subtraction error with float point numbers of fixed length.  

·        The Taylor Center utilizes the 80-bit float point type extended with 63-bit mantissa generic for the processor X87 FPU while other numeric programs (such as Mathematica®, Maple®, Matlab®, Simulink® with MathWorks®) use the 64-bit float type double with 52-bit mantissa for their generic float point computation. (The emulated unlimited precision in those system is several orders of magnitude slower than the generic float point processing).

To download the Demo for Windows click here, then download the file: Save, rather than Open it in your browser. Unzip it entirely (right click choosing Extract all) 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 2005DMS 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 evaluating the derivatives and integrating the Initial Value Problem.

To make inquiries, please contact Alexander Gofen at .

With the current version of the product you can:

In particular, the Demo includes fascinating examples of the so called choreography for the 3 body problem: 345 of them (courtesy of Dr. Carles Simò), plus 204 cases of periodic orbits of unusual shapes (courtesy of Ana Hudomal). Click here to learn more about the Choreographies of N-body problem and how to "feed" their ODEs and initial values into the Taylor Center, plot the curves and play the motion in the real-time mode: all in the same place. Similarly, click here for playing with the periodic orbits for the 3 body problem from the list of Ana Hudomal. These orbits are closed curves (as intuitively expected from periodic orbits). Here however you may see the newly discovered periodic orbits which are finite curved segments, whose extremes are resting points in the 3-body problem, so that the bodies periodically fulfill a free fall along these segmented orbits (the data was kindly provided by Xiaoming LI and Shijun LIAO). 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. Here is the list of brief explanations for all several hundreds samples pre-loaded with the program. 

As of now (and in foreseeable future), the Taylor Center will remain a 32-bit application run on the x87 FPU. This processor was designed to address only 32 bit address space, i.e. no more than 4 Gb memory, or 400 millions of variables and their expansions (as 10 byte float point numbers).


In order to randomly address a memory space larger than 4 Gb, Intel (and other companies) enhanced the x87 FPU processor by adding a set of instruction SSE capable to randomly access a 64 bit address space. However in doing so, Intel designed the SSE instructions to operate only on the 8 byte data types abandoning the 10 byte type extended. Though it is logically possible to expand the addressing ability of the x87 FPU by combining its 32-bit instructions with the SSE instructions to randomly address the 64 bit memory space, such a trick prohibitively slows the overall operation speed of x87 FPU on the 10 byte type.


For the program like the Taylor Center it's crucial to operated with the highest precision available in a processor. Therefore, with such a design blunder by Intel, it makes sense to maintain the Taylor Center only as a 32 bit application operating at the x87 FPU with the 10 byte type at the highest speed.  


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.

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.

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.

An accuracy benchmark for the so-called Pythagorean triangle 3-body problem (Pythagor.scr) demonstrates how successfully this software overcomes near catastrophic cancellation errors emerging during close approach of the bodies in this particular problem. While fixed order methods required the Levi-Civita's regularization of binary collisions in order to complete the integration, this software successfully integrated this problem without regularization as a general Newtonian 3-body problem.

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:

    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.

3  Bodies 2D

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 (and explanation of the optimal viewing conditions). 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 double Möbius surface outline:

This software offers several tricks for outlining 3D surfaces. The image above is based on visualizing one of two families of the internal coordinate lines of the Möbius surface. Another way of visualizing single Möbius surface is based on the sine wave with a short wavelength: here in 3D stereo or in tubular graphic below

One more approach of outlining is based on helix, say for outlining tori, or the Klein bottle below    

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 + 2x2 + xy - y2

y' = x + x2 - 3xy


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: