View of On some bivariate interpolation procedures

Rev. Anal. Num´er. Th´eor. Approx., vol. 33 (2004) no. 1, pp. 97–106 ictp.acad.ro/jnaat

ON SOME BIVARIATE INTERPOLATION PROCEDURES

DIMITRIE D. STANCU^∗and IOANA TAS¸CU^†

Abstract. In an important paper published in 1966 by the first author [10] was introduced and investigated a very general interpolation formula for univariate functions, which includes, as special cases, the classical interpolation formulae of Lagrange, Newton, Taylor and Hermite.

The purpose of the present paper is to extend that formula to the two- dimensional case. The remainders are expressed by means of partial divided differences and derivatives.

MSC 2000. 41A05, 41A10, 41A63.

Keywords. Bivariate interpolation, divided differences, bivariate polynomial interpolation formulas of Lagrange, Newton, Taylor, Hermite and Biermann.

1. INTRODUCTION

In this paper we start from a general decomposition formula of divided differences defined for a function f ∈C(D), D = [a, b]×[c, d], and for some groups of nodes from the rectangle D. We deduce a general interpolation formula for bivariate functions, corresponding to some general arrays of points.

As special cases we obtain several classical types of interpolation polynomials, including Lagrange, Newton, Biermann, Taylor and Hermite.

2. PRELIMINARIES

We first recall some of the principal results obtained in the paper [10].

Then we shall start from an array of M + 1 = p0+p1+· · ·+pm+m+ 1 points containing m+ 1 groups of nodes, denoted, by using subscripts and superscripts, by (a^k_i) i= 0, . . . , p_k,k= 0, . . . , m.

Let us use the following explicit notation for this array

(2.1) A=







a⁰₀ a⁰₁ . . . a⁰_p0 ... ... ... a^m₀ a^m₁ . . . a^m_pm





 .

We assume thata≤a^k₀ < a^k₁ <· · ·< a^k_p

k ≤b,k= 0, . . . , m.

The key role in deducing a general interpolation formula corresponding to the functionf ∈C[a, b] and the points (a^k_i),i= 0, . . . , p_k,k= 0, . . . , m, is the

∗“Babe¸s-Bolyai” University, Faculty of Mathematics and Computer Science, Str. M. Ko- g˘alniceanu 1, 400084 Cluj-Napoca, Romania, e-mail: [email protected].

†North University of Baia Mare, Victoriei 76, 430122, Baia Mare, Romania, e-mail:

[email protected].

following decomposition formula for divided differences on the distinct nodes a^k_i:

[a⁰₀, a⁰₁, . . . , a⁰_p0;a¹₀, a¹₁, . . . , a¹_p1;. . .;a^m₀ , a^m₁ , . . . , a^m_pm;f] = (2.2)

=

m

X

k=0

ha^k₀, a^k₁, . . . , a^k_p

k;_u^f(t)

k(t)

i,

where

uk(t) =γ0(x)γ1(x). . . γk−1(x)γk+1(x). . . γm(x), u0(x) = 1, and

γs(x) = (x−a^s₀)(x−a^s₁). . .(x−a^s_ps), s= 0, . . . , m.

The above brackets represent the symbol for divided differences.

If we introduce the node polynomial

(2.3) u(x) =

m

Y

s=0

(x−a^s₀)(x−a^s₁). . .(x−a^s_ps), then we can write: u_k(x) =u(x)/γ_k(x).

By using the decomposition formula (2.2) we can obtain the Stancu general interpolation formula

(2.4) f(x) = (SMf)(x) + (RMf)(x), where, in terms of divided differences, we have

(S_Mf)(x) =

m

X

k=0

u_k(x)

p_k

X

i=0

(x−a^k₀)(x−a^k₁). . .(x−a^k_i−1)(D^k_if), (2.5)

where

(D^k_if) =^ha^k₀, a^k₁, . . . , a^k_i;_u^f(t)

k(t)

i.

Obviously (SMf)(x) is a polynomial of degree not exceedingM =p0+p1+

· · ·+p_m+m.

