Solving linear and nonlinear problems using Taylor series method

Petr Veigend; Gabriela Nečasová; Václav Šátek

doi:10.1515/comp-2024-0005

Article Open Access

Solving linear and nonlinear problems using Taylor series method

Petr Veigend , Gabriela Nečasová and Václav Šátek

Published/Copyright: July 12, 2024

Published by

Become an author with De Gruyter Brill

Author Information Explore this Subject

From the journal Open Computer Science Volume 14 Issue 1

Abstract

The article deals with the solution of technical initial value problems. To solve such problems, an analytical or numerical approach is possible. The analytical approach can provide an accurate result; however, it is not available for all problems and it is not entirely suitable for calculation on a computer, due to the limited numerical accuracy. For this reason, the numerical approach is preferred. This approach uses ordinary differential equations to approximate the continuous behaviour of the real-world system. There are many known numerical methods for solving such equations, most of them limited in their accuracy, have a limited region of stability and can be quite slow to achieve the acceptable solution. The numerical method proposed in this article is based on the Taylor series and overcomes the biggest challenge, i.e. calculating higher derivatives. The aim of the article is therefore twofold: to introduce the method and show its properties, and to show its behaviour when solving problems composed of linear and nonlinear ordinary differential equations. Linear problems are modelled by partial differential equations and solved in parallel using the PETSc library. The parallel solution is demonstrated using the wave equation, which is transformed into the system of ordinary differential equations using the method of lines. The solution of nonlinear problems is introduced together with several optimisations that significantly increase the calculation speed. The improvements are demonstrated using several numerical examples that are solved using MATLAB software.

Keywords: initial value problems; Taylor series; MTSM; MATLAB; PETSc

1 Introduction

This article deals with the numerical solution of technical initial value problems (IVPs) described by the systems of ordinary differential equations (ODEs). The ODEs are solved using a variable-step variable-order numerical integration method, the modern Taylor series method (MTSM), and the solution is compared with state-of-the-art ODE solvers. The complicated calculation of higher-order derivatives does not need to be calculated because MTSM recurrently calculates the Taylor series terms in each time step [1]. This article compares the numerical results of MTSM with the results obtained by the state-of-the-art Runge-Kutta solvers implemented in MATLAB software and PETSc library and shows the advantages of the MTSM over Runge-Kutta-based solvers. The aim of the article is to show that the MTSM can solve both linear and nonlinear problems faster and more accurately than the state-of-the-art Runge-Kutta-based solvers.

The first implementation of the MTSM was TKSL/386 software (TKSL stands for Taylor-Kunovsky simulation language) [2]. Currently, the MTSM has been implemented and tested in MATLAB [3], in C/C++ languages (FOS [4] and TKSL/C software [5]). Additionally, the method can be effectively implemented in hardware [6]. 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 [7] and TAYLOR [8], which includes a detailed description of a variable-step-size version. Other implementations based on Taylor series include ATOMF [9], COSY INFINITY [10], and DAETS [11]. The variable step-size variable-order scheme is also described in previous studies [12–14], where simulations on a parallel computer are shown. An approach based on an approximate formulation of the Taylor methods can be found in the study by Baeza et al. [15].

The approach based on an approximate formulation of the Taylor methods can be found in the study by Baeza et al. [15]. For example, further research was performed by Amodio et al. [16], which describes the generalised implementation of the Taylor series-based method with its order limited to three.

The MTSM allows for computation with arbitrary accuracy and step size if variable-precision arithmetic and higher-order of method are used. The article by Dimova et al. [17] focuses on the open multi-processing (OpenMP) parallelisation of multiple precision Taylor series method using one computational node. The model problem is the chaotic dynamic system – the classical Lorenz system. The article also briefly mentions the clean numerical simulation (CNS) concept, originally published by Liao [18]. The CNS provides reliable chaotic trajectories in a long enough interval of time.

A hybrid message passing interface (MPI) with OpenMP parallelisation strategy for multiple precision Taylor series method with fixed step size and fixed order is discussed by Hristov et al. [19]. The hybrid strategy was used because OpenMP scalability is slightly better than MPI when using one computational node. The authors claimed that this hybrid strategy can be applied to a large class of chaotic dynamical systems.

The article by Hristov et al. [20] is based on the previous article by Hristov et al. [19] and introduces a modification of CNS with variable step size and fixed order. The order of the Taylor series method is higher than the fixed-order approach, but it results in a reduced number of integration steps thanks to the larger integration step size. Also, the higher-order method increases parallel efficiency.

The article consists of several sections. Section 2 introduces the Taylor-series-based numerical integration method called the modern Taylor series method. Section 3 focuses on a parallel approach of solving linear partial differential equations (PDEs). Section 4 shows the parallel solution of the wave equation, together with performance metrics and detailed numerical results. The formulation of the method for nonlinear problems is introduced in Section 5, including several optimisations. Numerical results for selected examples of nonlinear problems are summarised in Section 6.

2 Higher-order Taylor series method

This section introduces the developed method – MTSM. More information can be found in previous studies [1,21–23]. An ODE with an initial condition

(1) y ′ ( t ) = f ( t , y ( t ) ) , y ( 0 ) = y 0 ,

is called an initial value problem. The best-known and the most accurate method of the numerical solution of (1) is to construct the Taylor series in the form

(2) y i + 1 = y i + h ⋅ f ( t i , y i ) + h 2 2 ! ⋅ f ′ ( t i , y i ) + … + h n n ! ⋅ f [ n − 1 ] ( t i , y i ) ,

where h is the integration step size, y i ≈ y ( t i ) is the approximation of the current value, and y i + 1 ≈ y ( t i + h ) is the approximation of the next value of the function y ( t ) [24].

MTSM very effectively implements the variable-step-size, variable-order numerical solution of IVPs 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. Equation (2) can then be rewritten in the form

(3) y i + 1 = p ( 0 ) + p ( 1 ) + p ( 2 ) + … + p ( n ) ,

where p ( j ) , j = 0 … n denotes the Taylor series terms. The MTSM transforms the input problem into a system of autonomous ODEs, which allows the recurrent calculation of the Taylor series terms.

The Taylor series terms in each step are truncated when the stopping rule for the last σ number of Taylor series terms is met:

(4) ∑ j = n − σ n ∣ ∣ p ( j ) ∣ ∣ ≤ ε ,

where ε is the required accuracy of the calculation, which is chosen for each example. In this article, we consider σ = 3 .

An 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 P denote the function which changes during the computation and defines the number of Taylor series terms used in the current integration step ( P i + 1 = n ). More information about the method and additional comparison with state-of-the-art methods can be found in the study by Veigend et al. [21].

