Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials

Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials ∗

A. Messaoudi a , , R. Sadaka a , H. Sadok b a

Ecole Normale Supérieure de Rabat, Laboratory of Mathematics, Computing and Applications - Information Security (LabMiA-SI), B.P. 5118, Mohammed V University in Rabat, Morocco b Laboratoire de Mathématiques Pures et Appliquées, Université du Littoral Côte d’Opale, 50 Rue F. Buisson, BP 699 - 62228 Calais cedex, France

article

info

a b s t r a c t

Article history: Received 31 December 2018 Received in revised form 13 May 2019 MSC: 65F10 65F30 Keywords: Polynomial interpolation Hermite interpolation polynomials Matrix recursive interpolation algorithm Matrix recursive polynomial interpolation algorithm Generalized recursive polynomial interpolation algorithm

Let x0 , x1 , . . . , xn , be a set of n+1 distinct real numbers (i.e., xi ̸ = xj , for i ̸ = j) and ym,k , for m = 0, 1, . . . , n, and k = 0, 1, . . . , nm , with nm ∈ N, be given of real numbers, ∑n we know that there exists a unique polynomial pN −1 of degree N − 1 where N = i=0 (ni + 1), (k) such that pN −1 (xm ) = ym,k , for m = 0, 1, . . . , n and k = 0, 1, . . . , nm . pN −1 is the Hermite interpolation polynomial for the set {(xm , ym,k ), m = 0, 1, . . . , n, k = 0, 1, . . . , nm }. The polynomial pN −1 can be computed by using the Lagrange generalized polynomials. Recently Messaoudi et al. (2018) presented a new algorithm for computing the Hermite interpolation polynomials, for a general case, called Generalized Recursive Polynomial Interpolation Algorithm (GRPIA), this algorithm has been developed without using the Matrix Recursive Interpolation Algorithm (Jbilou and Messaoudi, 1999). Messaoudi et al. (2017) presented also a new algorithm called Matrix Recursive Polynomial Interpolation Algorithm (MRPIA), for a particular case where nm = µ = 1, for m = 0, 1, . . . , n. In this paper we will give the version of the MRPIA for a particular case nm = µ ≥ 0, for m = 0, 1, . . . , n. We will recall the result of the existence of the polynomial pN −1 for this case, some of its properties will also be given. Using the MRPIA, a method will be proposed for the general case, where nm , for some m, are different and some examples will also be given. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The general Hermite interpolation polynomial problem is defined as follows [1,2]: given that n is a non negative integer, let xm , m = 0, 1, . . . , n, be n + 1 distinct real numbers (i.e., xi ̸ = xj , for i ̸ = j) and ym,k , for m = 0, 1, . . . , n and k = 0, 1, . . . , nm , be given real numbers. The Hermite interpolation problem ∑n for this data consists of determining a polynomial pN −1 ∈ PN −1 whose degree does not exceed N − 1, where N = i=0 (ni + 1), and PN −1 denotes the set of all real-valued polynomials of degree ≤ N − 1, defined over the set R of real numbers, and which satisfies the following interpolation conditions (k)

pN −1 (xm ) = ym,k ,

m = 0, 1, . . . , n, k = 0, 1, . . . , nm .

(1)

∗ Corresponding author. E-mail addresses: [email protected] (A. Messaoudi), [email protected] (R. Sadaka), [email protected] (H. Sadok). https://doi.org/10.1016/j.cam.2019.112471 0377-0427/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

2

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