The remainder of formula (2.4) has the following expression (R_Mf)(x) =u(x)(D_p0,p1,...,pmf)(x), where

(D_p0,p1,...,pmf)(x) = [x, x⁰₀, a⁰₁, . . . , a⁰_p0;. . .;a^m₀ , a^m₁ , . . . , a^m_pm;f(t)]

is the divided difference on all the points from the tableA and x.

Now let us present three remarkable special cases of the above approxima- tion formula.

(i) Ifp0 =p1 =· · ·=pm = 0 then we have a single column in the array (2.1) and the Stancu approximation formula (2.4) reduces to the Lagrange interpo- lation formula corresponding to the function f and the nodesa⁰₀, a¹₀, . . . , a^m₀ .

(ii) In the case m= 0 then the arrayA reduces to the nodes from the first row and we obtain the Newton interpolation formula corresponding to the nodesa⁰₀, a⁰₁, . . . , a⁰_p0 and the functionf.

(iii) If we assume that the nodes from the group a^k₀, a^k₁, . . . , a^k_pk tend to the same value b_k,k= 0, . . . , m, then the polynomial (2.5) becomes the Hermite osculatory interpolation polynomial under the form given in 1931 by P. Jo- hansen [3]:

(HMf)(x) =

m

X

k=0

u_k(x)

pk

X

j=0 (x−b_k)^j

j!

_f(t)

uk(t)

(j) t=bk

,

where we use the notations u(x) =

m

Y

k=0

(x−b_k)^pk+1, u_k(x) = _(x−b^u(x)

k)^pk+1.

It should be noticed that this polynomial can be written also under the more explicit form, given in 1948, by W. Simonsen [5]:

(HMf)(x) =

m

X

k=0 pk

X

j=0

h_k,j(x)f^(j)(b_k),

where the basic osculatory interpolation polynomials h_k,j(x) satisfy the rela- tions h^(s)_k,j(b_ν) = 0 (ν 6= k, s = 0, . . . , p_ν) and h^(s)_k,j(b_k) = δ_j^s, s = 0, . . . , p_k, whereδ^s_j is the Kronecker delta.

(iv) In the casep_k=k,k= 0, . . . , m, the tableAfrom (2.1) leads us to the following triangular array ofm(m+ 1)/2 base points

T =

















 a⁰₀ a¹₀ a¹₁ a²₀ a²₁ a²₂ ... ... ... . ..

a^m₀ a^m₁ a^m₂ . . . a^m_m

and in the Stancu interpolation polynomial (2.5) we have to replace p_k = k and uk(x) =u(x)/γk(x), where