3 Linear problems and parallel solution

For linear systems of ODEs ( y ′ ( t ) = A y ( t ) + b ) , and (2) can be rewritten in matrix-vector notation as

(5) y i + 1 = y i + h ( A y i + b ) + h 2 2 ! A ( A y i + b ) + … + h n n ! A ( n − 1 ) ( A y i + b ) ,

where A is the constant Jacobian matrix and b is the constant right-hand side. Moreover, Taylor series terms in (3) can be computed recurrently using

(6) p ( 0 ) = y i , p ( 1 ) = h ( A y i + b ) , p ( j ) = h j A p ( j − 1 ) , j = 2 , … , n .

When solving linear problems in parallel, (5) can be rewritten as follows:

(7) y i + 1 = A y y i + A b b ,

where the matrices A y and A b are defined as follows:

(8) A y = ∑ k = 0 n h k k ! A k , A b = ∑ k = 1 n h k k ! A k − 1 .

Let A D l denote the submatrix of the matrix A decomposed by rows, where l = { 1 , 2 , … , n P } and n P is the number of processes. Let A be a matrix of size m × m . The number of rows of matrix A is evenly divided among processes. Matrices A y and A b are constant, and matrices are precalculated only once at the beginning of the calculation:

(9) A ˆ y D l = ∑ k = 0 n h k k ! A D l k , A ˆ b D l = ∑ k = 1 n h k k ! A D l k − 1 ,

where l = { 1 , 2 , … , n P } . After the parallel precalculation, the matrix A ˆ y G is obtained by gathering individual matrices A ˆ y D l . Similarly, the matrix A ˆ b G :

(10) A ˆ y G = ( A ˆ y D 1 , A ˆ y D 2 , … , A ˆ y D p ) T , A ˆ b G = ( A ˆ b D 1 , A ˆ b D 2 , … , A ˆ b D p ) T .

Final matrix A ˆ and vector b ˆ are calculated afterwards, and I is the identity matrix

(11) A ˆ = A ˆ y G A + I , b ˆ = A ˆ b G b .

Using (11), we can rewrite (7) and solve it in parallel using the row-wise decomposition of matrix A ˆ

(12) y i + 1 = A ˆ y i + b ˆ .

4 Linear problem example

The analytic notation of the wave equation is the second-order hyperbolic PDE [25] follows:

(13) ∂ 2 u ( x , t ) ∂ x 2 = ∂ 2 u ( x , t ) ∂ t 2 ( x , t ) ∈ ( 0 , L ) × ⟨ 0 , T ⟩ .

The homogeneous Dirichlet boundary conditions are set:

(14) u ( 0 , t ) = u ( L , t ) = 0 , 0 ≤ t ≤ T ,

where L is the string length and T is the maximum simulation time. The initial values follow:

(15) u ( x , 0 ) = sin ( π x ) ,

(16) ∂ u ( x , 0 ) ∂ t = 0 , 0 < x < L .

The wave equation describes the oscillations of an ideal string of a specified length. Both ends of the string are fixed in time (see the boundary conditions (14)). The initial velocity of the string is zero (16). The initial position of the string is modelled as a sine function (15). The resulting system of ODEs y ′ ( t ) = A y ( t ) + b , with initial conditions y ( 0 ) = y 0 arising from the Method Of Lines (MOL) is in the form:

(17) A = A 11 A 12 A 21 A 22 , y 0 = u ( x , 0 ) ∂ u ( x , 0 ) ∂ t , b = 0 ,

where A 12 is the spatial discretisation matrix, the three-point central difference formula (coefficients 1 , − 2 , 1 ) was used. The A 21 is the identity matrix, and matrices A 11 and A 22 are zero matrices. The size of the problem is denoted as S . The sparsity patterns of the matrices A and A ˆ precalculated using (12) for S = 100 are in Figure 1(a) and (b), respectively.

$Figure 1 Sparsity patterns of input matrices, wave equation, three-point central difference formula, S = 100 S=100 . (a) Matrix A {\bf{A}} and (b) matrix A ˆ \hat{{\bf{A}}} .$

Figure 1

Sparsity patterns of input matrices, wave equation, three-point central difference formula, S = 100 . (a) Matrix A and (b) matrix A ˆ .

The sparsity of the matrices A and A ˆ increases with the problem size. The sparsity of A ˆ decreases during precalculation as seen in Figure 1. Let us denote the size of the submatrix as m = S − 1 . The following text will elaborate number of nonzero elements (nnz) of matrices A and A ˆ in more detail. Total number of nonzero elements in matrix A :

(18) n n z A = n n z A 12 + m ,

(19) n n z A 12 = 3 ( m − 2 ) + 4 .

First, let us define the number of nonzero elements for matrix A 12 . Note that the first and last rows contain two elements (four overall), all other rows contain three elements. The last term, m , reflects the number of nonzero elements of the identity matrix A 21 . The term P max denotes the maximum order of MTSM_PRECALC, that is the number of Taylor series terms (Table 2). The total number of nonzero elements in matrix A ˆ for a given maximum order P max can be calculated as follows:

(20) n n z A ˆ = n n z A + n n z A ˆ 1 + n n z A ˆ 2 .

Table 2

Parameters for linear simulation experiments

Parameter	Value
Problem size	S = 256,000
Number of first-order ODEs	2 S ‒ 2
Length of the string	L = 25,600 ( mm )
Spatial step size	Δ x = L ∕ S = 0.1 ( mm )
Accuracy	rel TOL = abs TOL = 1 × 1 0 ‒ 10
Integration step size	h = 0.4 ( s )
Maximum simulation time	T = 10,000 ⋅ h ( s )
Maximum order of MTSM	P max = 64
Maximum order of MTSM_PRECALC	P max = 25

The expression n n z A ˆ 1 denotes the number of nonzero elements for odd orders of the MTSM_PRECALC:

(21) n n z A ˆ 1 = ∑ i = 1 P 1 2 ( m − i ) + 2 ( m − i − 1 ) ,

where 2 ( m − i ) denotes the number of nonzero elements in the submatrix A 21 and 2 ( m − i − 1 ) in the submatrix A 12 . If P max is odd, then P 1 = ( P max − 1 ) ∕ 2 , if P max is even, then P 1 = P max ∕ 2 − 1 .

The expression n n z A ˆ 2 denotes the number of nonzero elements for even orders of the MTSM_PRECALC:

(22) n n z A ˆ 2 = 2 n n z A 12 + ∑ i = 1 P 2 4 ( m − i − 1 ) ,