This problem differs from the usual interpolation problem for polynomials in that it prescribes at each support abscissa xm not only the value but also the first nm derivatives of the desired polynomial. There are exactly N conditions (1) for the N coefficients of the interpolating polynomial, leading us to expect that the Hermite interpolation problem can be solved uniquely. The case where nm = 0, for m = 0, 1, . . . , n, gives the Lagrange polynomial interpolation problem [2], for solving this problem a new algorithm has been proposed in [3], called the recursive polynomial interpolation algorithm (RPIA), which is deriving from the recursive interpolation algorithm (RIA), given by Brezinski in [4], the interpolation polynomial has been expressed as a Schur complement [5–7] and using the Sylvester identity [3] the algorithm is given. The case where nm = 1, for m = 0, 1, . . . , n, has been studied by Messaoudi et al. [8], a new algorithm has been developed for solving the Hermite polynomial interpolation problem. This algorithm is called the matrix recursive polynomial interpolation algorithm (MRPIA), which is deriving from the matrix recursive interpolation algorithm (MRIA), given in [9]. We will recall the new formulation of the general Hermite polynomial interpolation problem. A result of the existence of the Hermite polynomial interpolation, for a particular case nm = µ ≥ 0, for m = 0, 1, . . . , n, will be given. We will recall the MRPIA and some of its properties for solving this problem for this particular case. For the case where the derivation orders are different (ni ̸ = nj , for some i ̸ = j), a method will be proposed. Some examples will also be given. This paper will be organized as follows. In Section 2 we recall the formulation of the Hermite polynomial interpolation problem for the general case, we give a result for the existence of the polynomial pN −1 for a particular case nm = µ ≥ 0. In Section 3, we give the version of the MRPIA for this particular case, we also recall some of its properties. In Section 4, we show how to solve the Hermite polynomial interpolation problem for the general case where nm are different for some m. Section 5 is concerned with some examples. 2. General Hermite polynomial interpolation We recall the new formulation of the general Hermite polynomial interpolation problem [10]. We assume that x0 , x1 , . . . , xn are given distinct real numbers and ym,k , for m = 0, 1, . . . , n and k = 0, 1, . . . , nm , are given of real numbers. 2.1. The formulation of the Hermite interpolation The general Hermite polynomial interpolation problem can be defined as follows: Let qm , for m = 0, 1, . . . , N − 1, be polynomials of PN −1 of degree m, we assume that q0 (x) ̸ = 0. Let cm , for m = 0, 1, . . . , n, be the functional defined by cm (g (k) (x)) = g (k) (xm ) for any given function g, we assume that the kth derivative of this function exists (with the convention g (0) (x) = g(x), for all x). Find the polynomial pN −1 ∈ PN −1 , of degree N − 1, such that pN −1 (x) =

N −1 ∑

αi qi (x),

(2)

i=0

and for m = 0, 1, . . . , n and k = 0, 1, . . . , nm (k)

(k)

cm (pN −1 (x)) = pN −1 (xm ) = ym,k .

(3)

Remark 2.1. If nm = 0, for m = 0, 1, . . . , n, then the Hermite polynomial interpolation will be the well known the Lagrange polynomial interpolation. The real numbers ym,k , for m = 0, 1, . . . , n and k = 0, 1, . . . , nm , can be given as the values of a real-valued function f and its derivative f (k) , defined on a closed real interval [a, b], at the distinct interpolation points xm ∈ [a, b], m = 0, 1, . . . , n, then pN −1 , defined by (2) and (3) is the Hermite interpolation polynomial of degree N − 1 for f . 2.2. The case where nm = µ ≥ 0 Now we will consider the particular case nm = µ ≥ 0, for m = 0, 1, . . . , n. For this case we have N −1=

n ∑

(µ + 1) − 1 = (µ + 1)n + µ.

i=0

So for this case we have to find the polynomial pN −1 ∈ PN −1 , of degree N − 1, such that pN −1 (x) =

N −1 ∑

αi qi (x),

(4)

i=0

and for m = 0, 1, . . . , n, and k = 0, . . . , µ (k)

(k)

cm (pN −1 (x)) = pN −1 (xm ) = ym,k .

(5)

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

3

Now we will show that the solution of the problem (4)–(5) exists and is unique. From the relations (4) and (5) we get the following system q0 (x0 )

⎛

qµ (x0 )

··· .. .

⎜ .. ⎜ . ⎜ (µ) ⎜q0 (x0 ) ⎜ ⎜ .. ⎜ . ⎜ ⎜ q0 (xn ) ⎜ ⎜ .. ⎝ .

.. .

(µ) qµ (x0 )

··· .. . ··· .. .

(µ)

q0 (xn )

.. .

qµ (xn )

.. .

q(µµ) (xn )

···

··· .. .

qN −µ−1 (x0 )

··· .. . ··· .. .

qN −µ−1 (x0 )

.. .

··· .. . ··· .. .

···

qN −µ−1 (xn )

···

(µ)