u(x) =(x−a⁰₀(x−a¹₀)(x−a¹₁). . .(x−a^m₀ )(x−a^m₁ ). . .(x−a^m_m), u₀(x) = 1, γ_k(x) =(x−a^k₀)(x−a^k₁). . .(x−a^k_k).

3. SOME BIVARIATE INTERPOLATION FORMULAS

Let C(D) be the space of all real-valued functions continuous on the rec- tangleD= [a, b]×[c, d].

Besides the M+ 1 distinct points (a^k_i) from the interval [a, b], we consider also theN + 1 =s₀+s₁+· · ·+sn+n+ 1 distinct points (b^r_j), j= 0, . . . , sr, r = 0, . . . , n, from the interval [c, d], assuming that b^r₀ < b^r₁ < · · · < b^r_sr,

r= 0, . . . , n. We denote by B the table of these points

(3.1) B =























b⁰₀ b⁰₁ . . . b⁰_s0 b¹₀ b¹₁ . . . b¹_s1 ... ... ... b^r₀ b^r₁ . . . b^r_ss

... ... ... bⁿ₀ bⁿ₁ . . . bⁿ_sn





















 .

We want to approximate the function f ∈C(D) by an interpolation poly- nomial (T f)(x, y)(x, y) relative to the grid of nodes fromD:

G={(a^k_i, b^r_j), i= 0, . . . , p_k, j = 0, . . . , s_r, k= 0, . . . , m, r= 0, . . . , n}.

In order to find the expression of this bivariate interpolation polynomial using these nodes, we first apply formula presented at (2.4), with respect to the first variable and we obtain

(3.2) f(x, y) = (Sf)(x;y) + (Rf)(x;y), where

(Sf)(x;y) =

m

X

k=0

uk(x)

pk

X

i=0

ω0,i−1(x)(D^k_if)(y) + (Rf)(x;y), u_k(x) =γ₀(x)γ₁(x). . . γk−1(x)γ_k+1(x). . . γ_m(x), γ_k(x) =(x−a^k₀)(x−a^k₁). . .(x−a^k_p

k)

ω_0,i−1^k (x) =(x−a^k₀)(x−a^k₁). . .(x−a^k_i−1), ω^k_0,−1(x) = 1, (D^k_if)(y) =^ha^k₀, a^k₁, . . . , a^k_i;^f(t,y)_u

k(t)

i. The remainder has the following expression

(Rf)(x;y) =u(x)[x, x⁰₀, . . . , a⁰_p0;. . .;a^m₀ , . . . , a^m_pm;f(t, y)], where

u(x) =γ0(x)γ1(x). . . γm(x) =uk(x)γk(x).

Then we apply the above result with respect to the variableyand the points b^r_j from the array B given at (3.1). We obtain

(3.3) (T f)(x, y) =

m

X

k=0 n

X

r=0

u_k(x)v_r(y)

pk

X

i=0 sr

X

j=0

ω_0,i−1^k (x)γ_0,j−1^r (y)(D^k,r_i,jf), where

γ_0,j−1^r (y) =(y−b^r₀)(y−b^r₁). . .(y−b^r_j−1), γ_0,−1^r (y) = 1, vr(y) =δ0(y)δ1(y). . . δr−1(y)δr+1(y). . . δn(y), v0(y) = 1,

δt(y) =(y−b^t₀)(y−b^t₁). . .(y−b^t_st), k= 0, . . . , n, v(y) =vr(y)δr(y).

On the other hand we used the bidimensional divided difference D_i,j^k,rf =

"

a^k₀, a^k₁, . . . , a^k_i

b^r₀, b^r₁, . . . , b^r_j ; f(t, z) u_k(t)v_r(z)

# .

The remainder of the interpolation formula (3.2) has the following expres- sion

(Rf)(x, y) =

=u(x)[x, a⁰₀, a⁰₁, . . . , a⁰_p0;. . .;a^m₀ , a^m₁ , . . . a^m_pm;f(t, y)]

+v(y)

m

X

k=0

u_k(x)

pk

X

i=0

ω0,i−1(x)

a^k₀, a^k₁, . . . , a^k_i

y, b⁰₀, . . . , b⁰_s0, . . . , bⁿ₀, . . . , bⁿ_sn ;f(t, z) uk(t)

.

4. INTERPOLATION FORMULAS USING A RECTANGULAR OR A TRIANGULAR GRID OF NODES

A) In the special cases p₀ =p₁ =· · ·=p_m = 0, s₀ =s₁ =· · ·=s_m = 0 in the tables (2.1) and (3.1) remain only the first columns (a^k₀),k= 0, . . . , m, and (b^r₀), r = 0, . . . , n and the nodes will be the points M_k,r(a^k₀, b^r₀) which are at the intersections of the vertical linesx=a^k₀,k= 0, . . . , m, with the horizontal linesy =b^r₀,r= 0, . . . , n, in the plane.

In this case the interpolation polynomial (3.3) becomes (4.1) (Lm,nf)(x, y) =

m

X

k=0 n

X

r=0

uk(x)vr(y)

uk(a^k₀)vr(b^r₀)f(a^k₀, b^r₀), where

uk(x) =(x−a⁰₀)(x−a¹₀). . .(x−a^k−1₀ )(x−a^k+1₀ ). . .(x−a^m₀ ), v_r(y) =(y−b⁰₀)(y−b¹₀). . .(y−b^r−1₀ )(y−b^r+1₀ ). . .(y−bⁿ₀).