where 2 n n z A 12 denotes the number of nonzero elements for P = 2 and 4 ( m − i − 1 ) is number of nonzero elements for even P > 2 . If P max is odd, then P 2 = ( P max − 1 ) ∕ 2 , if P max is even, then P 2 = P max ∕ 2 .

The density ( d ) of a matrix (percentage of nonzero elements) can be expressed as follows:

(23) d [ % ] = n n z ( 2 S − 2 ) 2 ⋅ 100 ,

where nnz is number of nonzero elements in matrix A or A ˆ and ( 2 S − 2 ) 2 is the total number of elements. The density of the matrix A for a given problem size S = 256,000 is 3.91 × 1 0 − 4 % . The matrix A ˆ precalculated to the order 25 is approximately 25.5 times less sparse, so the density is 9.96 × 1 0 − 3 % .

The numerical solution of a given PDE is obtained by the MOL, which discretises the problem in spatial dimension and integrates the semi-discrete system of ODEs. The resulting system of ODEs is solved as IVP in the time domain using numerical integration methods. Table 1 shows the PETSc solvers used for the numerical experiments. The results obtained by Taylor series-based solvers were compared with two selected Runge-Kutta solvers. The TSRK5DP implements the fifth-order Dormand-Prince 5(4) method with a fourth-order embedded method (known as ode45 in MATLAB). The TSRK8VR is the eighth-order robust Verner scheme with a seventh-order embedded method with thirteen stages. All Taylor series-based solvers were implemented using PETSc library routines.

Table 1

PETSc numerical integration methods

Solver	Method
TSRK5DP	Dormand-Prince [26] (equivalent with ode45 in MATLAB)
TSRK8VR	Verner Runge-Kutta [27]
MTSM	Linear MTSM (5)
MTSM_PRECALC	MTSM with the A ˆ precalculation (11) and (12)

The simulation experiments were performed on the Barbora supercomputer cluster, IT4Innovations National Supercomputing Center, Ostrava, Czech Republic [28]. The 32 compute nodes were utilised, each running 36 processes, resulting in a total count of 1152 MPI processes. The PETSc library [29–32] leverages the MPI standard for message-passing communication and offers data structures along with numerical methods that incorporate automatic step size control. It is explicitly designed for the parallel solution of scientific applications modelled by PDEs. All Taylor series-based solvers were implemented using PETSc library routines. Walltime was set to 15 min for all experiments. The parameters for the simulation experiments are shown in Table 2. Note that the maximum order for MTSM is 64, and the maximum order for MTSM_PRECALC is 25 (matrix A ˆ was precalculated to the order of 25). Three-point central difference formula is used to discretise the spatial domain. The integration step size was chosen based on the numerical stability analysis of Taylor-series-based solvers [33,34].

4.1 Performance metrics

The performance metrics evaluate, characterise, diagnose, and tune parallel performance [35]. Numerical results are presented using the following performance metrics: average time, speedup, speedup against the TSRK5DP solver, parallel efficiency, and parallel cost metrics.

The average computation time ( t a ) is defined as follows:

(24) t a [ s ] = ∑ i = 1 n R t n R ,

where n R denotes the total number of runs.

The speedup ( s ) of the solver is defined as follows:

(25) s = t 1 t N ,

where t 1 denotes the serial computation time of one compute node and t N denotes the parallel computation time of n N compute nodes. When s > 1 , the speedup metric conveys performance improvement. On the contrary, when s < 1 , the speedup metric conveys performance degradation. The ideal speedup (ideal scaling) is defined as s = n N . The superlinear speedup is achieved when s > n N .

The speedup-against ( s a ) ratio of the solver is defined as follows:

(26) s a = t a 1 t a 2 ,

where t a 1 denotes the average time of computation of the reference solver TSRK5DP against t a 2 which denotes the calculation time using one of the solvers in Table 1.

s a ≫ 1 indicates a significantly faster computation time using the given solver.

The parallel efficiency ( e N ) is defined as follows:

(27) e N [ % ] = s n N ⋅ 100 = t 1 t N n N ⋅ 100 ,

The parallel cost ( c N ) of an algorithm is defined as follows:

(28) c N [ node-hours ] = t a 3,600 ⋅ n N ,

where t a ∕ 3,600 is the average time in hours.

The parallel cost ratio ( c R ) is defined as follows:

(29) c R = c N 1 c N 2 ,

where c N 1 is the parallel cost for one compute node and c N 2 is the parallel cost for more than one compute node. Parallel cost is calculated using (28).

The parallel speedup-cost ratio ( s R ) is defined as follows:

(30) s R = s c R ,

where s is calculated by (25) and c R by (29).

4.2 Numerical results

Selected problem sizes (denoted as S ) are 64,000, 128,000, 256,000, 512,000, and 1,024,000 second-order ODEs. In this section, results for S = 256,000 will be introduced in detail. Since the wave equation is a second-order PDE, it is necessary to reduce the order of derivatives to obtain the system of the first-order ODEs. Therefore, the number of first-order ODEs in (17) is two times bigger, that is, 2 S − 2 (Table 2). For simplicity, n e = 2 S denotes the problem size including two Dirichlet boundary conditions (14).

In the final part of this section, results for other problem sizes will be discussed. The average order of MTSM is 19. The density of matrices A and A ˆ (23) is 3.91 × 1 0 − 4 % and 9.96 × 1 0 − 3 % , respectively.

The number of integration steps and the average step sizes for each solver are shown in Table 3. The MTSM uses a step size approximately 7.8 times larger than the TSRK5DP solver and approximately 3.2 times larger than the TSRK8VR solver.

Table 3

Number of integration steps, average step sizes

Solver	# steps	Average h
MTSM_PRECALC	10,000	4.00 × 1 0 − 1
MTSM	10,000	4.00 × 1 0 − 1
TSRK5DP	78,238	5.11 × 1 0 − 2
TSRK8VR	32,000	1.26 × 1 0 − 1

Figure 2 visualize the performance metrics introduced in Section 4.1 for all selected solvers. The results for MTSM_PRECALC are marked in red, for MTSM in blue, for TSRK5DP in green, and for TSRK8VR in black. Figure 2(a) shows the average computation time (24) where the dashed lines show the ideal average times for the given number of processes. Figure 2(b) shows the parallel efficiency (27). Figure 2(c) depicts the speedup for each solver (25), and the dashed line shows the ideal speedup for the given number of processes. Figure 2(d) shows the speedup ratio with respect to the TSRK5DP solver (26), and s ≫ 1 indicates significantly faster computation using the given solver defined in Table 1. We can clearly see that the MTSM_PRECALC solver overperforms the state-of-the-art solvers. Only the MTSM and MTSM_PRECALC solvers can provide results for all 1–32 compute nodes.