.. . .. .

qN −µ−1 (xn ) (µ)

qN −1 (x0 )

··· .. .

⎞ ⎛ ⎞ α0 y0,0 ⎟⎜ . ⎟ ⎜ . ⎟ ⎟ ⎜ .. ⎟ ⎜ .. ⎟ ⎟⎜ (µ) ⎟ ⎜ ⎟ qN −1 (x0 )⎟ ⎜ αµ ⎟ ⎜y0,µ ⎟ ⎟⎜ ⎟ ⎜ ⎟ ⎜ .. ⎟ ⎜ .. ⎟ .. ⎟⎜ . ⎟ = ⎜ . ⎟ ⎟. . ⎟⎜ ⎟ ⎜ ⎟ qN −1 (xn )⎟ α y ⎜ n ,0 ⎟ N −µ−1 ⎟ ⎟⎜ ⎜ ⎟ ⎜ ⎟ ⎝ .. ⎠ ⎝ .. ⎟ .. ⎠ ⎠ . . . (µ) αN −1 yn,µ qN −1 (xn ) .. .

⎞⎛

(6)

If the matrix of (6), which we denote by VN −1 , is nonsingular, then pN −1 exists and is unique. Now we will give a result about the existence of this polynomial. 2.3. The existence of pN −1 We set for m = 0, . . . , n, and for k = 0, 1, . . . , µ, q0 (x0 )

⎛

V(µ+1)m+k

⎜ .. ⎜ . ⎜ (µ) ⎜q0 (x0 ) ⎜ ⎜ = ⎜ ... ⎜ ⎜ q0 (xm ) ⎜ ⎜ .. ⎝ .

(k) q0 (xm )

qµ (x0 )

··· .. .

··· .. .

q(µ+1)m (x0 ) q(µ+1)m (x0 )

.. .

··· .. . ··· .. .

q(k) µ (xm )

···

.. .

(µ) qµ (x0 )

··· .. . ··· .. .

.. .

qn0 (xm )

···

(µ)

··· .. .

.. .

··· .. . ··· .. .

.. .

q(µ+1)m (xm )

.. .

(k) q(µ+1)m (xm )

···

q(µ+1)m+k (x0 )

⎞

.. .

⎟ ⎟ ⎟ (µ) q(µ+1)m+k (x0 ) ⎟ ⎟ ⎟ .. ⎟. . ⎟ q(µ+1)m+k (xm )⎟ ⎟ ⎟ .. ⎠ .

(7)

(k)

q(µ+1)m+k (xm )

We have the following result, which is a particular case of the result given in [10]. Lemma 2.1.

For m = 1, . . . , n, and k = 0, . . . , µ, we have m−1

|V(µ+1)m+k | = β(µ+1)m+k

∏

(xm − xj )(k+1)(µ+1) ×

∏

(xi − xj )(µ+1) , 2

0≤j
j=0

with

β(µ+1)m+k = k!a(µ+1)m+k β(µ+1)m+k−1 , where ai is the coefficient of xi in qi

βµ = |Vµ | =

µ ∏

k!ak ,

k=0

and where |Z | denotes the determinant of the square matrix Z . Remark 2.2. As can be seen VN −1 is a strongly nonsingular matrix because xi ̸ = xj for i ̸ = j and |Vi | ̸ = 0, for i = 0, 1, . . . , N − 1, see [10]. 3. The MRPIA for the case where nm = µ ≥ 0 So from (4) and (5) and for k = 0, 1, . . . , µ, the kth derivative of the Hermite interpolation polynomial pN −1 ∈ PN −1 , can be written as follows

⎡

(k)

[

(k)

pN −1 (x) = q0 (x)

···

(k) qµ (x)

···

y0,0

⎤

⎢ .. ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢y0,µ ⎥ ] ⎢ . ⎥ (k) −1 ⎥ qN −1 (x) VN −1 ⎢ ⎢ .. ⎥ , ⎢ ⎥ ⎢ yn,0 ⎥ ⎢ . ⎥ ⎣ . ⎦ .

(8)

yn,µ Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

4

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx (k)

We see that pN −1 (x) can be expressed as a Schur complement [8]

⎛⎡

(k)

q0 (x) q0 (x0 )

0