At (4.1) we have the bivariate Lagrange interpolation polynomial corre- sponding to the function f ∈C(D) and to the nodes M_k,r from the rectangle D= [a, b]×[c, d].

The corresponding remainder of the interpolation formula (4.2) f(x, y) = (L_m,nf)(x, y) + (R_m,nf)(x, y) can be expressed by means of the nodal polynomials

(4.3) u_m(x) =

m

Y

k=0

(x−a^k₀), v_n(y) =

n

Y

r=0

(y−b^r₀) and the partial divided differences, namely

(R_m,nf)(x, y) =u_m(x)[x, a⁰₀, a¹₀, . . . , a^m₀ ;f(t, y)]

+vn(y)[y, b⁰₀, b¹₀, . . . , bⁿ₀;f(x, z)]

−u_m(x)v_n(y)

x, a⁰₀, a¹₀, . . . , a^m₀

y, b⁰₀, b¹₀, . . . , bⁿ₀ ;f(t, z)

.

If we now assume that the function f has continuous partial derivatives f^(p,q)(x, y) on the rectangle D then this remainder can be expressed in the following form

(R_m,nf)(x, y) =_(m+1)!^um(x)f^(m+1,0)(ξ, y) +_(n+1)!^vn(y)f^(0,n+1)(x, η) (4.4)

−(m+1)!(n+1)!^um(x)vn(y) f^(m+1,n+1)(ξ, η),

whereξ∈(a, b) andη∈(c, d) are the same in both terms in which they occur.

B) If we assume that m=n= 0 then in the tables (2.1) and (3.1) remain only the first rows: a⁰₀, a⁰₁, . . . , a⁰_p0 andb⁰₀, b⁰₁, . . . , b⁰_s0 and the formula (3.3) leads us to the Newton bivariate interpolation formula

(4.5) f(x, y) = (N_p0,s0f)(x, y) + (R_p0,s0f)(x, y), where we have

(N_p0,s0f)(x, y) = (4.6)

=

p0

P

i=0 s0

P

j=0

(x−a⁰₀)(x−a⁰_n). . .(x−a⁰_i−1)(y−b⁰₀)(y−b⁰₁). . .(y−b⁰_j−1)(D⁰_i,jf)(t, z) and

(D⁰_i,jf)(x, y) =

a⁰₀, a⁰₁, . . . , a⁰_i

b⁰₀, b⁰₁, . . . , b⁰_j ;f(t, z)

is the bidimensional divided difference of the functionf on the indicated nodes.

The remainder of the interpolation formula (4.5) has the following expres- sion, in terms of partial divided differences,

(Rp0,s0f)(x, y) =up0(x)[x, a⁰₀, a⁰₁, . . . , a⁰_p0;f(t, y)]

+v_s0(y)[y, b⁰₀, b⁰₁, . . . , b⁰_s0;f(x, z)]

−up0(x)vs0(y)

x, a⁰₀, a⁰₁, . . . , a⁰_p0

y, b⁰₀, b⁰₁, . . . , b⁰_s0 ;f(t, z)

.

If f ∈ C^p0,s0(D) then we can obtain the following estimation for this re- mainder

(R_p0,s0f)(x, y) =_(p^up0(x)

0+1)!f^(p0+1,0)(ξ, y) +_(s^vs0(y)

0+1)!f^(0,s0+1)(x, η) (4.7)

−_(p^up0(x)vs0(y)

0+1)!(s0+1)!f^(p0+1,s0+1)(ξ, η), whereξ ∈(a, b) andη ∈(c, d).

C) If we use the notations p0 = p, a⁰_i = xi, b⁰_j = yj and assume that s₀ = p −i, i = 0, . . . , m; j = 0, . . . , n, then we arrive from (4.6) at the Biermann interpolation polynomial [1], [9]:

(4.8) (B_pf)(x, y) =

p

X

i=0 p−i

X

j=0

(x−x₀). . .(x−xi−1)(y−y₀). . .(y−yj−1)D_i,j(f),

where

D_i,j(f) =

x₀, x₁, . . . , x_i y0, y1, . . . , yj