$Figure 2 Performance metrics: average time, parallel efficiency, parallel speedup, and speedup against the TSRK5DP solver, n e = 512,000 ne=\hspace{0.1em}\text{512,000}\hspace{0.1em} . (a) Average time and (b) parallel efficiency, (c) parallel speedup, and (d) speedup against the TSRK5DP solver.$

Figure 2

Performance metrics: average time, parallel efficiency, parallel speedup, and speedup against the TSRK5DP solver, n e = 512,000 . (a) Average time and (b) parallel efficiency, (c) parallel speedup, and (d) speedup against the TSRK5DP solver.

Figure 3 shows the parallel cost metrics for each solver. Each subfigure contains three curves. The curve labelled with the solver’s name shows the parallel speedup of a given solver (25). The magenta curve represents the parallel cost ratio (29), and the cyan curve shows the parallel speedup-cost tradeoff (30). Figure 3(a) shows that the ideal number of compute nodes for the MTSM_PRECALC solver is 32 (1152 MPI processes) because the speedup (red curve) increases proportionally with the parallel cost ratio (magenta curve). Therefore, the speedup-cost ratio (cyan curve) also increases. MTSM, TSRK5DP, and TSRK8VR solvers are ideal to use 16 nodes (576 MPI processes), see Figures 3(b), (c), and (d), respectively.

$Figure 3 Parallel cost metrics: speedup, parallel cost and parallel speedup-cost ratio for each solver, n e = 512,000 ne=\hspace{0.1em}\text{512,000}\hspace{0.1em} . (a) Parallel speedup-cost tradeoff, MTSM_PRECALC, (b) parallel speedup-cost tradeoff, MTSM, (c) parallel speedup-cost tradeoff, TSRK5DP, and (d) parallel speedup-cost tradeoff, TSRK8VR.$

Figure 3

Parallel cost metrics: speedup, parallel cost and parallel speedup-cost ratio for each solver, n e = 512,000 . (a) Parallel speedup-cost tradeoff, MTSM_PRECALC, (b) parallel speedup-cost tradeoff, MTSM, (c) parallel speedup-cost tradeoff, TSRK5DP, and (d) parallel speedup-cost tradeoff, TSRK8VR.

Tables 4, 5, 6, and 7 summarise results for wave equation discretised in the spatial domain using the three-point central difference formula for each problem size n e , namely, 128,000, 256,000, 512,000, 1,024,000, and 2,048,000 ODEs. These results are averages of values for 1–32 compute nodes.

Table 4

Average time, n e = 512,000

Solver	# compute nodes ( t a ) [s]
	1	4	8	12	16	20	24	28	32
MTSM_PRECALC	16.78	3.49	1.67	0.88	0.78	0.77	0.62	0.52	0.53
MTSM	40.96	25.58	16.19	12.82	11.05	10.75	9.83	10.13	9.19
TSRK5DP	141.11	52.95	29.36	21.63	17.36	15.74	13.64	—	—
TSRK8VR	189.06	59.53	32.87	23.51	18.81	15.98	—	—	—

Table 5

Yes/No table

Solver	Problem size ( n e )
	128,000	256,000	512,000	1,024,000	2,048,000
MTSM_PRECALC	Yes	Yes	Yes	Yes	Yes
MTSM	Yes	Yes	Yes	Yes	No (6)
TSRK5DP	Yes	Yes	No (27)	No (4)	—
TSRK8VR	Yes	Yes	No (21)	No (2)	—

Table 6

Average efficiency (e_N) comparison for 1–32 nodes

Solver	Problem size ( n e )
	128,000	256,000	512,000	1,024,000	2,048,000
MTSM_PRECALC	51.10 ̲	77.41 ̲	127.64 ̲	136.43 ̲	119.31 ̲
MTSM	15.30	19.77	25.62	33.67	—
TSRK5DP	25.91	29.86	—	—	—
TSRK8VR	29.76	40.24	—	—	—

The underlined values indicate much faster calculation than the state-of-the-art numerical solvers.

Table 7

Average speedup against TSRK5DP comparison for 1–32 nodes

Solver	Problem size ( n e )
	128,000	256,000
MTSM_PRECALC	16.04 ̲	20.90 ̲
MTSM	1.04	1.63
TSRK8VR	0.97	1.23

The underlined values indicate much faster calculation than the state-of-the-art numerical solvers.

Table 4 shows the average computation times for a selected number of processes, that is, for 1, 4, 8, 12, 16, 20, 24, 28, and 32 compute nodes. The MTSM_PRECALC is the fastest of all solvers.

Table 5 indicates whether a solver calculates (“Yes”) or does not calculate (“No”) the result for 1–32 nodes for a given problem size. The notation “No ( X )” implies that a given solver calculated the result for the maximum number of compute nodes denoted as X .

Notice that for the problem size greater than n e = 256,000 , the TSRK5D and TSRK8VR solvers did not calculate the result for all 1–32 compute nodes because the maximum walltime (15 min) was exceeded. For example, for n e = 1,024,000 , the TSRK5DP and TSRK8VR did not compute the results for all 32 compute nodes and were able to use four and two nodes maximally, respectively. Similarly, the MTSM did not calculate results for problem sizes greater than n e = 1,024,000 . Only the MTSM_PRECALC solver calculated results for all problem sizes for all 1–32 nodes.

The cells in Table 6 with average parallel efficiency (27) greater than or equal to 50% are in bold. The MTSM_PRECALC solver offers an efficiency greater than 50% for all problem sizes.

The cells in Table 7 showing the speedup ratio (26) with respect to the TSRK5DP solver where s a ≫ 1 is also marked in bold. The MTSM_PRECALC solver is always faster than the TSRK5DP solver.

Runge-Kutta solvers did not provide results for n e > 256,000 ODEs for all 1–32 compute nodes (Table 5), for that reason, Table 7 has only two columns.

5 Nonlinear problems and optimisations

For nonlinear problems, IVP (1) has the following form:

(31) y ′ = A y + B 1 y j k + B 2 y j k l + … + b , y ( 0 ) = y 0 ,

where A ∈ R n e × n e is the constant matrix for the linear part of the system (Section 3), matrices B 1 ∈ R n e × n m j k , B 2 ∈ R n e × n m j k l are the constant matrices for nonlinear part of the system. The vector b ∈ R n e is the right-hand side for the forces incoming to the system, y 0 is a vector of the initial conditions, and symbol n e stands for the number of equations of the system of ODEs. Symbols n m j k and n m j k l represent the number of two and three-functions multiplications, respectively.