⎜⎢ y ⎜ ⎢ 0 ,0 ⎜⎢ . ⎜⎢ .. ⎜⎢ ⎜⎢ ⎜⎢ y0,µ (k) ⎢ pN −1 (x) = − ⎜ ⎜⎢ .. ⎜⎢ . ⎜⎢ ⎜⎢ yn,0 ⎜⎢ . ⎜⎢ . ⎝⎣ .

.. .

(µ)

q0 (x0 )

.. .

··· .. . ··· .. .

q0 (xn )

···

.. .

q0 (xn ) (µ)

yn,µ

q(k) µ (x) qµ (x0 )

··· ··· .. .

··· ···

.. .

qµ (xn )

.. .

q(µµ) (xn )

⎞

··· .. . ··· .. .

⎟ ⎟ ⎥ ⎟ ⎥ ⎟ ⎥ ⎟ ⎥ ⎟ (µ) qN −1 (x0 ) ⎥ ⎟ ⎥ ⎟ .. ⎥ /VN −1 ⎟ . ⎥ ⎟ . ⎥ ⎟ ⎟ qN −1 (xn ) ⎥ ⎥ ⎟ .. ⎥ ⎟ ⎦ ⎠ .

···

qN −1 (xn )

.. .

(µ) qµ (x0 )

⎤

(k)

qN −1 (x) qN −1 (x0 ) ⎥ ⎥

(9)

(µ)

3.1. The algorithm Now we will recall how the MRPIA has been developed in [8]. We set

⎧ ( )T 0 · · · 0 ∈ Rµ+1 , and for m = 0, 1, . . . , n 0= 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ( )T ⎪ ⎪ ⎪ ym,1 · · · ym,µ ∈ Rµ+1 , Ym = ym,0 ⎪ ⎪ ⎪ ⎪ ⎪ ( )T ⎪ ⎪ (µ) ⎪ ′ ⎪ Pm (x) = p(µ+1)m+µ (x) , p (x) · · · p (x) ⎪ (µ+1)m+µ (µ+1)m+µ ⎪ ⎪ ⎪ ⎨ ⎛ ⎞ q(µ+1)m (x) · · · q(µ+1)m+µ (x) ⎪ ⎪ ⎜ ⎟ ⎪ ⎪ ⎜ ′ ⎟ ′ ⎪ ⎪ q (x) · · · q (x) ⎜ ⎟ ⎪Q (x) = ⎜ (µ+1)m (µ+1)m+µ ⎪ ⎟, m ⎪ . . . ⎪ ⎜ ⎟ ⎪ .. .. .. ⎪ ⎝ ⎠ ⎪ ⎪ ⎪ (µ) (µ) ⎪ ⎪ q (x) · · · q (x) (µ+1)m (µ+1)m+µ ⎪ ⎪ ⎪ ⎪ ⎪ ( ) ( ) ⎩ Dm = V(µ+1)m+µ = ci Qj (x) 0≤i,j≤m = Qj (xi ) 0≤i,j≤m . We remark that (8) can be written as follows 0 ⎜⎢ Y0 ⎜⎢ Y ⎢ 1 Pn (x) = − ⎜ ⎜⎢

Q0 (x) Q0 (x0 ) Q0 (x1 )

Yn

Q0 (xn )

⎛⎡

⎝⎣ ...

Q1 (x) Q1 (x0 ) Q1 (x1 )

.. .

.. .

··· ··· ··· .. .

Qn (x) Qn (x0 ) ⎥ ⎟ ⎟ ⎥ Qn (x1 ) ⎥ /D ⎟ . ⎥ n⎟

Q1 (xn )

···

Qn (xn )

.. .

⎤

⎞

⎦

⎠

(10)

For computing recursively Pn , we set for m = 0, 1, . . . , n, 0 ⎜⎢ Y 0 ⎜⎢ Y ⎢ 1 Pm (x) = − ⎜ ⎜⎢

Q0 (x) Q0 (x0 ) Q0 (x1 )

Ym

Q0 (xm )

⎛⎡

⎝⎣ ...

.. .

Q1 (x) Q1 (x0 ) Q1 (x1 )

.. .

··· ··· ··· .. .

