Piecewise interpolation solution of ordinary differential equations with application to numerical modeling problems

Piecewise interpolation approximation of functions of one real variable, derivatives and integrals is constructed with the help of Lagrange and Newton interpolation polynomials. The polynomials are transformed to the form of algebraic polynomials with numerical coefficients by means of restoring the polynomial coefficients by its roots. Formulas different from Vieta’s formulas are applied. In the resulting form, the polynomials interpolate the right-hand sides of ordinary differential equations, the expression of the antiderivative ones is used to approximate the solution. Iterative refinement is performed. Error estimates and results of numerical experiments for stiff and non-stiff problems in physical, chemical models and other processes are presented.


Introduction and formulation of the question
An approximate solution of the Cauchy problem for ordinary differential equations (ODEs) is constructed on the basis of Lagrange and Newton interpolation polynomials. Both polynomials are equivalently transformed to the form of algebraic polynomials with numerical coefficients. In this case, the analytical expression of the derivative and the antiderivative is obtained. With the help of transformed polynomials, the right-hand sides of each equation of the ODE system are approximated. The antiderivative of the polynomial with the corresponding value of the constant -the coordinate of the approximate solution, is substituted into the right -hand side of the function. The process is iteratively repeated, and as a result, the solution is approximated with relatively high accuracy. The task of the message is to present a mathematical description of the method, its algorithmization and software implementation. The method is combined with piecewise interpolation. Numerical stability, error bounds, and time complexity for stiff and non-stiff problems are analyzed. A comparison is made with the methods of Runge-Kutta, Butcher and Dormand-Prince [1,2]. The aim of the work is to show the distinctive quality of the proposed method, which consists in a limited accumulation of error while maintaining an acceptable labor intensity, which is in demand in problems of physical and chemical processes' modeling [3] and in problems of planetary astronomy [4]. A numerical experiment is described, and the possibility of reducing the error without automatic selection of parameters described in [5] is investigated. Together with the simplification of the method, it is necessary to reduce its labor intensity.
. When n k = the left-hand side parts of (2) will coincide with the desired values of the polynomial coefficients (1). The algorithm is stored for complex roots and coefficients, and is programmed in the form of a double cycle [6]. Relations (2) are used to transform interpolation polynomials. For the interpolation of the function , the Lagrange polynomial has the form: where r x are interpolation nods, based on numerical experiment, the computer implementation gives a lower error in the case of equidistant nodes.