The unknown function y j k ∈ R n m j k represents the vector of two-terms multiplications y j ⊙ y k and similarly y j k l ∈ R n m j k l represents the vector of three-terms multiplications y j j ⊙ y k k ⊙ y l l , where indices j , k , j j , k k , l l ∈ ( 1 , … , n e ) come from multiplications terms in (31). The operation ⊙ stands for element-by-element multiplication, i.e. y j ⊙ y k is a vector ( y j 1 y k 1 , y j 2 y k 2 , … , y j n m j k y k n m j k ) T . For simplification, the matrices A , B 1 , B 2 , … and the vector b are constant. The higher derivatives of the terms B 1 y j k , B 2 y j k l in (31) can be included in a recurrent calculation of the Taylor series terms p B 1 and p B 2

(32) p ( 1 ) A = h ( A y i + b ) , p ( 1 ) B 1 = h ( B 1 y j k ) , p ( 1 ) B 2 = h ( B 2 y j k l ) , p ( r ) A = h r A p ( r − 1 ) , p ( r ) B 1 = h r B 1 ∑ a = 1 r p ( a − 1 ) j ⊙ p ( r − a ) k , p ( r ) B 2 = h r B 2 ∑ a = 0 r − 1 p ( a ) j j ⊙ ∑ b = 1 r − a p ( b − 1 ) k k ⊙ p ( r − a − b ) l l ,

where r = 2 , … , n . Finally, the Taylor series terms are calculated as a sum of linear and nonlinear terms

(33) p ( s ) = p ( s ) A + p ( s ) B 1 + p ( s ) B 2 , s = 1 , … , n ,

where r and s are the current indexes of the Taylor series terms, a and b are the auxiliary indexes for the summation of two and three-terms multiplications in the nonlinear part of the Taylor series, and p ( s ) A is the linear term computed using the recurrent calculation for linear systems. The next value of the function can be calculated using

(34) y i + 1 = p ( 0 ) i + p ( 1 ) i + p ( 2 ) i + … + p ( n i ) ,

where p ( 0 ) i is the value of the function y i , p ( 1 ) i , …, p ( n ) i are the Taylor series terms calculated using (33). Multiplication terms of the Taylor series for more multiplications p B 3 , p B 4 , … can be calculated recurrently in a similar way. More information can be found in the study by Veigend et al. [21].

The important fact to consider is the number of element-by-element multiplications to calculate the solution. The number of element-by-element multiplications used in (32) is visualised in Figure 4.

Figure 4

Number of multiplications.

Due to the fact that the number of element-by-element multiplications increases rapidly when more functions are multiplied, the calculation of higher derivatives has to be optimised. Several optimisations are discussed in this article.

5.1 Auxiliary generating equations

This optimisation is used automatically when using the presented method. However, the equations that contain terms with many function multiplications (the example behaviour for five multiplication is in Figure 5(b)) can be simplified by applying the generating process again. The newly generated equation contains several terms and therefore decreases the number of function multiplications in the original term.

5.2 Partly calculated results of the previous term

This optimisation is universal and can be used for any system. Consider Table 8, which contains the Taylor series terms used when calculating higher derivatives. The coloured areas when calculating the Taylor series terms can be saved and used when calculating further terms.

Table 8

Three function multiplications

The optimisation becomes more pronounced for more function multiplications, and it can be performed multiple times. This is shown in Table 9, which shows partial results being saved and used. Both partial (darker colours) and full results (lighter colours) are saved and used.

Table 9

Five function multiplications – partial and full results

The number of performed function multiplications with and without optimisations is shown in Figure 5(a) and (b).

$Figure 5 Comparison of a number of element-by-element multiplications ( ⊙ \odot ) with and without optimisations. (a) Three function multiplications and (b) Five function multiplications.$

Figure 5

Comparison of a number of element-by-element multiplications ( ⊙ ) with and without optimisations. (a) Three function multiplications and (b) Five function multiplications.

The figures shows that the number of operations decreases quite dramatically compared to the case without any optimisations (Figure 4), and it is almost linear.

5.3 Variable step size

When solving nonlinear problems, the P function can oscillate quite rapidly. For some nonlinear problems, however, it rapidly changes near the beginning of the calculation and then stabilises. When this happens, the method often calculates with a relatively small integration step. Due to the fact that the method can automatically adjust the value of P in the current step based on the size of the step, the size can be dynamically increased to decrease the overall number of operations the method has to perform.

To work with this optimisation, the constant h scale is defined as the scaling factor for the size of integration step h

h new = h scale ⋅ h .

The new value for the size of the integration step h new is used until the calculation ends. The scaling factor is only applied when

∑ a = i − 3 i P ( a ) ≤ P min ⋅ 3 , i = 4 … n T ,

where n T is the number of time steps and P min is the value of the P function that has to be kept for three integration steps. This approach is useful for problems where the previous approach cannot be used (i.e. for systems that only contain two function multiplications). These optimisations and their advantages and disadvantages are discussed on a set of benchmarks in Section 6.

5.4 Further optimisations

Additional optimisations are possible and are being actively studied. One of the approaches that is currently being tested limits the number of multiplications to two by expressing the Taylor series terms using systems of differential-algebraic equations (DAEs). The results are promising and are going to be published at a later date.

6 Nonlinear problem examples

The previously mentioned optimisations are tested using the selected set of nonlinear nonstiff problems presented in [36]. The selected benchmark problems are a subset of a full benchmark collection, which is intended to show the potential strengths and weaknesses of newly developed numerical methods.

The selected problems are going to be analysed and then solved using the nonlinear MTSM solver with or without the optimisations.

Note that the systems of ODEs have to be always transformed into the autonomous system of homogeneous ODEs without division and with just basic arithmetic operations [21].

Experiments were repeated 100 times. The median computation time ( t m ) was calculated from the computation times during the individual runs of the solvers. The speedup-against ( s a ) ratio of the solver is defined as follows:

(35) s a = t m 1 t m 2 ,

where t m 1 denotes the median calculation time using one of the MATLAB ode solvers from Table 10 against t m 2 which denotes the median calculation time using of the MTSM solvers listed in the same table.

Table 10

MATLAB numerical methods

Solver	Method
ode23	Bogacki-Shampine [37]
ode45	Dormand-Prince [26]
ode113	Adams-Bashfort method with predictor-corrector PECE scheme
MTSM orig	Nonlinear MTSM solver without optimisations
MTSM opt	Nonlinear MTSM solver with optimisations from Section 5