Qm (x) Qm (x0 ) ⎥ ⎟ ⎥ ⎟ Qm (x1 ) ⎥ /D ⎟ , ⎥ m⎟

Q1 (xm )

···

Qm (xm )

.. .

⎤

⎞

⎦

⎠

(11)

with P−1 (x) = 0, and for i = 0, 1, . . . , n − 1, and m > i, we need the following auxiliary matrix polynomials Qm (x) ⎜⎢ Qm (x0 ) ⎜⎢ Q (x ) ⎢ m 1 Gi,m (x) = ⎜ ⎜⎢

⎛⎡

⎝⎣

.. .

Qm (xi )

Q0 (x) Q0 (x0 ) Q0 (x1 )

.. .

Q0 (xi )

Q1 (x) Q1 (x0 ) Q1 (x1 )

.. .

··· ··· ··· .. .

Qi (x) Qi (x0 ) ⎥ ⎟ ⎥ ⎟ Qi (x1 ) ⎥ /D ⎟ , i ⎥ ⎟ .

Q1 (xi )

···

Qi (xi )

..

⎤

⎞

⎦

⎠

(12)

with G−1,m (x) = Qm (x). From Lemma 2.1 we remark that Dn = V(µ+1)n+µ = (ci Qj (x))0≤i,j≤n is a strongly nonsingular matrix, so as we know it is also a block strongly nonsingular matrix [8]. Then the Schur complements are well defined. So using the Sylvester identity and some properties of the Schur complement, see [8], we get the Matrix Recursive Polynomial Interpolation Algorithm (MRPIA). Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

5

Algorithm 1 : the MRPIA

• P−1 (x) = 0; • for m = 0, 1, · · · , n 1. G−1,m (x) = Qm (x); 2. if m ≥ 1 for i = 0, 1, · · · , m − 1 [ ]−1 Gi,m (x) = Gi−1,m (x) − Gi−1,i (x) Gi−1,i (xi ) Gi−1,m (xi ); end i, [ ]−1 3. Pm (x) = Pm−1 (x) + Gm−1,m (x) Gm−1,m (xm ) (Ym − Pm−1 (xm ));

• end m.

3.2. Some properties of the MRPIA In this section we will recall some important properties of the MRPIA, which have been studied in [8], for a particular case µ = 1. The proofs of those properties are the same as in [8]. Proposition 3.1.

We have

(1) The MRPIA is well defined (i.e. Gm−1,m (xm ) are nonsingular). (2) cm (Gj,i (x)) = Gj,i (xm ) = 0(µ+1) , for i > j ≥ m. (3) cm (Pi (x)) = cm (Pn (x)) = Ym , for i = 0, . . . , n and m = 0, . . . , i. (4) Gj,i (x) is a matrix combination of the matrix polynomials Q0 (x), . . . , Qj (x), Qi (x). Now we will recall some other properties of the MRPIA, given in [8]. For m = 0, . . . , n, and for any real function g(x) we denote by cm ∗ g(x) = cm (g(x)) = g(xm ), we also set

[

Cm = c0 ∗

c1 ∗

[

Qm (x) = Q0 (x)

cm ∗ ,

]

··· Q1 (x)

···

Qm (x) ,

]

T T Sm = Qm (x)[Cm Qm (x)]−1 Cm .

Remark 3.1. Dm =

It is easy to see that, for m = 0, . . . , n, we have

,

T Cm Qm (x)

(13)

Pm (x) = Sm Pn (x),

(14)

Gm,i (x) = (I − Sm )Qi (x), Proposition 3.2. (1)

2 Sm

for i > m.

(15)

We have

= Sm ,

(2) Sm Si = Si Sm = Si , if m ≥ i. Now we set, for m = 0, . . . , n

[

Gm (x) = G−1,0 (x)

G0,1 (x)

···

Gm−1,m (x) ,

˜ m = [ci (Gj−1,j (x))]0≤i,j≤m = CmT Gm (x), D T T 1 T S˜ m = Gm (x)[Cm Gm (x)]−1 Cm = Gm (x)D˜ − m Cm . Proposition 3.3.

]

(16) (17) (18)

˜ n = CnT Gn (x) is a block lower triangular strongly nonsingular matrix and we have for m = 0, 1, . . . , n D

