Numerical Quadratures Using the Interpolation Method of Hurwitz-Radon Matrices

Mathematics and computer sciences need suitable methods for numerical calculations of integrals. Classical methods, based on polynomial interpolation, have many weak sides: they are useless to interpolate the function that fails to be differentiable at one point or differs from the shape of polynomials considerably. We cannot forget about the Runge’s phenomenon. To deal with numerical interpolation and integration dedicated methods should be constructed. One of them, called by author the method of Hurwitz-Radon Matrices (MHR), can be used in reconstruction and interpolation of curves in the plane. This novel method is based on a family of Hurwitz-Radon (HR) matrices. The matrices are skew-symmetric and possess columns composed of orthogonal vectors. The operator of Hurwitz-Radon (OHR), built from that matrices, is described. It is shown how to create the orthogonal and discrete OHR and how to use it in a process of function interpolation and numerical integration. Created from the family of N-1 HR matrices and completed with the identical matrix, system of matrices is orthogonal only for vector spaces of dimensions N = 2, 4 or 8. Orthogonality of columns and rows is very significant for stability and high precision of calculations. MHR method is interpolating the curve point by point without using any formula of function. Main features of MHR method are: accuracy of curve reconstruction depending on number of nodes and method of choosing nodes; interpolation of L points of the curve is connected with the computational cost of rank O(L); MHR interpolation is not a linear interpolation.


Introduction
The following question is important in mathematics and computer sciences: is it possible to find a method of function interpolation and numerical integration without building the interpolation polynomials?Our paper aims at giving the positive answer to this question.This work is the effect of author studies on doctoral thesis.Current methods of numerical calculation of integrals are mainly based on classical polynomial interpolation: Newton, Lagrange or Hermite polynomials and spline curves which are piecewise polynomials [1] [2].Classical methods are useless to interpolate the function that fails to be differentiable at one point, for example the absolute value function ( ) at x = 0.If point (0;0) is one of the interpolation nodes, then precise polynomial interpolation of the absolute value function is impossible.Also when the graph of interpolated function differs from the shape of polynomials considerably, for example f(x) = 1/x, interpolation is very hard because of existing local extrema of polynomial.We cannot forget about the Runge's phenomenon: when interpolation nodes are equidistance then high-order polynomial oscillates toward the end of the interval, for example close to −1 and 1 with function f(x) = 1/(1 + 25x 2 ) [3].
This paper deals with the problem of interpolation [4] [5] and numerical quadrature without computing the polynomial or any fixed function.Values of nodes are used to build the orthogonal matrix operators and a linear (convex) combination of OHR operators leads to calculation of curve points.Main idea of MHR method is that the curve is interpolated point by point and computing the unknown coordinates of the points.The only significant factors in MHR method are choosing the interpolation nodes and fixing the dimension of HR matrix (N = 2, 4 or 8).Other characteristic features of function, such as shape or similarity to polynomials, derivative or Runge's phenomenon, are not important in the process of MHR interpolation.The curve is parameterized for value α ∈ [0;1] in the range of two successive interpolation nodes.In this paper computational algorithm is considered and then we have to talk about time.Complexity of calculations for one unknown point in Lagrange or Newton interpolation based on n nodes is connected with the computational cost of rank O(n 2 ).Complexity of calculations for L unknown points in MHR interpolation based on n nodes is connected with the computational cost of rank O(L).This is a very important feature of MHR method.

The Interpolation Method of Hurwitz-Radon Matrices
are called a family of Hurwitz-Radon matrices.A family of HR matrices in Equation ( 1) has important features: HR matrices are skew-symmetric ( ) 4 or 8 the family of Hurwitz-Radon matrices consists of N-1 matrices [6].
For N = 2 we have one matrix: For N = 4 there are three matrices with integer entries: For N = 8 we have seven matrices with elements 0, ±1 [7].Let's assume there is given a finite set of points of the function, called further nodes (x i , y i ) ∈ R 2 such as: 1) nodes (interpolation points) are settled at local extrema (maximum or minimum) and at least one point between two successive local extrema; 2) there are five nodes or more.Assume that the nodes belong to a curve in the plane.How the whole curve could be reconstructed using this discrete set of nodes?Proposed method [8] [9] is based on local, orthogonal matrix operators.Values of nodes' coordinates (x i , y i ) are connected with HR matrices [10] build on N dimensional vector space.It is important that HR matrices are skew-symmetric and only for dimension N = 2, 4 or 8 columns and rows of HR matrices are orthogonal [11].
If one curve is described by a set of nodes ( )