The s a ≫ 1 indicates a significantly faster computation time using the MTSM solvers.

The error (err) is calculated as norm of the difference between numerical solution MTSM solvers and MATLAB ode solvers at the time T

err = ‖ y M ( T ) − y O ( T ) ‖ ,

where y M ( T ) is the value calculated by the MTSM solver at the time T and y O ( T ) is the value calculated by the ode solver at the time T . Table 10 lists the state-of-the-art solvers that are part of MATLAB and implemented MTSM solvers for MATLAB used in this section.

The parameters for all used solvers are summarised in Table 11.

Table 11

Parameters for nonlinear simulation experiments

Parameter	Value
MTSM accuracy	ε = 1 × 1 0 ‒ 9
Relative, absolute tolerances	rel TOL = abs TOL = 1 × 1 0 ‒ 9
Maximum simulation time	T = 20 (s)
Maximum order for MTSM solvers	P max = 64

6.1 Problem A2

Problem A2 is defined as follows:

y ′ = − y 3 2 = − 0.5 y 3 y ( 0 ) = 1 .

To solve this nonlinear system using MTSM, the term y 3 has to be replaced by a set of auxiliary ODEs

y 1 ′ = − 0.5 y 1 3 y 1 ( 0 ) = 1 y 2 = y 1 3 y 2 ′ = 3 y 1 2 y 1 ′ = 3 y 1 2 ( − 0.5 y 1 3 ) = − 1.5 y 1 2 y 2 y 2 ( 0 ) = y 1 ( 0 ) 3 y 3 = y 1 2 y 3 ′ = 2 y 1 y 1 ′ = 2 y 1 ( − 0.5 y 1 3 ) = − y 1 y 2 y 3 ( 0 ) = y 1 ( 0 ) 2 ,

so the final system becomes

y 1 ′ = − 0.5 y 2 y 1 ( 0 ) = 1 y 2 ′ = − 1.5 y 2 y 3 y 2 ( 0 ) = y 1 ( 0 ) 3 y 3 ′ = − y 1 y 2 y 3 ( 0 ) = y 1 ( 0 ) 2 ,

which can be transformed into a matrix-vector representation (31)

A = 0 − 0.5 0 0 0 0 0 0 0 B 1 = 0 0 − 1.5 0 0 − 1 y j k = y 2 y 3 y 1 y 2 .

The numerical results are presented in Table 12. MTSM solvers outperform the state-of-the-art solvers.

Table 12

Results for problem A2, h = 0.5 a , h scale = 5 , P min = 9

			MTSM orig		MTSM opt
Solver	# of steps	Time of calculation [s]	err	Ratio ( s a )	e r r	Ratio ( s a )
MTSM orig	40	0.00152015	—	—	—	—
MTSM opt	16	0.0007533	—	—	—	—
ode23	28,550	0.168798	1.67406 × 1 0 ‒ 7	114.25	2.19986 × 1 0 ‒ 7	230.1
ode45	2,177	0.0040288	1.67407 × 1 0 ‒ 7	2.73	2.19987 × 1 0 ‒ 7	5.49
ode113	232	0.0031412	1.67406 × 1 0 ‒ 7	2.13	2.19986 × 1 0 ‒ 7	4.28

Optimisations in Sections 5.1 and 5.3 are used. The automatic change in step size at t = 5 s is shown in Figure 6.

$Figure 6 Order function P P for problem A2, h = 0.5 s h=0.5\hspace{0.33em}{\rm{s}} , h scale = 5 {h}_{{\rm{scale}}}=5 , P min = 9 {P}_{\min }=9 .$

Figure 6

Order function P for problem A2, h = 0.5 s , h scale = 5 , P min = 9 .

6.2 Problem B4

This problem is interesting because it contains nontrivial mathematical operations that have to be replaced using auxiliary equations. The problem is defined as follows:

y 1 ′ = − y 2 − y 1 y 3 y 1 2 + y 2 2 y 1 ( 0 ) = 3 y 2 ′ = y 1 − y 2 y 3 y 1 2 + y 2 2 y 2 ( 0 ) = 0 y 3 ′ = y 1 y 1 2 + y 2 2 y 3 ( 0 ) = 0 .

After generating auxiliary equations and decreasing the number of function multiplications (optimisation from Section 5.1), the system of ODEs can be written as follows:

(36) y 1 ′ = − y 2 − y 11 y 2 ′ = y 1 − y 2 y 3 y 5 y 3 ′ = y 1 y 5 y 4 ′ = − 2 y 3 y 8 y 9 − 2 y 3 y 8 y 10 y 5 ′ = 2 y 3 y 6 y 9 + 2 y 3 y 6 y 10 y 6 ′ = 8 y 3 y 6 y 7 y 9 + 8 y 3 y 6 y 7 y 10 y 7 ′ = 6 y 3 y 6 y 8 y 9 + 6 y 3 y 6 y 8 y 10 y 8 ′ = 4 y 3 y 5 y 6 y 9 + 4 y 3 y 5 y 6 y 10 y 9 ′ = − 2 y 9 − 2 y 9 y 11 y 10 ′ = 2 y 1 y 2 − 2 y 3 y 5 y 10 y 11 ′ = − y 2 y 3 y 5 − y 3 y 5 y 11 + y 8 y 9 + 2 y 1 y 6 y 9 y 12 + 2 y 1 y 6 y 10 y 12 y 12 ′ = 2 y 11

with initial conditions y ( 0 ) = ( 3 , 0 , 0 , y 1 ( 0 ) 2 + y 2 ( 0 ) 2 , 1 y 4 ( 0 ) , y 5 ( 0 ) 4 , y 5 ( 0 ) 3 , y 5 ( 0 ) 2 , y 1 ( 0 ) 2 , y 2 ( 0 ) 2 , y 1 ( 0 ) y 3 ( 0 ) y 5 ( 0 ) , y 3 ( 0 ) 2 ) T The matrix-vector representation of system of ODEs (36) using the notation from (32) follows:

the 12 × 12 matrix A has five nonzero elements A ( 1 , 2 ) = − 1 , A ( 1 , 11 ) = − 1 , A ( 2 , 1 ) = 1 , A ( 9 , 9 ) = − 2 and A ( 12 , 11 ) = 2 ,
the 12 × 4 matrix B 1 has four nonzero elements B 1 ( 3 , 1 ) = 1 , B 1 ( 9 , 2 ) = − 2 , B 1 ( 10 , 3 ) = 2 , B 1 ( 11 , 4 ) = 1 ,
the 12 × 8 matrix B 2 has eight nonzero elements B 2 ( 2 , 1 ) = − 1 , B 2 ( 4 , 2 ) = − 2 , B 2 ( 4 , 3 ) = − 2 , B 2 ( 5 , 4 ) = 2 , B 2 ( 6 , 5 ) = 2 , B 2 ( 10 , 6 ) = − 2 , B 2 ( 11 , 7 ) = − 1 and B 2 ( 11 , 8 ) = − 1 ,
the 12 × 8 matrix B 3 has eight nonzero elements B 3 ( 6 , 1 ) = 8 , B 3 ( 6 , 2 ) = 8 , B 3 ( 7 , 3 ) = 6 , B 3 ( 7 , 4 ) = 6 , B 3 ( 8 , 5 ) = 4 , B 3 ( 8 , 6 ) = 4 , B 3 ( 11 , 7 ) = 2 and B 3 ( 11 , 8 ) = 2 ,
two-term multiplications y j k = ( y 1 y 5 , y 9 y 11 , y 1 y 2 , y 8 y 9 ) T ,
three-term multiplications y j k l = ( y 2 y 3 y 5 , y 3 y 8 y 9 , y 3 y 8 y 10 , y 3 y 6 y 9 , y 3 y 6 y 10 , y 3 y 5 y 10 , y 2 y 3 y 5 , y 3 y 5 y 11 ) T ,
four-term multiplications y j k l m = ( y 3 y 6 y 7 y 9 , y 3 y 6 y 7 y 10 , y 3 y 6 y 8 y 9 , y 3 y 6 y 8 y 10 , y 3 y 5 y 6 y 9 , y 3 y 5 y 6 y 10 , y 1 y 6 y 9 y 12 , y 1 y 6 y 10 y 12 ) T ,
vector representing the right-hand side of the system b = 0 .

The numerical results are in Table 13. Note that optimisations 5.2 and 5.3 are also used.

Table 13

Results for problem B4, h = 0.5 s , h scale = 2.5 , and P min = 12

			MTSM orig		MTSM opt
Solver	# of steps	Time [s]	err	Ratio ( s a )	err	Ratio ( s a )
MTSM orig	40	9.8158 × 1 0 ‒ 3	—	—	—	—
MTSM opt	19	6.6631 × 1 0 ‒ 3	—	—	—	—
ode23	282,708	1.76957	7.45931 × 1 0 ‒ 7	180.28	8.39871 × 1 0 ‒ 7	265.58
ode45	11,773	2.31782 × 1 0 ‒ 2	7.45942 × 1 0 ‒ 7	2.36	8.39869 × 1 0 ‒ 7	3.48
ode113	468	6.847 × 1 0 ‒ 3	7.45943 × 1 0 ‒ 7	0.7	8.39869 × 1 0 ‒ 7	1.03

The plot of the order function P with step size scaling factor set to h scale = 2.5 and P min = 12 is in Figure 7. The automatic change in step size is visible at t = 2 s .

$Figure 7 Order function P P for benchmark problem B4, h = 0.5 s h=0.5\hspace{0.33em}{\rm{s}} , h scale = 2.5 {h}_{{\rm{scale}}}=2.5 , P min = 12 {P}_{\min }=12 .$

Figure 7

Order function P for benchmark problem B4, h = 0.5 s , h scale = 2.5 , P min = 12 .

The results show that the optimisations have a profound impact on the performance of the method and are beneficial. Presented optimisations can be freely combined to improve performance without any loss in accuracy.

7 Conclusion

The article dealt with the numerical solution of both linear and nonlinear IVPs represented by a system of ODEs. The MTSM outperforms the state-of-the-art MATLAB and PETSc ODE solvers in all cases. For linear problems, the wave equation, represented by large systems of ODEs, was solved in parallel using the method of lines. For nonlinear problems, optimisations were presented and demonstrated on the set of benchmark problems.

Our research is going to focus on effective solution of nonlinear systems of DAEs using the presented method.

Acknowledgement

The article includes the results of the internal BUT FIT project FIT-S-23-8151. This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254). The work extends the conference article [23] (DOI: 10.1109/Informatics57926.2022.10083462).

Funding information: The authors state that no funding is involved.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: Authors state no conflict of interest.
Data availability statement: Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

References

[1] J. Kunovský, Modern Taylor Series Method, Habilitation work, Faculty of Electrical Engineering and Computer Science, Brno University of Technology, 1994. Search in Google Scholar

[2] J. Kunovský, Tksl. online, https://www.fit.vut.cz/research/product/51/.en. Search in Google Scholar

[3] Matlab and simulink software, 2017, http://www.mathworks.com. Search in Google Scholar

[4] F. Kocina, Fos. online, http://www.fit.vutbr.cz/iveigend/fos. Search in Google Scholar

[5] V. Šátek, High performance computing research group, online, https://www.fit.vut.cz/research/group/hpc/.en. Search in Google Scholar

[6] F. Kocina, J. Kunovský, G. Nečasová, V. Šátek, and P. Veigend, “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, pp. 302–309, 2016. 10.1109/HPCSim.2016.7568350Search in Google Scholar

[7] R. Barrio, M. Rodríguez, A. Abad, and F. Blesa, “TIDES: A free software based on the Taylor series method,” Monografías de la Real Academia de Ciencias de Zaragoza, vol. 35, pp. 83–95, 2011. Search in Google Scholar

[8] A. Jorba and M. Zou, “A software package for the numerical integration of ODE by means of high-order Taylor methods,” Exp. Math., vol. 14, pp. 99–117, 2005. 10.1080/10586458.2005.10128904Search in Google Scholar

[9] Y. F. Chang and G. Corliss, “ATOMF: Solving ODEs and DAEs using Taylor series,” Computers & Mathematics with Applications, vol. 28, pp. 209–233, 1994. 10.1016/0898-1221(94)00193-6Search in Google Scholar

[10] M. Berz, COSY infinity version 8 Reference Manual, Technical Report MSUCL-1088, National Superconducting Cyclotron Lab., Michigan State University, East Lansing, Mich., 1997. Search in Google Scholar

[11] N. S. Nedialkov and J. Pryce, “Solving differential algebraic equations by Taylor series III. the DAETS code,” JNAIAM J. Numer. Anal. Ind. Appl. Math, vol. 3, pp. 61–80, 2008. Search in Google Scholar

[12] R. Barrio, F. Blesa, and M. Lara, “VSVO formulation of the Taylor method for the numerical solution of ODEs,” Computers and Mathematics with Applications, vol. 50, pp. 93–111, 2005. 10.1016/j.camwa.2005.02.010Search in Google Scholar