Transformation of the Lagrange interpolation polynomial
segment equidistant nodes for the interpolation polynomial (3) are taken: A similar transformation is given in [7]. The difference will be in the conversion of the numerators from (4) to the form of polynomials with integer coefficients. As a result, analytical expression of the antiderivatives, derivatives and organization of the iterations in the right-hand sides of the ODE will The coefficients of the polynomial (6) will be integers. They can be calculated a priori and stored in the computer memory as constants that do not depend on the interpolated function ) (x f , on the range and location of its argument, and this can be done for any j and all degrees of the n polynomial in a priori given boundaries. The denominator in (4) differs from the numerator in that it has j t = in it. As a result the numerator and denominator in (7) are conveniently calculated according to the Horner's method (scheme) In this notation If we collect the coefficients with equal degrees of terms then where  D is the result of this type. Further, the approximation ) ( is performed with the help of (8), (9). Hence two varieties of the derivative approximation: In practice, (10) is somewhat more accurate than (11).

Piecewise interpolation
If (8) -(11) are applied on small subintervals of equal length with common partition boundaries then the accuracy of the approximation will increase: the function and integral will be approximated with accuracy to   [7]. The relation [8] is valid for 0 ≥ ∀ k : where it is assumed to be that ) (x f is defined and continuously differentiable at is the Newton interpolation polynomial for forward interpolation on the subinterval . Both these statements and (13) as well are preserved for the considered transformations of the Lagrange interpolation polynomial.

Adding (16) over all subintervals implies
If on each segment in the right-hand side (14) ) (x f approximates with the estimate (13), also valid for the considered transformation of the Lagrange polynomial, where The estimation (18) takes place under the conditions described above, under which (13) is fulfilled.
. Formula (17) can be interpreted as an analytical version of the Newton-Cotes formulas [7].

Piecewise interpolation solution of the Cauchy problem for ODEs with iterative refinement
Let the Cauchy problem be considered in the domain is defined, it is continuously differentiable (at points a -in the right-hand side, b -in the left-hand side) and satisfies the Lipschitz condition: . It is assumed that in R the solution of problem (19) exists and is unique. For simplicity of notation, a , b are the same as in (12). To interpolate the right-hand side of (19), an approximate value y is substituted in is approximated by polynomials of the form in (8) , is taken as an approximation of the solution: , and for the same n , on the same segment, an interpolation polynomial of the form (9) is again constructed: . Again, the antiderivative with the same value of the constant , which is then similarly interpolated: continue up to a given boundary const = ≤ q  , abstractly their number is not limited. The formula was taken for the above value is taken for i y 0 . Error estimates will be performed below under additional assumptions. Given that for equidistant nodes on the same segment, when interpolating the same function, the Lagrange and Newton interpolation polynomials coincide to the same degree in all the considered forms, to simplify the notation, it is assumed that the interpolation is performed by a polynomial of the form (3). More precisely, the interpolation is assumed without transition to a variable t , while the polynomial has fully reduced numerical coefficients: such a Lagrange polynomial is written as The interpolation error on ] , [ . It will be assumed that the solution on this segment corresponds not to the initial conditions given in (19), but to the initial conditions on this particular segment: In this case, the solution ) (x y and the have the same initial values on the given segment. Therefore, (20) is taken into account, so that no replacement of the variable is required. First, it is assumed that on a sufficiently large segment of the sequence of numbers  taking into account the same initial values of the exact solution and its approximation, the absolute error of the  iteration will take the By construction )) ( , ( ) ( , the interpolation error is denoted by k i с , in this notation Applying the Lipschitz condition, arbitrariness in both parts of the inequality, you can shift to the maximum: Let the left-hand side of inequality (25) be denoted by  i ε , then (26) From (26), (28) The inequality (29) is a consequence of (25) and is true for those consecutive ones   , 2 , 1 , 0 = for which (21) is not violated. Without detracting from generality, we can assume that 1 ) . Therefore, inequality (21) will prove to be violated. Violation of (21) means that (31) Hence, taking into account the implementation of (29) for 1 0 + =   and the fact that the integrand in (31) can be replaced by the corresponding expression in (24): But for the right-hand side (32), after shifting to the maximum in both parts of the inequality, the estimate (29) is preserved, at which the left-hand side of the inequality (29) already violates (31) and does not exceed k i c . Therefore, the left-hand side of (32) will not exceed k i c : By analogy, , therefore, the integrand on the righthand side of (34) can be replaced by the corresponding integrand from (24), which is the same as the integrand from (32): ]. , Repeating the previous reasoning entails By obvious induction, inequalities (33) and (35) will become inequality where r is an arbitrary natural number.
From the stated above we derive Lemma 1.  (12), there is a number 0 r such that the iterative refinement will satisfy the inequality 0 ) ( (36) For numbers from (36), (23) is preserved, whence by applying the Lipschitz condition (37) The estimate (38) means that the application of iterative refinement allows us to approximate the solution on the entire segment (12) with an absolute error of piecewise interpolation on one, separately taken, subinterval from (12), up to a constant multiplier ) ( a b N − . Substituting in (38) k i с from (13) implies The above arguments and estimates are transferred to the case of the Newton interpolation polynomial, represented in the form of an algebraic polynomial with numerical coefficients.
The theorem and corollaries give a formal estimate of the error under abstract conditions involving the existence of 1 + n derivative of the (19) function on the right-hand side. However, interpolation itself is possible under the broadest conditions, so the method can always be applied in conditions of unique existence. In practice, much is determined not only by the smallness of subinterval in (39), but by a view of the right-hand side of (19), the stability of the solution in the sense of Lyapunov, stiffness or non-stiffness class of problems. Nevertheless, in all such experimentally considered cases, the proposed method has a lower error than the well-known methods, differing in a limited accumulation of error with an increase in the interval of the approximate solution.

Numerical Experiment
Part of the experiment is given on the examples of problems with analytical solutions. The results of chemical processes modeling and problems of celestial mechanics are also discussed. The experiment was conducted using a PC based on an Intel(R) Core(TM) i5-2500.    Table 2. Error and number of calls to the right-hand side when solving the problem of Example 2 3 × ≈ fc [1]. The Gauss-Everhart integrator of the 19th order -GAUSS_32 [9,10], implemented in the Delphi environment, reaches the lowest error limit, of the order 16 10 − , on the same segment, among the studied methods of numerical approximation. The method is adapted for solving problems of celestial mechanics, in particular, the mechanism for choosing the size of the integration step is implemented taking into account the specific nature of the planar two-body problem [10]. The piecewise interpolation solution of the problem (40)  with variations in the degree of the interpolation polynomial and the number of subintervals of the partition are presented in [5]. Fixing the method parameters and excluding additional clarifying procedures of the program allowed to reduce the time for solving the problem with the error bound 14 10 − from min 11 ≈ t to ms 677 14s

Conclusion
An approximate solution of the Cauchy problem for ODE is constructed by applying piecewise interpolation based on Lagrange and Newton interpolation polynomials, equivalently transformed to the form of algebraic polynomials with numerical coefficients. With the help of transformed polynomials, the right-hand sides of the equations of the system are approximated. The antiderivatives of polynomials are the coordinates of the approximate solutions. The process is iteratively refined, and