;f(t, z)

.

The Biermann polynomial is of total (global) degree p inx and y and uses a triangular array of base nodes (xi, yi), i= 0, . . . , p,j= 0, . . . ,(p−i).

D) When the elements of the rows from the array (2.1) tend respectively to the same values, that is a^k_i → c_k, i = 0, . . . , p_k, k = 0, . . . , m, where c₀, c₁, . . . , c_m are distinct numbers, while the elements of the rows from the array (3.1) tend also to distinct values, that is b^r_j → dr, j = 0, . . . , sr, r = 0, . . . , n, then we can write the nodal polynomials

u(x) =

m

Y

i=0

(x−ci)^ri+1, v(y) =

n

Y

j=0

(y−dj)^sj+1 and

u_k(x) =u(x)/(x−c_k)^pk+1, v_r(y) =v(y)/(y−d_r)^sr+1. Because we have

ω_0,i−1^k (x) = (x−x_k)ⁱ, γ_j−1^r (y) = (y−d_r)^j,

"

a^k₀, a^k₁, . . . , a^k_i

b^r₀, b^r₁, . . . , b^r_j ; f(t, z) u_k(t)vr(z)

#

=

ck, ck, . . . , ck

d_r, d_r, . . . , d_r ; f(t, z) u_k(t)vr(z)

= 1 i!j!

f(t, z) u_k(t)v_r(z)

(i,j) ck,dr

,

it follows that the polynomial (3.3) may be expressed in the following form (4.9)

(HM,Nf)(x, y) =

m

X

k=0 n

X

r=0

u_k(x)vr(y)

pk

X

i=0 sr

X

j=0

(x−c_k)ⁱ(y−d_r)^j i!j!

_f(t,z)

uk(t)vr(z)

(i,j) ck,dr

,

which is the Hermite osculatory bivariate interpolation polynomial of degree (M, N), where M =p₀+p₁+· · ·+p_m+m and N =s₀+s₁+· · ·+s_n+n, which enjoys the following properties:

(H_M,N^(ν,µ))(c_k, d_r) =f^(ν,µ)(c_k, d_r),

whereν = 0, . . . , p_k,µ= 0, . . . , s_r,k= 0, . . . , mand r= 0, . . . , n.

The remainder of the approximation formula of the function f by this in- terpolation polynomial has the following expression

(R_M,Nf)(x, y) =

=u(x) x

1 , c0

p₀+ 1 , c1

p₁+ 1 , . . . , cm

p_m+ 1 ;f(t, y)

+v(y) y

1 , d0

s₀+ 1 , d1

s₁+ 1 , . . . , dn

s_n+ 1 ;f(x, z)

−u(x)v(y)·

· x

1, c₀

p0+ 1, c₁

p1+ 1, . . . , c_m

pm+ 1; d₀

s0+ 1, d₁

s1+ 1, . . . , d_n

sn+ 1;f(t, z)

In the brackets, in the first row, there are the coordinates of the nodes and in the second are indicated their corresponding orders of multiplicities.

Using the partial derivatives of the function f, this remainder can be esti- mated by the following formula

(RM,Nf)(x, y) =_(M+1)!^u(x) f^(M+1,0)(ξ, y) +_(N+1)!^v(y) f^(0,N+1)(x, η)

−(M+1)!(N+1)!^u(x)v(y) f^(M+1,N+1)(ξ, η), whereξ ∈(a, b),η∈(c, d).

Let us now mention that the Hermite interpolation polynomial (4.9) can be also expressed in a more explicit form

(H_M,Nf)(x, y) =

m

X

k=0 n

X

r=0 pk

X

i=0 sr

X

j=0

g_k,i(x)h_r,j(y)f^(i,j)(c_k, d_r), where

gk,i(x) = ^(x−c_i!^k)ih

pk−i

X

ν=0

(x−c_k)^ν ν!

1 u_k(t)

(ν) ck

i uk(x) and

h_r,j(y) = ^(y−d_j!^r)jh

sr−j

X

r=0