{ }
, , 1, 2, , i i x y i n =  then HR matrices combined with identity matrix are used to build an orthogonal and discrete Hurwitz-Radon Operator (OHR).For nodes (x 1 , y 1 ), (x 2 , y 2 ) OHR of dimension N = 2 is constructed: For nodes (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ), (x 4 , y 4 ) OHR of dimension N = 4 is constructed: where For nodes ( ) ( ) ( ) , , , , , ,  x y x y x y  OHR M of dimension N = 8 is built [8] similarly as Equation ( 3): The components of the vector ( ) , appearing in the matrix M in Equation ( 4), are defined by Equation ( 5) in the similar way to Equations ( 2) and (3) but in terms of the coordinates of the above 8 nodes.Note that OHR operators in Equations ( 2)-( 4) satisfy the condition of interpolation , , , , , , How can we compute coordinates of points settled between interpolation nodes?On a segment of a line every number "c" situated between "a" and "b" is described by a linear (convex) combination ( ) Average OHR operator M 2 of dimension N = 2, 4 or 8 is constructed as follows: ( ) with the operator M 0 built in Equations ( 2)-( 4) by "odd" nodes ( ) ( ) ( ) , , , , , , built in Equations ( 2)-( 4) by "even" nodes ( ) ( ) ( ) x b y x y x y =  .Notice that having the operator M 2 for coordinates x i < x i+1 it is possible to reconstruct the second coordinates of points (x, y) in terms of the vector C defined with ( ) as . The required formula is adequate to Equation ( 6): in which components of vector Y(C) give us the second coordinate of the points (x, y) corresponding to the first coordinate, given in terms of components in Equation ( 9) of the vector C.
After computing in Equations ( 7)- (10) for any α ∈ [0;1], we have a half of reconstructed (j = 1 in Algorithm 1).Now it is necessary to find second half of unknown coordinates (j = 2 in Algorithm 1) for ( ) There is no need to build the OHR for nodes ( ) ( ) ( ) x a y x y x y =  because we just find M 1 .This operator will play as role as M 0 in Equation ( 8).New M 1 must be computed for nodes ( ) ( ) x b y x y , . As we see the minimum number of interpolation nodes n = 2N + 1 = 5, 9 or 17 using OHR operators of dimension N = 2, 4 or 8 respectively.If there is more nodes than 2N + 1, the same calculations in Equations ( 7)- (11) have to be done for next range(s) or last range of 2N + 1 nodes.For example, if n = 9 then we can use OHR operators of dimension N = 4 or OHR operators of dimension N = 2 for two subsets of nodes: ( ) ( ) { } , , , , x y x y  and ( ) ( ) , , , , x y x y  .We summarize this section in the following algorithm of points reconstruction for 2N+1 = 5, 9 or 17 successive nodes.
Step 3. Determine the number of points to be reconstructed K j > 0 between two successive nodes, let k = 1.
The number of reconstructed points in Algorithm 1 is . If there is more nodes than 2N + 1 = 5, 9 or 17, Algorithm 1 has to be done for next range(s) or last range of 2N + 1 nodes.Reconstruction of curve points using Algorithm 1 is called by author the method of Hurwitz-Radon Matrices (MHR).

