This article deals with a high order integration method based on the Taylor series. The paper shows many positive properties of this method on a set of technical initial value problems. These problems can be transformed into the autonomous systems of ordinary differential equations for both linear and nonlinear problems. The MATLAB implementation of the method is compared with state-of-the-art MATLAB solvers.
The article deals with the solution of technical initial value problems (IVPs). The solution of the technical IVPs means to find the numerical solution of the system of ordinary differential equations (ODEs). In the article, ODEs are solved using a high order variable-step variable-order numerical integration method – Modern Taylor Series Method (abbreviated as MTSM in the article) .
The article consists of several sections. In Section 2, MTSM and its properties is introduced. Subsections 2.1 and 2.2 focus on the solution of linear and nonlinear problems using MTSM, respectively. Several examples of linear and nonlinear problems are introduced in Section 3; namely, movement of a charged particle (Section 3.1), calculation of Fourier coefficients (Section 3.2) and Kepler problem (Section 3.3). The results obtained by MTSM are compared with state-of-the-art MATLAB solvers. This article extends the paper . The new nonlinear problem (Kepler problem) has been added and other experiments were revised. Nonlinear problems and their solutions using the MTSM are analyzed in detail, including step-by-step simplification and transformations of the problem into several autonomous systems that are compared.
The article follows the previously published article , which focused on linear systems of ODEs and their solutions.
2 The Modern Taylor Series Method
with initial condition
is called an initial value problem. The best-known and the most accurate method of calculating a new value of the numerical solution of an ODE is to construct the Taylor series in the form
where h is the size of the integration step, yi = y(ti) is the previous value and yi+1 = y(ti + h) is the next value of the function y(t) .
The MTSM very effectively implements the variable-step-size, variable-order numerical calculation of differential equations using the Taylor series’ . It is based on a recurrent calculation of the Taylor series’ terms for each integration step. Therefore, the complicated calculation of higher-order derivatives does not need to be performed, but rather the value of each Taylor series term can be numerically calculated .
Equation (3) can then be rewritten in the form
where DY denotes the Taylor series’ terms.
The MTSM automatically transforms the input problem into a system of autonomous ODEs, which allows the recurrent calculation of terms of the Taylor series. Another important part of the method is an automatic integration order setting, i.e. using as many Taylor series terms as the defined accuracy requires. Let us denote as ORD the function which changes during the computation and defines the number of Taylor series terms used in the current integration step yi+1 (ORDi+1 = n). The MTSM allows for computation with arbitrary accuracy and step-size if variable-precision arithmetic and higher order of the method is used.
The first implementation of the MTSM was TKSL/386 (TKSL stands for Taylor-Kunovsky Simulation Language) . Currently, the MTSM has been implemented and tested in MATLAB , in C/C++ languages (FOS  and TKSL/C software ). Additionally, the method can be effectively implemented in hardware . Several other implementations of the Taylor series method in a variable order and variable step context were presented by different authors e.g. TIDES software  and TAYLOR , which includes a detailed description of a variable-step-size version. Other implementations based on Taylor series include ATOMF , COSY INFINITY , and DAETS . The variable step-size variable-order scheme is also described in [17, 18], and , where simulations on a parallel computer are shown. An approach based on an approximate formulation of the Taylor methods can be found in .
2.1 Linear Problems
For linear systems of ODEs, Equation (1) is in the form
where A is the constant Jacobian matrix and b is the constant right-hand side. Moreover, Taylor series terms DY in (4) can be computed recurrently using
2.2 Nonlinear Problems
The effective solution of nonlinear systems of ODEs is described. For such a system, a Taylor series based numerical method can be implemented in a very effective way.
Equation (1) for nonlinear systems of ODEs can be rewritten as
where A ∈ ℝne×ne is the matrix representing the linear part of the system and B1 ∈ ℝne×nmjk, B2 ∈ ℝne×nmjkl are the matrices representing the nonlinear part of the system, b ∈ ℝne is the right-hand side for the forces incoming to the system, y0 is a vector of the initial conditions, the symbol ne stands for the number of equations of the system of ODEs and the symbols nmjk, nmjkl mean the number of multiplications. The unknown function yjk ∈ ℝnmjk represents the vector of multiplications yj. * yk and similarly yjkl ∈ ℝnmjkl represent the vector of multiplications yjj. * ykk. * yll, where indexes j, k, jj, kk, ll ∈ (1, . . . , ne) come from multiplications terms in the system (8). The operation ‘.*’ stands for element-by-element multiplication, i.e. yj. * yk is a vector (yj1 yk1, yj2 yk2, . . . , yjnmjk yknmjk)T.
For simplification, the matrices A, B1, B2, . . . and the vector b are constant. The higher derivatives of the terms B1yjk, B2yjkl in (8) can be included in a recurrent calculation of Taylor series’ terms DYB1 and DYB2
where the linear term DYAr−1 is computed using (7). The Taylor series higher order multiplication terms DYB3, DYB4, . . . can be calculated recurrently in a similar way.
3 Selected examples
In this section, several examples of linear and nonlinear problems are presented and solved using the MTSM. The MTSM was implemented in MATLAB software  using vectorization. The results are compared with state-of-the-art non-stiff solvers in MATLAB, specifically ode23, ode45 and ode113 .
All experiments were performed using MATLAB 2015a on the Salomon supercomputer . Note that columns labeled as Time of calculation and Ratio are taken as a median value from 100 runs. Ratios of computation times ratio = ode/MTSM ≫ 1 indicate significantly faster computation using the MTSM.
3.1 Movement of a charged particle
where e = (ex, ey, ez)T is the electric field, b = (bx, by, bz)T is the magnetic field, vector v is instantaneous velocity of the charged particle with charge q. The force vector f can be substituted as ma,
We can further suppose that acceleration a in x, y and z axis (x″, y″, z″) can be calculated as
By substituting y = (x′, y′, z′, x, y, z)T, the problem can be transformed into a system of linear ODEs in a matrix-vector representation (5),
where initial conditions are y0 = (−8 · 107, 0, 0, 0, 0, 0)T. Constants are set to m = 9.10938356 · 10−31 [kg], q = −1.6 · 10−19 [C], bx = by = 0, bz = 0.5 [T] and ex = ey = ez = 0 [T]. Further, the right-hand side . The behavior of the electron is in Figure 1.
Results of simulations using IVP (13) for the given parameters of the electromagnetic and electrostatic fields and time of simulation t ∈ (0, 10−8) are in Table 1. The MTSM (with fixed integration step size h = 10−10) calculates the problem faster than the state-of-the-art MATLAB ode solvers. The order of the method is ORD = 51 ± 1.
|Solver||Time of calculation [s]||Ratio||Number of steps|
3.2 Fourier coefficients
A Fourier series is the limit of the sequence of trigonometric polynomials that have a cosine part and a sine part. They are mainly used in the study of phenomena with a periodic character. The advantage of the Fourier series is that the requirements for its convergence to the developed function is weaker than in the case of the Taylor series. For example, the existence of derivatives of all orders of a given function at a given point is not required. The calculation of coefficients can be (especially when using numerical methods) easier than for the Taylor series.
Any periodic function f (t) can be written as a sum of harmonics (14) called the Fourier series
where T is the period of the function f (t). The calculation of the Fourier series’ coefficients requires the calculation of definite integrals. The definite integral
can be transformed into the differential equation
where y(T) = Y and T denotes the maximum time of simulation t ∈ (0, T). Note that the same approach can be used for the calculation of multiple definite integrals. More details about transformation and numerical solution of definite integrals using the MTSM can be found in .
For example, the analytic calculation of coefficient A2 using the definite integral
can be transformed into the IVP
As an example, consider the simple harmonic function
with parameters , T = 2 [s]. The coefficients of the Fourier series of (22) given by
can be calculated analytically
as A0 = 1 and . The coefficients can be calculated numerically using the following system of ODEs
The solution of the IVP in t ∈ (0, T), where T = 2 [s] is in Figure 2. The final values of functions a0(t) and a2(t) show the calculated values of Fourier coefficients a0(T) = 1, a2(T) = −0.5.
The IVP (25) can be transformed into different autonomous systems with auxiliary generating equations. The first simple substitution y = (a0, a2, sin(ωt), cos(ωt), sin(2ωt), cos(2ωt))T leads to the following nonlinear system of ODEs
The first experiments with (26) with multiplication of 3 terms ( ) in the second equation were performed in . The more effective way of solving the IVP (25) is to transform it into the autonomous system with fewer multiplications or even better to the linear system of ODEs. The next substitution y = (a0, a2, sin2(ωt), sin(ωt) cos(ωt), cos2(ωt), cos(2ωt), sin(2ωt))T leads to the following nonlinear system of ODEs with at most two term multiplication
The last transformation of the IVP (25) using substitution y = (a0, a2, sin2(ωt), sin(ωt) cos(ωt), cos2(ωt), sin2(ωt) cos(2ωt), sin2(ωt) sin(2ωt), sin(ωt) cos(ωt) cos(2ωt), sin(ωt) cos(ωt) sin(2ωt), cos2(ωt) cos(2ωt), cos2(ωt) sin(2ωt))T leads to the linear IVP
Results of the calculation of the Fourier coefficients using IVPs (27) and (28) are in Tables 2 and Table 3, respectively. The tolerance for all solvers was set to 10−10 with fixed integration step size for the MTSM solver set at h = 0.4 in all cases. The column labeled ||error|| shows the norm of the error between analytic (A0, A2) and numerical solution (a0(T), a2(T)). The results show that the MTSM is faster and more accurate than all tested state-of-the-art solvers. The average order used by the MTSM method is ORD = 25 ± 1.
|Solver||Time of calculation [s]||Ratio||||error|||
|Solver||Time of calculation [s]||Ratio||||error|||
3.3 Kepler problem
The Kepler problem describes the motion of a single planet around a heavy sun . The problem is simplified so that the planet does not gravitationally influence the sun and the sun is treated as stationary. The motion of the planet is also limited to the plane (it does not move up or down).
Let y1(t) and y2(t) denote rectangular coordinates centered at the sun, specifying at time t the position of the planet. Also let y3(t) and y4(t) denote a components of velocity in the y1 and y2 directions. If M denotes the mass of the sun, γ the gravitational constant and m mass of the planet, then the attractive force will have the magnitude
The accelerations of the planet in y1 and y2 directions can be calculated as follows:
where the negative sign denotes the inward direction of the acceleration. By adjusting the scales of the variables, the factor γM can be removed. The system of ODEs representing the motion is as follows:
where e is an eccentricity of a rotating body, and . The solution of this system is known and can be either elipses, parabollas or hyperbolas.
The analytical solution of the Kepler problem (30) is defined by the following ellipse
The trajectory for e = 0.75 is in Figure 6 (right). The change in y1 and y2 coordinates in time for different values of e is in Figure 3. The first experiments with solution of the system (30) using the MTSM were performed in .
The term has to be replaced by the equivalent generating ODEs that are added to (30). First, the term r3 is substituted to give
Then the substitution of the term is introduced to give
After the previous substitution, the terms that contain operation division ( and ) can be removed from the system using the following auxiliary ODEs
The new autonomous system of ODEs with the substituted equations without operation division is obtained to give
Note that the new IVP (32) is an autonomous system of ODEs with operation multiplication in the form (8), and such system can be solved using the MTSM with recurrent calculation of the Taylor series’ terms using (9).
the 8 × 8 matrix A has two non-zero elements A(1, 3) = 1, A(2, 4) = 1,
the 2 × 8 matrix B1 has two nonzero elements B1(3, 1) = −1, B1(4, 2) = −1,
the 8 × 4 matrix B2 has four nonzero elements B2(5, 1) = 3, B2(5, 2) = 3, B2(6, 3) = 1, B2(6, 4) = 1,
the 8 × 4 matrix B4 is defined with four nonzero elements B4(7, 1) = −3, B4(7, 2) = −3, B4(8, 3) = −1, B4(8, 4) = −1.
The system of ODEs (32) can be further simplified by reducing the number of multiplications. By substituting (y1y3 + y2y4) we obtain another auxiliary ODE:
Adding the new equation into (32) decreases the number of multiplications to four. We obtain the following:
the 9 × 9 matrix A has two non-zero elements A(1, 3) = 1, A(2, 4) = 1,
the 9 × 9 matrix B1 has eight nonzero elements B1(3, 1) = −1, B1(4, 2) = −1, B1(5, 3) = 3, B1(6, 4) = 1, B1(9, 5) = 1, B1(9, 6) = 1, B1(9, 7) = −1, B1(9, 8) = −1,
the 9 × 2 matrix B2 has two nonzero elements B2(9, 1) = −1, B2(9, 2) = −1,
the 9 × 2 matrix B3 has two nonzero elements B3(7, 1) = −3, B3(8, 2) = −1.
The number of multiplications can be decreased even further by simplifying the remaining ODEs with four multiplications. The multiplication y6y9 can be substituted by an auxiliary ODE:
A new multiplication (y8y9) in the auxiliary ODE can also be substituted, to give
Finally, the bracket can be substituted to give
Adding the three new auxiliary ODEs into (33) decreases the number of multiplications to three
the 12 × 12 matrix A has four non-zero elements A(1, 3) = 1, A(2, 4) = 1, A(5, 10) = 3, A(6, 11) = 1,
the 12 × 8 matrix B1 has eight nonzero elements B1(3, 1) = −1, B1(4, 2) = −1, B1(9, 3) = 1, B1(9, 4) = 1, B1(9, 5) = −1, B1(10, 6) = 1, B1(12, 7) = 2, B1(12, 8) = 2,
the 12 × 9 matrix B2 has nine nonzero elements B2(7, 1) = −3, B2(8, 2) = −1, B2(10, 3) = 1, B2(10, 4) = 1, B2(10, 5) = −1, B2(11, 6) = −1, B2(11, 7) = 1, B2(11, 8) = 1, B2(11, 9) = −1.
The vectors on the right hand side of (8) are defined as follows:
The benchmark results of (32), (33), (34) for the MTSM with fixed integration step size h0.25 = π/25, h0.5 = π/50, h0.75 = π/100 and e = 0.25, e = 0.5, e = 0.75 are in Tables 4, 5 and 6, respectively. The tolerances for all solvers were set to obtain the ||error|| = 10−8, where the ||error|| is computed as the norm of (31). The maximum time of the simulation is set to 2 cycles, T = 4 · π [s]. The obtained ratios indicate faster computation of the MTSM in most cases (ratio > 1).
|e||ode23 ratio||ode45 ratio||ode113 ratio||expTay [s]|
|e||ode23 ratio||ode45 ratio||ode113 ratio||expTay [s]|
|e||ode23 ratio||ode45 ratio||ode113 ratio||expTay [s]|
The autonomous system with low number of multiplications (34) is the fastest one, see Table 6. The greater number of multiplications slows down the calculation considerably. The ratios between the slowest system (32) and the systems of ODEs (33) and (34) for different values of e are in Table 7. The table shows that the system with three term multiplication (34) is approximately two times faster than the system with five term multiplication (32). The difference is caused by the large number of multiplications of higher order terms in the Taylor series. Currently, the Taylor series’ terms are computed directly using (9) without any optimization to decrease the number of multiplications.
Table 8 shows the number of integration steps for different values of e. The number of integration steps grows with increasing values of the constant e. The Kepler problem becomes stiff for a larger constant e. Explicit numerical methods have difficulties with such systems.
|e||number of steps|
The last experiment shows that the solution of the Kepler problem becomes unstable using the state-of-the-art solvers. Results for T = 200 · π [s] and e = 0.75 are in Figure 6. The default error control tolerance Tol = 10−6 was set for all solvers. Only the MTSM provides a stable solution, see Figure 6 (right).
The article dealt with the numerical solution of both linear and nonlinear IVPs represented by a system of ODEs. The results showed that the MTSM overperforms state-of-the-art MATLAB ode solvers in most cases and that the method is more stable in some cases. Future research will focus on further improvements and optimizations of the nonlinear MTSM solver.
This research was financially supported by the Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project “IT4Innovations excellence in science - LQ1602”. The paper also includes the results of the internal BUT FIT project FIT-S-20-6427.
 Kunovský J., Modern Taylor Series Method, Habilitation work, Faculty of Electrical Engineering and Computer Science, Brno University of Technology, 1994Search in Google Scholar
 Veigend P., Nečasová G., Šátek V., High Order Numerical Integration Method and its Applications – The First 36 Years of MTSM, in 2019 IEEE 15th International Scientific Conference on Informatics, 2019, 288–29310.1109/Informatics47936.2019.9119258Search in Google Scholar
 Veigend P., Nečasová G., Šátek V., Model of the telegraph line and its numerical solution, Open Computer Science, 8, in press, 2018, 1–8, https://doi.org/10.1515/comp-2018-000210.1515/comp-2018-0002Search in Google Scholar
 Kocina F., Nečasová G., Veigend P., Šátek V., Kunovský J., Parallel Solution of Higher Order Differential Equations, in 14th International Conference on High Performance Computing & Simulation, Institute of Electrical and Electronics Engineers, 2016, ISBN 978-1-5090-2088-1Search in Google Scholar
 Hairer E., Nørsett S.P., Wanner G., Solving Ordinary Differential Equations I, vol. Nonstiff Problems. Springer-Verlag Berlin Heidelberg, 1987, ISBN 3-540-56670-810.1007/978-3-662-12607-3Search in Google Scholar
 Kocina F., FOS WEBGUI, http://www.fit.vutbr.cz/~iveigend/fos/index_en.phpSearch in Google Scholar
 Kocina F., Kunovský J., Nečasová G., Šátek V., Veigend P., Parallel solution of higher order differential equations, in Proceedings of the 2016 International Conference on High Performance Computing & Simulation (HPCS 2016), Institute of Electrical and Electronics Engineers, 2016, 302–30910.1109/HPCSim.2016.7568350Search in Google Scholar
 Barrio R., Rodríguez M., Abad A., Blesa F., TIDES: A free software based on the Taylor series method, Monografías de la Real Academia de Ciencias de Zaragoza, 35, 2011, 83–95Search in Google Scholar
 Jorba A., Zou M., A software package for the numerical integration of ODE by means of high-order Taylor methods, in Exp. Math., volume 14, 2005, 99–11710.1080/10586458.2005.10128904Search in Google Scholar
 Berz M., COSY INFINITY version 8 reference manual, Technical Report MSUCL–1088, National Superconducting Cyclotron Lab., Michigan State University, East Lansing, Mich., 1997Search in Google Scholar
 Nedialkov N.S., Pryce J., Solving differential algebraic equations by Taylor series III. The DAETS code, JNAIAM J. Numer. Anal. Ind. Appl. Math., 3, 2008, 61–80Search in Google Scholar
 Barrio R., Blesa F., Lara M., VSVO Formulation of the Taylor Method for the Numerical Solution of ODEs, in Computers and Mathematics with Applications, volume 50, 2005, 93–11110.1016/j.camwa.2005.02.010Search in Google Scholar
 Barrio R., Performance of the Taylor series method for ODEs/DAEs, in Applied Mathematics and Computation, volume 163, 2005, 525–545, ISSN 0096300310.1016/j.amc.2004.02.015Search in Google Scholar
 Mohazzabi P., L. Becker J., Numerical Solution of Differential Equations by Direct Taylor Expansion, Journal of Applied Mathematics and Physics, 05(03), 2017, 623–630, 10.4236/jamp.2017.53053Search in Google Scholar
 Baeza A., Boscarino S., Mulet P., Russo G., Zorío D., Approximate Taylor methods for ODEs, Computers & Fluids, 159, 2017, 156 – 166, https://doi.org/10.1016/j.compfluid.2017.10.00110.1016/j.compfluid.2017.10.001Search in Google Scholar
 MATLAB: Choose an ODE Solver, https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.htmlSearch in Google Scholar
 Richard P., Feynman L., Robert B., Sands M., Feynman Lectures on Physics: The New Millennium Edition, Basic Books, 2015Search in Google Scholar
 Chaloupka J., Kocina F., Veigend P., Nečasová G., Šátek V., Kunovský J., Multiple Integral Computations, in 14th International Conference of Numerical Analysis and Applied Mathematics, 1863, American Institute of Physics, 2017, 1–4, 10.1063/1.4992650Search in Google Scholar
 Sehnalová P., Šátek V., 35 Years of Taylor-Kunovsky Simulation Language, in 16th International Conference of Numerical Analysis and Applied Mathematics, volume 2116, American Institute of Physics, 2019, 10.1063/1.5114335Search in Google Scholar
© 2021 Petr Veigend et al., published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.