S˜ m = Sm ,

(19)

Pm (x) = S˜ m Pn (x),

(20)

Gm,i (x) = (I − S˜ m )Qi (x),

for i > m.

(21)

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

6

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Remark 3.2. Pn (x)

Using the relations (18) and (20), we get

=

S˜ n Pn (x)

=

1 T ˜− Gn (x)D n Cn Pn (x)

=

n ∑

Gi−1,i (x)ωi ,

i=0

˜ n is a block lower triangular matrix and ωi , for i = 0, 1, . . . , n, are (µ + 1) × 1 vectors solution of the block linear where D ˜ n Ωn = Y˜n , where system D

⎡ ⎤ ω0 ⎢ .. ⎥ Ωn = ⎣ . ⎦ , ωn and

⎡ ⎤ Y0

⎢ Y˜n = CnT Pn (x) = ⎣

.. ⎥ . .⎦

Yn 4. The general case We consider the general case of the Hermite polynomial interpolation problem, i.e. nm are different for some m. For this case the problem will be as follows: Find pN −1 ∈ PN −1 such that pN −1 (x) =

N −1 ∑

αi qi (x),

(22)

i=0

and for m = 0, 1, . . . , n, and k = 0, . . . , nm (k)

(k)

cm (pN −1 (x)) = pN −1 (xm ) = ym,k ,

(23)

where N =

n ∑

(nm + 1).

m=0

For this case we will apply the MRPIA, given by Algorithm 1, by choosing

µ = max{nm , m ∈ [0 : n]}. So for m = 0, 1, . . . , n we have nm ≤ µ. We set N′ =

n ∑

(µ + 1) = (µ + 1)(n + 1) = (µ + 1)n + µ + 1,

i=0

Iµ = {m : m ∈ [0 : n], nm < µ}, and SU = {zm,k : m ∈ Iµ , k = nm + 1, . . . , µ}, the set of unknowns zm,k to be determined, such that (k)

pN −1 (xm ) = zm,k ,

for m ∈ Iµ

and k = nm + 1, . . . , µ.

We also set N ′ −1

p˜ N ′ −1 (x) =

∑

βi qi (x),

i=0

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

7

So as xm , for m = 1, . . . , n, are distinct and using Lemma 2.1, we know that the interpolation polynomial p˜ N ′ −1 is unique and pN −1 is also unique then pN −1 (x) =

N −1 ∑

αi qi (x),

i=0 N ′ −1

p˜ N ′ −1 (x) =

∑

βi qi (x) =

i=0

N −1 ∑

N ′ −1

βi qi (x) +

∑

βi qi (x).

i=N

i=0

We must have pN −1 (x) = p˜ N ′ −1 (x), so we have to solve the following N ′ − N square linear system SS : βi = 0,

for i = N , N + 1, . . . , N ′ − 1.

The solutions of this system are the values of the unknowns zm,k of SU. So by replacing the values of zm,k in βi , for i = 0, 1, . . . , N − 1, we get p˜ N ′ −1 (x) =

N −1 ∑

βi qi (x) =

i=0

N −1 ∑

αi qi (x) = pN −1 (x).

i=0

5. Some examples We will give four examples for illustrating our algorithm for the two cases nm = µ and nm different for some m. For two examples we take the canonical basis of PN −1 , qi (x) = xi , and for the other two examples we choose qi (x) arbitrary. Example 5.1. This example is given for illustrating∑ the case where nm = µ. We choose n = 2, nm = µ = 2, then 8 i i N − 1 = (µ + 1)n + µ = 3 ∗ 2 + 2 = 8, and p8 (x) = i=0 αi x . We also choose qi (x) = x , for i = 0, 1, . . . , 8, and the set {(xm , ym,k ) : m = 0, 1, 2, k = 0, 1, 2} is given by m

xm 3

0

ym,0

2 −1 1

2

ym,2

1

1

3 0 1

3

4

−1

−

1

ym,1 1

−1 −

2

−1 1

Then applying the MRPIA, we get for p8 p8 (x)

=

1901 1944

x8 +

119951 62208

17069 3888

x3 +

x7 +

8927 6912

5743

x2 +

972

x6 +

2369 2304

x−

1763 7776

x5 −

145787 31104

x4 −