MHR Calculations and Graphical Examples
In this section we consider the number of multiplications and divisions for MHR method during reconstruction of K = L -n points having n interpolation nodes of the curve consists of L points.First we present a formula for computing one unknown coordinate of a single point.Assume there are given four nodes (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ) and (x 4 , y 4 ).OHR operators of dimension N = 2 are built in Equation (2) as follows: x y x y x y x y M x y x y x y x y x x Let first coordinate c 1 of reconstructed point is situated between x 1 and x 2 : Compute second coordinate of reconstructed point y(c 1 ) for Y(C) = [y(c 1 ), y(c 2 )] T from Equation ( 10): ( ) ( ) ( ) After calculation by Equation ( 13): ( ) ( ) ( ) . y c y y x x y x x y x x y x x y x x x x y x x y x x y x x y x x So each point of the curve P = (c 1 , y(c 1 )) settled between nodes (x 1 , y 1 ) and (x 2 , y 2 ) is parameterized by P(α) for Equations ( 12) and ( 14) and α ∈ [0;1].
If nodes (x i , y i ) are equidistance in coordinate x i , then parameterization of unknown coordinate ( 14) is simpler.Let four successive nodes (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ) and (x 4 , y 4 ) are equidistance in coordinate x i and a = x 1 , h/2 = x i+1 -x i = const.Calculations in Equations ( 13) and ( 14) are done for c 1 in Equation ( 12): ( ) As we can see in Equations ( 15) and ( 16), MHR interpolation is not a linear interpolation.It is possible to estimate the in terpolation error of MHR method (Algorithm 1) for the class of linear function f: Notice that estimation from Equation ( 17) has the biggest value 0.25s for β = α = 0.5, when c 1 is situated in the middle between x 1 and x 2 .The goal of this paper is not a reconstruction of single point, like for example Equations ( 14) and ( 15), but interpolation of curve consists of L points.If we have n interpolation nodes, then there is K = L -n points to find using Algorithm 1 and MHR method.Now we consider the complexity of MHR calculations.
Lemma 1.Let n = 5, 9 or 17 is the number of interpolation nodes, let MHR method (Algorithm 1) is done for reconstruction of the curve consists of L points.Then MHR method is connected with the computational cost of rank O(L).
Proof. □ The lowest computational cost appears in MHR method with five nodes and OHR operators of dimension N = 2. Therefore whole set of n nodes can be divided into subsets of five nodes.Then whole curve is to be reconstructed by Algorithm 1 with all subsets of five nodes: ( ) ( ) { }  in Algorithm 1. Function f(x) = 1/x is an example when the graph of interpolated function differs from the shape of polynomials considerably.Then classical interpolation is very hard because of existing local extrema of polynomial (Figure 1).Here is the application of Algorithm 1 for this function and five nodes.
Figure 2 contains not many interpolated points (twenty six) and minimal number of nodes (five), so numerical calculations of integral (precise value I = 2.196) by trapezoidal rule I 1 = 2.213 are not always satisfying.Greater number of nodes and interpolated points gives us more accurate value of quadrature.Other examples of MHR interpolation and numerical quadratures: Figure 3 contains minimal number of nodes (five) and not many interpolated points (twenty two), so numerical calculations of integral (precise value I = 0.549) are not always satisfying: 1) trapezoid method: I 1 = 0.534; 2) Simpson's rule: I 2 = 0.538.Greater number of nodes and interpolated points gives us more accurate value of quadrature.Figure 4 contains minimal number of nodes (five) and not many interpolated points (thirty six), but numerical calculations of integral (precise value I = 1.029) are interesting: 1) trapezoid method: I 1 = 1.000; 2) Simpson's rule: I 2 = 0.999.After computing of K points for interpolated function (algorithm 1), it is possible to calculate the integral by trapezoid method or Simpson's rule in range of each pair of successive interpolated points.Greater number of nodes and interpolated points gives us more accurate value of quadrature.
Figures 5-8 are the examples of interpolation for such functions as polynomial, logarithmic and exponential.

Summary
The method of Hurwitz-Radon Matrices (MHR-Algorithm 1) leads to curve interpolation depending on the number of nodes and location of nodes.No characteristic features of curve are important in MHR method: failing to be differentiable at any point, the Runge's phenomenon or differences from the shape of polynomials.These features are very significant for classical polynomial interpolations.MHR method gives the possibility of curve reconstruction and then numerical calculations of integrals for interpolated function via trapezoid method or Simpson's rule for example.The only condition is to have a set of nodes according to assumptions in Algorithm 1. Curve modeling [12] by MHR method is connected with possibility of changing the nodes coordinates   and reconstruction of new curve or contour for new set of nodes, no matter what shape of curve or contour is to be reconstructed.Main features of MHR method are: 1) accuracy of curve reconstruction depending on number of nodes and method of choosing nodes; 2) Reconstruction of curve consisting of L points is connected with the computational cost of rank O(L); 3) Algorithm 1 is dealing with local operators: average OHR operator M 2 (8) is built by successive 4, 8 or 16 nodes, what is connected with smaller computational costs then uses all nodes.Future works are connected with: geometrical transformations of curve (translations, rotations, scaling)only nodes are transformed and new curve (for example contour of the object) for new nodes is reconstructed; estimation of curve length [13]; possibility to apply MHR method to three-dimensional curves [8] [9]; object recognition [14], shape representation [15] [16] and parameterization [17] in image processing; curve extrapolation [18].

Figure 1 .
Figure 1.Lagrange interpolation polynomial of function f = 1/x differs extremely from the original.

Figure 2 .
Figure 2.Twenty six interpolated points of function f = 1/x.Lagrange interpolation polynomial for function f(x) = 1/x and nodes (5;0.2),(5/3;0.6),(1;1), (5/7;1.4),(5/9;1.8)has one local minimum and two roots.Other examples of MHR interpolation and numerical quadratures: Figure3contains minimal number of nodes (five) and not many interpolated points (twenty two), so numerical calculations of integral (precise value I = 0.549) are not always satisfying:1) trapezoid method: I 1 = 0.534; 2) Simpson's rule: I 2 = 0.538.Greater number of nodes and interpolated points gives us more accurate value of quadrature.Figure4contains minimal number of nodes (five) and not many interpolated points (thirty six), but numerical calculations of integral (precise value I = 1.029) are interesting:1) trapezoid method: I 1 = 1.000; 2) Simpson's rule: I 2 = 0.999.After computing of K points for interpolated function (algorithm 1), it is possible to calculate the integral by trapezoid method or Simpson's rule in range of each pair of successive interpolated points.Greater number of nodes and interpolated points gives us more accurate value of quadrature.Figures5-8are the examples of interpolation for such functions as polynomial, logarithmic and exponential.