Abstract
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.
1 Introduction
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) [1].
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 [2]. 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 [3], which focused on linear systems of ODEs and their solutions.
2 The Modern Taylor Series Method
This section of the article briefly introduces the developed method – MTSM. It compiles information from previously published works, mainly [1, 4, 5] and other works cited throughout.
An ODE
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) [6].
The MTSM very effectively implements the variable-step-size, variable-order numerical calculation of differential equations using the Taylor series’ [1]. 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 [1].
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) [7]. Currently, the MTSM has been implemented and tested in MATLAB [8], in C/C++ languages (FOS [9] and TKSL/C software [10]). Additionally, the method can be effectively implemented in hardware [11]. 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 [12] and TAYLOR [13], which includes a detailed description of a variable-step-size version. Other implementations based on Taylor series include ATOMF [14], COSY INFINITY [15], and DAETS [16]. The variable step-size variable-order scheme is also described in [17, 18], and [19], where simulations on a parallel computer are shown. An approach based on an approximate formulation of the Taylor methods can be found in [20].
2.1 Linear Problems
For linear systems of ODEs, Equation (1) is in the form
For a linear system of ODEs (5), the Taylor series (3) can be rewritten in matrix-vector notation as
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 [8] using vectorization. The results are compared with state-of-the-art non-stiff solvers in MATLAB, specifically ode23, ode45 and ode113 [21].
All experiments were performed using MATLAB 2015a on the Salomon supercomputer [22]. 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
Movement of a charged particle in a electromagnetic and electrostatic field is given by Lorenz (electromagnetic) force [1, 23]
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
Position of a particle (upper) and trajectory of the particle (lower).
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.
Results of calculations for movement of a particle.
| Solver | Time of calculation [s] | Ratio | Number of steps |
|---|---|---|---|
| MTSM | 0.017836 | – | 100 |
| ode23 | 193.097 | 10826.2 | 2235540 |
| ode45 | 4.83689 | 2712.2 | 213925 |
| ode113 | 0.78659 | 44.1 | 9446 |
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
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 [24].
For example, the analytic calculation of coefficient A2 using the definite integral
can be transformed into the IVP
The solution of the IVP (21) at the maximum time T represents the calculated value of the definite integral (20) (A0 ≈ a0(T), A2 ≈ a2(T)).
As an example, consider the simple harmonic function
with parameters
can be calculated analytically
as A0 = 1 and
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.
Fourier coefficients of f (t) = sin2(ωt).
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 (
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.
Results of calculations for Fourier coefficients, Eq. (27).
| Solver | Time of calculation [s] | Ratio | ||error|| |
|---|---|---|---|
| MTSM | 0.00135245 | – | 2.7e−15 |
| ode23 | 0.392118 | 289.9 | 7.2e−10 |
| ode45 | 0.0294324 | 21.8 | 8.3e−11 |
| ode113 | 0.0122999 | 9.1 | 1.4e−11 |
Results of calculations for Fourier coefficients, Eq. (28).
| Solver | Time of calculation [s] | Ratio | ||error|| |
|---|---|---|---|
| MTSM | 0.0006809 | – | 1.1e−11 |
| ode23 | 0.533741 | 783.9 | 2.9e−13 |
| ode45 | 0.0456895 | 67.1 | 1.8e−12 |
| ode113 | 0.0159332 | 23.4 | 1.6e−10 |
3.3 Kepler problem
The Kepler problem describes the motion of a single planet around a heavy sun [25]. 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 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 [26].
Kepler problem, solution: e = 0.25 (left), e = 0.5 (middle), e = 0.75 (right).
To solve the system (30) using the MTSM with recurrent calculation of the Taylor series’ terms (9) the system of ODEs has to be transformed to the new autonomous equivalent system (8).
The term
Then the substitution of the term
After the previous substitution, the terms that contain operation division (
The new autonomous system of ODEs with the substituted equations without operation division is obtained to give
The correct initial values for all auxiliary ODEs follow:
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 matrix-vector representation (8) of the system (32) is as follows:
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:
Again, the correct initial values for all auxiliary ODEs are added:
The matrix-vector representation of the ODEs system (33) using the notation of (8) is as follows:
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
Adding the three new auxiliary ODEs into (33) decreases the number of multiplications to three
The correct initial values to all auxiliary ODEs are added as follows:
The matrix-vector representation of the ODEs system (34) in the (8) notation follows:
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).
Results for Kepler problem, Eq. (32).
| e | ode23 ratio | ode45 ratio | ode113 ratio | expTay [s] |
|---|---|---|---|---|
| 0.25 | 146.4 | 4.7 | 1.3 | 0.0293 |
| 0.5 | 122.8 | 4.1 | 1.0 | 0.0536 |
| 0.75 | 90.4 | 3.0 | 0.7 | 0.1158 |
Results for Kepler problem, Eq. (33).
| e | ode23 ratio | ode45 ratio | ode113 ratio | expTay [s] |
|---|---|---|---|---|
| 0.25 | 213.9 | 6.8 | 1.8 | 0.0205 |
| 0.5 | 174.4 | 5.7 | 1.5 | 0.0377 |
| 0.75 | 145.1 | 4.8 | 1.2 | 0.073 |
Results for Kepler problem, Eq. (34).
| e | ode23 ratio | ode45 ratio | ode113 ratio | expTay [s] |
|---|---|---|---|---|
| 0.25 | 312.2 | 12.1 | 2.5 | 0.0144 |
| 0.5 | 262.5 | 9.1 | 2.1 | 0.0263 |
| 0.75 | 221.3 | 7.4 | 1.7 | 0.0509 |
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.
Number of integration steps, Kepler problem (34).
| e | number of steps | |||
|---|---|---|---|---|
| ode23 | ode45 | ode113 | expTay | |
| 0.25 | 105180 | 8916 | 419 | 100 |
| 0.5 | 161323 | 12148 | 654 | 200 |
| 0.75 | 262631 | 19288 | 955 | 400 |
The stiffness indicator σ for the time normalized Kepler problem with different values of constant e is in Figure 5 [27].
The MTSM automatically detects the stiffness in the system and uses larger order (more terms of the Taylor series) in rapidly changing parts of solutions, see Figures 3 and Figure 4.
Kepler problem, MTSM ORD function: e = 0.25 (left), e = 0.5. (middle), e = 0.75 (right).
Kepler problem, stiffness indicator σ: e = 0.25 (left), e = 0.5 (middle), e = 0.75 (right).
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).
Kepler problem, circuit test for the solvers (from the left): ode23, ode45, ode113, MTSM.
4 Conclusion
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.
Acknowledgement
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.
References
[1] Kunovský J., Modern Taylor Series Method, Habilitation work, Faculty of Electrical Engineering and Computer Science, Brno University of Technology, 1994Search in Google Scholar
[2] 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
[3] 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
[4] 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
[5] Veigend P., Šátek V., Nečasová G., Model of the Telegraph line and its Numerical Solution, Open Computer Science, 8(1), 2018, 10–17, 10.1515/comp-2018-0002Search in Google Scholar
[6] 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
[7] Kunovský J., TKSL, https://www.fit.vut.cz/research/product/51/.enSearch in Google Scholar
[8] MATLAB and Simulink software, http://www.mathworks.comSearch in Google Scholar
[9] Kocina F., FOS WEBGUI, http://www.fit.vutbr.cz/~iveigend/fos/index_en.phpSearch in Google Scholar
[10] Šátek V., High Performance Computing Research Group, https://www.fit.vut.cz/research/group/hpc/.enSearch in Google Scholar
[11] 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
[12] 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
[13] 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
[14] Chang Y.F., Corliss G., ATOMF: solving ODEsand DAEs using Taylor series, Computers Math. Applic., 28, 1994, 209–23310.1016/0898-1221(94)00193-6Search in Google Scholar
[15] 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
[16] 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
[17] 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
[18] 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
[19] 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
[20] 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
[21] MATLAB: Choose an ODE Solver, https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.htmlSearch in Google Scholar
[22] IT4I, SOLOMON Supercomputer, https://docs.it4i.cz/salomon/introduction/Search in Google Scholar
[23] Richard P., Feynman L., Robert B., Sands M., Feynman Lectures on Physics: The New Millennium Edition, Basic Books, 2015Search in Google Scholar
[24] 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
[25] Butcher J.C., Numerical methods for ordinary differential equations, John Wiley & Sons, 201610.1002/9781119121534Search in Google Scholar
[26] 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
[27] Söderlind G., Jay L., Calvo M., Stiffness 1952–2012: Sixty years in search of a definition, Bit Numer Math, 55, 2015, 531–558, 10.1007/s10543-014-0503-3Search 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.