(y−dr)^µ µ!

1 vr(z)

(µ) dr

iv_r(y).

In the particular casem=n= 0,p0=p,s0 =swe obtain the Taylor-type bivariate formula

(4.10) f(x, y) =

p

X

k=0 s

X

r=0

(x−c0)^k(y−d0)^r

k!r! f^(k,r)(c₀, d₀) + (Rf)(x, y), where

(Rf)(x, y) =^(x−c_(p+1)!^0)p+1f^(p+1,0)(ξ, y) +^(y−d_(s+1)!^0)s+1f^(0,s+1)(x, η)

−^(x−c(p+1)!(s+1)!^{0)p+1(y−d0)s+1}f^(p+1,s+1)(ξ, η), ξ and η being the same in the terms in which they occur.

By using the integration by parts, the first author has obtained in [8] the following integral representation for the remainder of the Taylor-type for- mula (4.10)

(Rf)(x, y) = Z x

c0

(x−t)^p

p! f^(p+1,0)(t, y)dt+ Z y

d0

(y−z)^s

s! f^(0,s+1)(x, z)dz

− Z x

c0

Z y d0

(x−t)^p(y−z)^s

p!s! f^(p+1,s+1)(t, z)dtdz.

It should be further noted that employing the Biermann interpolation poly- nomial given at (4.8) we can obtain as a limit case the Taylor bivariate poly- nomial of total degreem:

(4.11) (Tpf)(x, y) =

p

X

i=0 p−i

X

j=0

(x−c)ⁱ(y−d)^j

i!j! f^(i,j)(c, d).

If we assume that f belongs to the class C^p+1 of functions having contin- uous all the partial derivatives of orders (p+ 1−i, i), (i= 0,1, . . . , p+ 1) in a neighborhood Ec,d of the point (c, d), then the remainder Rpf of the ap- proximation formula of f by the bivariate Taylor polynomial (4.11) can be represented under the following form

(Rpf)(x, y) =

= _p!¹ Z 1

0

(1−u)ⁿ(x−c)_∂x^∂ + (y−d)_∂y^∂(p+1)f c+ (x−c)u, d+ (y−d)udu, whenever the point (x, y) belongs to E_c,d. This formula was deduced in the paper [8] of the first author.

REFERENCES

[1] Biermann, O., Uber n¨¨ aherungsweise Kubaturen, Monatsh. Math. Phys., 14, pp. 211–225, 1903.

[2] Hermite, C.,Sur la formule d’interpolation de Lagrange, J. Reine Angew. Math.,84, pp. 70–79, 1878.

[3] Johansen, P.,Uber osculierende Interpolation, Skand. Actuarietidskr,¨ 14, pp. 231–237, 1931.

[4] Salzer, A.,A multi-point generalization of Newton’s divided difference formula, Proc.

Amer. Math. Soc.,13, pp. 210–212, 1962.

[5] Simonsen, W.,On divided differences and osculatory interpolation, Skand. Akturietid- skr.,31, pp. 157–164, 1948.

[6] Stancu, D. D.,On the interpolation formula of Hermite and on some applications of it(Romanian), Acad. R. P. Rom. Fil. Cluj, Stud. Cerc. Mat.,8, pp. 339–355, 1957.

[7] Stancu, D. D.,Considerations on the polynomial interpolation formulas for functions of several variables, Bul. Univ. Babe¸s-Bolyai, Ser. St. Nat. (Romanian), 1, pp. 43–82, 1957.

[8] Stancu, D. D.,On the integral representation of the remainder in a Taylor formula of two variables, Acad. R. P. Rom. Fil. Cluj, Stud. Cerc. Mat.,13, pp. 175–182, 1962.

[9] Stancu, D. D., The remainder of certain linear approximation formulas in two vari- ables, J. SIAM Numer. Anal. Ser. B,1, pp. 137–143, 1964.

[10] Stancu, D. D.,On Hermite’s osculatory interpolation formula and on some generali- zations of it, Mathematica (Cluj),8(31), pp. 373–391, 1966.

Received by the editors: January 15, 2004.