1783 2304

Example 5.2. This example is given for illustrating the case where ∑8 nm = µ. We choose n = 2, nm = µ = 2, for m = 0, 1, 2. Then N − 1 = (µ + 1)n + µ = 3 ∗ 2 + 2 = 8, and p8 (x) = i=0 αi qi (x), where the polynomials qi (x), for i = 0, 1, . . . , 8, are given by q0 (x) = 2 q3 (x) = x3 + 2x q6 (x) = x6 + 3x

q1 (x) = x − 2 q4 (x) = x4 + 2x − 4 q7 (x) = x7 + x2 + 2

q2 (x) = 2x2 + x q5 (x) = x5 + 2x q8 (x) = x8 + x3

and the set {(xm , ym,k ) : m = 0, 1, 2, k = 0, 1, 2} is given by m 0 1 2

xm 1 2 3

−

2 1

ym,0

ym,1

ym,2

1

0

−1

−1

2

3

−2

−1

1

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

8

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Then applying the MRPIA, we get for p8 p8 (x)

=

−

140526 3125

3317801 50000

108331

x8 −

x3 +

3125

x7 +

10989921 50000

657338 3125

x2 −

x6 +

6885053 50000

389939 6250

x+

x5 −

368399 1000

x4 +

1242073 50000

Example 5.3. This example is given ∑2 for illustrating the case where nm are different for some m. We choose n = 2, n0 = 0, n1 = 2, n2 = 1, then N − 1 = i=0 (ni + 1) − 1 = 5, µ = max(n0 , n1 , n2 ) = n1 = 2, Iµ = {m, : m ∈ [0 : n], nm < µ} = {0, 2}, N ′ − 1 = (µ + 1)n + µ = 8. The polynomials are the elements of the canonical basis qi (x) = xi , for i = 0, 1, . . . , 8. The set {(xm , ym,k ) : m = 0, 1, 2, k = 0, 1, 2} is given by m 0

xm 1

ym,0

ym,1

ym,2

3 2

−1

z0,1

z0,2

1

−

1

3

2

1

−

2

0

2 2

3

1

5

z2,2

4

∑8SU =i {z0,1 , z0,2 , z2,2 }, is the set of the unknowns to be determined. Then applying the MRPIA, we get p˜ 8 (x) = i=0 βi x and for SS ⎧ 3 243 1627393 441 ⎪ ⎪ z0,1 + z0,2 + z2,2 − =0 ⎪β6 = ⎪ 32 2 16 640 ⎪ ⎪ ⎪ ⎪ ⎨ 2349 117 2187 13683273 SS : β7 = − z0,1 − z0,2 − z2,2 + =0 ⎪ 256 128 128 5120 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 567 27 729 4336731 ⎪ ⎩β8 = z0,1 + z0,2 + z2,2 − =0 256

128

128

5120

The solution of SS is z0,1 =

16439 60

,

z0,2 = −

59953 30

,

z2,2 =

10457 90

.

Then we get

α0 =

1849

α3 =

9423

90

20

,

α1 = −

2119

,

α4 = −

10359

45

,

α2 = − ,

20

α5 =

10229 90

933 5

,

.

So p5 (x) is given as follows p5 (x) =

933 5

x5 −

10359 20

x4 +

9423 20

x3 −

10229 90

x2 −

2119 45

x+

1849 90

.

Example 5.4. This example is given ∑2 for illustrating the case where nm are different for some m. We choose n = 2, n0 = 1, n1 = 0, n2 = 2, then N − 1 = i=0 (ni + 1) − 1 = 5, µ = max(n0 , n1 , n2 ) = n2 = 2, Iµ = {m, : m ∈ [0 : n], nm < µ} = {0, 1}, N ′ − 1 = (µ + 1)n + µ = 8. The polynomials qi , for i = 0, 1, . . . , 8, are chosen as follows q0 (x) = 3 q3 (x) = x3 + x2 − 7x q6 (x) = x6 + x3 + x

q1 (x) = x − 8 q4 (x) = x4 − 2x2 + 2x q7 (x) = x7 + 2x4 − x2 + 6

q2 (x) = x2 + x − 3 q5 (x) = x5 − x3 + x + 2 q8 (x) = x8 + x5 + x2 − x