[13] R. Barrio, “Performance of the Taylor series method for ODEs/DAEs,” In: Applied Mathematics and Computation, vol. 163, pp. 525–545, 2005. ISSN 00963003. 10.1016/j.amc.2004.02.015Search in Google Scholar

[14] P. Mohazzabi and J. L. Becker, “Numerical solution of differential equations by Direct Taylor expansion,” Journal of Applied Mathematics and Physics, vol. 5, no. 3, 623–630, 2017. 10.4236/jamp.2017.53053Search in Google Scholar

[15] A. Baeza, S. Boscarino, P. Mulet, G. Russo, and D. Zorío, “Approximate Taylor methods for ODEs,” Computers & Fluids, vol. 159, pp. 156–166, 2017. 10.1016/j.compfluid.2017.10.001Search in Google Scholar

[16] P. Amodio, F. Iavernaro, F. Mazzia, M. S. Mukhametzhanov, and Y. D. Sergeyev, “A generalized Taylor method of order three for the solution of initial value problems in standard and infinity floating-point arithmeticm,” Mathematics and Computers in Simulation, vol. 141, pp. 24–39, 2017. 10.1016/j.matcom.2016.03.007Search in Google Scholar

[17] S. Dimova, I. Hristov, R. Hristova, I. Puzynin, T. Puzynina, Z. Sharipov, et al., OpenMP parallelization of multiple precision Taylor series method, 2019. https://arxiv.org/abs/1908.09301.Search in Google Scholar

[18] S. Liao, “On the reliability of computed chaotic solutions of non-linear differential equations,” Tellus A: Dynamic Meteorology and Oceanography, vol. 61, no. 4, pp. 550–564, 2008. 10.1111/j.1600-0870.2009.00402.xSearch in Google Scholar

[19] I. Hristov, R. Hristova, S. Dimova, P. Armyanov, N. Hegunov, I. Puzynin, et al., “Parallelizing multiple precision Taylor series method forï£¡integrating the Lorenz system,” In: I. Georgiev, H. Kostadinov, and E. Lilkova, editors, Advanced Computing in Industrial Mathematics, Springer International Publishing, Cham, pp. 56–66, 2023. 10.1007/978-3-031-20951-2_6Search in Google Scholar

[20] I. Hristov, R. Hristova, S. Dimova, P. Armyanov, N. Shegunov, I. Puzynin, et al., “On the efficient parallel computing of long term reliable trajectories for the Lorenz system,” arXiv:2101.06682v1, 2021. Search in Google Scholar

[21] P. Veigend, G. Nečcasová, and V. Šátek, “Taylor series based numerical integration method,” Open Computer Science, vol. 11, no. 1, pp. 60–69, 2021. 10.1515/comp-2020-0163Search in Google Scholar

[22] P. Veigend, G. Nečasová, and V. Šátek, “High order numerical integration method and its applications - the first 36 years of MTSM,” In: 2019 IEEE 15th International Scientific Conference on Informatics, Institute of Electrical and Electronics Engineers, pp. 25–30, 2019. 10.1109/Informatics47936.2019.9119258Search in Google Scholar

[23] P. Veigend, G. Nečasová, and V. Šátek, “Taylor series method in numerical integration: Linear and nonlinear problems,” In: 2022 IEEE 16th International Scientific Conference on Informatics, Informatics 2022 - Proceedings, IEEE, pp. 239–244. Communications Society, 2023. 10.1109/Informatics57926.2022.10083462Search in Google Scholar

[24] E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems. Springer-Verlag, Berlin Heidelberg, 1987. ISBN 3-540-56670-8. 10.1007/978-3-662-12607-3Search in Google Scholar

[25] G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, UK, 1985. Search in Google Scholar

[26] J. R. Dormand and P. J. Prince, “A family of embedded runge-kutta formulae,” Journal of Computational and Applied Mathematics, vol. 6, no. 1, pp. 19–26, 1980. 10.1016/0771-050X(80)90013-3Search in Google Scholar

[27] J. H. Verner, “Explicit Runge-Kutta methods with estimates of the local truncation error,” Society for Industrial and Applied Mathematics, vol. 15, pp. 772–790, 1978. 10.1137/0715051Search in Google Scholar

[28] IT4Innovations Documentation. Introduction. online, 2022, https://docs.it4i.cz/barbora/introduction/. Search in Google Scholar

[29] S. Abhyankar, J. Brown, E. M. Constantinescu, D. Ghosh, B. F. Smith, and H. Zhang. PETSc/TS: A Modern Scalable ODE/DAE Solver Library, 2018. Search in Google Scholar

[30] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, et al. PETSc Web page, 2024. https://www.mcs.anl.gov/petsc. Search in Google Scholar

[31] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, et al. PETSc Users Manual, Technical Report ANL-95/11 - Revision 3.12, Argonne National Laboratory, 2024, https://www.mcs.anl.gov/petsc. Search in Google Scholar

[32] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, “Efficient management of parallelism in object-oriented numerical software libraries,” In: E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, Birkhauser Press, Basel, pp. 163–202, 1997. 10.1007/978-1-4612-1986-6_8Search in Google Scholar

[33] S. Hamdi, W. E. Schiesser, and G. W. Griffiths, “Method of lines,” Scholarpedia, vol. 2, no. 7, pp. 1–11, 2007. 10.4249/scholarpedia.2859Search in Google Scholar

[34] R. Vichnevetsky, Numerical stability of methods of lines for partial differential equations, Rutgers University, New Jersey, 1971. 10.1177/003754977101600403Search in Google Scholar

[35] A. D. Malony, Metrics, Springer US, Boston, MA, pp. 1124–1130, 2011. Search in Google Scholar

[36] W. H. Enright and J. D. Pryce, “Two fortran packages for assessing initial value methods,” In: ACM Transactions on Mathematical Software, ACM, Mar, vol. 13, pp. 1–27. 1987. 10.1145/23002.27645Search in Google Scholar

[37] L. F. Shampine and M. W. Reichelt, “The matlab ode suite,” Journal of Numerical Analysis and Applied Mathematics, vol. 18, no. 1, pp. 1–22, 1997. 10.1137/S1064827594276424Search in Google Scholar

Received: 2024-01-15

Revised: 2024-03-09

Accepted: 2024-04-02

Published Online: 2024-07-12

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/comp-2024-0005

Keywords for this article

initial value problems; Taylor series; MTSM; MATLAB; PETSc

Creative Commons

BY 4.0