and the set {(xm , ym,k ) : m = 0, 1, 2, k = 0, 1, 2} is chosen as follows m 0 1 2

xm

−1 1 2 3 2

ym,0 1

ym,1 −1

ym,2 z0,2

−1

z1,1

z1,2

−2

−1

1

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.

A. Messaoudi, R. Sadaka and H. Sadok / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

SU = {z0,2 , and for SS

z1,1 ,

z1,2 }, is the set of the unknowns to be determined. Then applying the MRPIA, we get p˜ 8 (x) =

9

∑8

i=0

β i xi

⎧ 32 8 286708 104 ⎪ ⎪ z0,2 + z1,1 + z1,2 + =0 β6 = ⎪ ⎪ 3375 27 27 253125 ⎪ ⎪ ⎪ ⎪ ⎨ 128 4 10 50174 SS : β7 = − z0,2 + z1,1 + z1,2 + =0 ⎪ 3375 9 27 253125 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 8 4 52436 32 ⎪ ⎩β8 = z0,2 − z1,1 − z1,2 − =0 3375

The solution of SS is 1667 z0,2 = − , 450 Then we get

α0 = −

1217

α3 = −

1433

2500

,

27

z1,1 = −

α1 = − ,

11250 So p5 is given as follows p5 (x) =

674 5625

x5 −

1451 5625

3379 3750

766 625

α4 = −

27

,

,

z1,2 =

α2 =

1451 5625

x4 −

253125

,

1279 2500

α5 =

1433 11250

938 5625

,

674 5625

x3 +

.

.

1279 2500

x2 −

766 625

x−

1217 2500

.

6. Conclusion In this work we recalled the new formulation of the Hermite polynomial interpolation problem and gave the version of the Matrix Recursive Polynomial Interpolation Algorithm (MRPIA) for computing the Hermite polynomial algorithm for a particular case: the derivation orders are all equal to µ ≥ 0 (i.e., nm = µ ≥ 0, for m = 0, 1, . . . , n), some of its properties have also been given. For the general case, where the derivation orders, nm for m = 0, 1, . . . , n, are different for some m, we presented a method for solving this problem: we do the MRPIA for µ = max(n0 , n1 , . . . , nn ) and we solve a linear system. The study of the cost, the storage and the stability of this algorithm are under investigation and will be included in a forthcoming work. Acknowledgments We would like to thank the referees for their helpful comments and valuable suggestions. References [1] [2] [3] [4] [5] [6] [7] [8]

J. Stoer, R. Bulirsh, Introduction to Numerical Analysis, second ed., Springer-Verlag, New York, 1991. E. Süli, D. Mayers, An Introduction to Numerical Analysis, Cambridge University Press, 2003. A. Messaoudi, H. Sadok, Recursive polynomial interpolation algorithm (RPIA), Numer. Algorithms 76 (3) (2017) 675–694. C. Brezinski, Recursive interpolation, extrapolation and projection, J. Comput. Appl. Math. 9 (1983) 369–376. C. Brezinski, Other manifestations of the Schur complement, Linear Algebra Appl. 111 (1988) 231–247. R.W. Cottle, Manifestations of the Schur complement, Linear Algebra Appl. 8 (1974) 189–211. I. Schur, Potenzreihn im innern des einheitskreises, J. Reine Angew. Math. 147 (1917) 205–232. A. Messaoudi, R. Sadaka, H. Sadok, New algorithm for computing the Hermite interpolation polynomial, Numer. Algorithms 77 (4) (2017) 1069–1092. [9] K. Jbilou, A. Messaoudi, Matrix recursive interpolation algorithm for block linear systems, Direct methods, Linear Algebra Appl. 294 (1999) 137–154. [10] A. Messaoudi, M. Errachid, K. Jbilou, H. Sadok, GRPIA : a new algorithm for computing interpolation polynomials, Numer. Algorithms 80 (1) (2018) 253–278.

Please cite this article as: A. Messaoudi, R. Sadaka and H. Sadok, Matrix recursive polynomial interpolation algorithm: An algorithm for computing the interpolation polynomials, Journal of Computational and Applied Mathematics (2019) 112471, https://doi.org/10.1016/j.cam.2019.112471.