A New Algorithm to Approximate Bivariate Matrix Function via Newton-Thiele Type Formula

Abstract

A new method for computing the approximation of bivariate matrix function is introduced. It uses the construction of bivariate Newton-Thiele type matrix rational interpolants on a rectangular grid. The rational interpolant is of the form motivated by Tan and Fang (2000), which is combined by Newton interpolant and branched continued fractions, with scalar denominator. The matrix quotients are based on the generalized inverse for a matrix which is introduced by C. Gu the author of this paper, and it is effective in continued fraction interpolation. The algorithm and some other important conclusions such as divisibility and characterization are given. In the end, two examples are also given to show the effectiveness of the algorithm. The numerical results of the second example show that the algorithm of this paper is better than the method of Thieletype matrix-valued rational interpolant in Gu (1997).

1. Introduction

Matrix-valued rational interpolation and approximation theory have further practical application in many fields, such as in automatic control theory, computer science, and elementary particle physics [1]. Kuchminska et al. proved an analog of the van Vlerk theorem and constructed an interpolation formular of the Newton-Thiele type in [2]. Tan and Fang in [3] put emphasis on the study of Newton-Thiele bivariate rational interpolants. Graves-Morris provided a practical Thiele-fraction method for rational interpolation of vectors, based on the Samelson inverse in [4]. Gu et al. generalized the definition of Samelson inverse to the case of matrices and applied it to deal with the problems of rational interpolation of matrices [5–9]. Bose and Basu gave the existence, nonuniqueness, and recursive computation of the two-dimensional matrix Padé approximants in [10].

In this paper, we introduce a bivariate matrix-valued rational interpolant, with scalar denominator, in Newton-Thiele form motivated by [3]. The construction of the interpolant is combined by the classic methods: Newton interpolant and branched continued fractions. As we know, branched continued fraction has been studied by Cuyt and Verdonk and Siemaszko [11, 12] and many other authors. The interpolant of this paper is based on using the generalized inverse of matrices. A new algorithm to approximate bivariate matrix function is given and some examples are also put to show the effectiveness of the algorithm which is much better than the method mentioned in [9].

First, we will give the definition of the so-called generalized inverse of a matrix. Let 𝐶^𝑢×𝑣 consists of all 𝑢×𝑣 matrices with their elements in the complex plan 𝐶, and let 𝐴=(𝑎_𝑖𝑗),𝐵=(𝑏_𝑖𝑗),and𝐴,𝐵∈𝐶^𝑢×𝑣.

Definition 1 (see [8]). The scalar product of matrices 𝐴 and 𝐵 is defined by 𝐴⋅𝐵=^𝑢𝑣_𝑖=1𝑗=1𝑎_𝑖𝑗𝑏_𝑖𝑗=tr𝐴𝐵^∗, (1) where 𝐵^∗ denotes the transpose of 𝐵 and the Euclidean norm of 𝐴 is given by ‖𝐴‖=^𝑢𝑣_𝑖=1𝑗=1𝑎_𝑖𝑗^21/2. (2) It follows from Definition 1 and (2) that 𝐴⋅𝐴=^𝑢𝑣_𝑖=1𝑗=1𝑎_𝑖𝑗²=‖𝐴‖²=tr𝐴𝐴^∗, (3) where 𝐴 denotes the complex conjugate matrix of 𝐴 and 𝐴^∗ denotes the complex conjugate transpose matrix of 𝐴. On the basis of (2) and (3), the generalized inverse of the matrix 𝐴 is defined as 𝐴_𝑟⁻¹=1𝐴=𝐴‖𝐴‖²,𝐴≠0,𝐴∈𝐶^𝑢×𝑣, (4) in particular, for 𝐴∈𝑅^𝑢×𝑣,𝐴_𝑟⁻¹=1/𝐴=𝐴/‖𝐴‖²,𝐴≠0.

By means of generalized inverse 𝐴_𝑟⁻¹ for matrices, we want to define bivariate Newton-Thiele rational interpolants and give the algorithm and some properties on a rectangular grid.

2. Newton-Thiele Interpolation Formula

Let Λ_𝑚,𝑛={(𝑥_𝑙,𝑦_𝑘)∶𝑙=0,1,…,𝑚,𝑘=0,1,…,𝑛,(𝑥_𝑙,𝑦_𝑘)∈ℝ²}.

A matrix-valued set Π^𝑚,𝑛=𝑌_𝑙,𝑘∶𝑌_𝑙,𝑘=𝑌𝑥_𝑙,𝑦_𝑘∈𝐶^𝑢×𝑣,𝑥_𝑙,𝑦_𝑘∈Λ_𝑚,𝑛. (5)

We need to find a bivariate matrix rational function 𝑅_𝑚,𝑛=𝑁𝑥,𝑦𝑥,𝑦𝐷𝑥,𝑦, (6) whose numerator 𝑁(𝑥,𝑦) is a complex or real polynomial matrix and denominator 𝐷(𝑥,𝑦) is a real polynomial, and 𝑅_𝑚,𝑛𝑥_𝑙,𝑦_𝑘=𝑁𝑥_𝑙,𝑦_𝑘𝐷𝑥_𝑙,𝑦_𝑘=𝑌_𝑙,𝑘,𝑥_𝑙,𝑦_𝑘∈Λ_𝑚,𝑛. (7)

First, some notations need to be given as follows: 𝜑_0,0𝑥_𝑙,𝑦_𝑘=𝑌_𝑙,𝑘,∀𝑥_𝑙,𝑦_𝑘∈Λ_𝑚,𝑛,𝜑_1,0𝑥_𝑖,𝑥_𝑗,𝑦_𝑘=𝜑_0,0𝑥_𝑗,𝑦_𝑘−𝜑_0,0𝑥_𝑖,𝑦_𝑘𝑥_𝑗−𝑥_𝑖,𝜑_𝑙,0𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑘=𝜑_𝑙−1,0𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑗,𝑦_𝑘−𝜑_𝑙−1,0𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑦_𝑘×𝑥_𝑗−𝑥_𝑖⁻¹,𝜑_𝑙,1𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑟,𝑦_𝑠=𝑦_𝑠−𝑦_𝑟×𝜑_𝑙,0𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑠−𝜑_𝑙,0𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑟⁻¹,𝜑_𝑙,𝑘𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑎,…,𝑦_𝑏,𝑦_𝑟,𝑦_𝑠=𝑦_𝑠−𝑦_𝑟×𝜑_{𝑙,𝑘−1}𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑎,…,𝑦_𝑏,𝑦_𝑠−𝜑_{𝑙,𝑘−1}𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗,𝑦_𝑎,…,𝑦_𝑏,𝑦_𝑟⁻¹, (8) where the first subscript 𝑙 of 𝜑 means the number of nodes 𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗 minus 1, and the second subscript 𝑘 of 𝜑 means the number of nodes 𝑦_𝑎,…,𝑦_𝑏,𝑦_𝑟,𝑦_𝑠 minus 1.

For simplicity, we let the nodes 𝑥_𝑝,…,𝑥_𝑞,𝑥_𝑖,𝑥_𝑗 be replaced by 𝑥₀,…,𝑥_𝑙, and the nodes 𝑦_𝑎,…,𝑦_𝑏,𝑦_𝑟,𝑦_𝑠 be replaced by 𝑦₀,…,𝑦_𝑘, then we can get the following definition.

Definition 2. By means of generalized inverse (4), we define bivariate Newton-Thiele type matrix blending differences as follows: 𝜑_0,0𝑥_𝑙,𝑦_𝑘=𝑌_𝑙,𝑘,𝜑_1,0𝑥₀,𝑥₁,𝑦_𝑘=𝜑_0,0𝑥₁,𝑦_𝑘−𝜑_0,0𝑥₀,𝑦_𝑘𝑥₁−𝑥₀,𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦_𝑘=𝜑_𝑙−1,0𝑥₀,…,𝑥_𝑙−2,𝑥_𝑙,𝑦_𝑘−𝜑_𝑙−1,0𝑥₀,…,𝑥_𝑙−2,𝑥_𝑙−1,𝑦_𝑘×𝑥_𝑙−𝑥_𝑙−1⁻¹,𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁=𝑦₁−𝑦₀𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₁−𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀,𝜑_0,𝑘𝑥_𝑙,𝑦₀,…,𝑦_𝑘=𝑦_𝑘−𝑦_𝑘−1×𝜑_0,𝑘−1𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−𝜑_0,𝑘−1𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−1⁻¹,𝜑_𝑙,𝑘𝑥₀,𝑥₁,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘=𝑦_𝑘−𝑦_𝑘−1×𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−1⁻¹. (9) We assume that for all 𝑙,𝑥_𝑙≠𝑥_𝑙−1 and 𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘≠𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−1,∀𝑙,𝑘. (10)
From (9) matrix-valued Newton-Thiele type continued fractions for two-variable function can be constructed as follows: 𝑅_𝑚,𝑛𝑥,𝑦=𝑈₀𝑦+𝑈₁𝑦𝑥−𝑥₀+⋯+𝑈_𝑚𝑦𝑥−𝑥₀×𝑥−𝑥₁⋯𝑥−𝑥_𝑚−1, (11) where for 𝑙=0,1,…,𝑚𝑈_𝑙𝑦=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦−𝑦₀𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁+⋯+𝑦−𝑦_𝑛−1∣∣𝜑_𝑙,𝑛𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑛. (12)

Remark 3. Because of the construction, 𝐷(𝑥,𝑦) in 𝑅_𝑚,𝑛(𝑥,𝑦) mentioned as (6) is actually 𝐷(𝑦) and each 𝑈_𝑙(𝑦) depends on 𝑛.

We can also construct the antithetical form of bivariate matrix continued fraction in light of Definition 4.

Definition 4. By means of generalized inverse (4), we define bivariate antithetical Newton-Thiele type matrix blending differences 𝜓_0,0𝑥_𝑙,𝑦_𝑘=𝑌_𝑙,𝑘,𝜓_0,1𝑥_𝑙,𝑦₀,𝑦₁=𝜓_0,0𝑥_𝑙,𝑦₁−𝜓_0,0𝑥_𝑙,𝑦₀𝑦₁−𝑦₀𝜓_0,𝑘𝑥_𝑙,𝑦₀,…,𝑦_𝑘=𝜓_0,𝑘−1𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−𝜓_0,𝑘−1𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−1×𝑦_𝑘−𝑦_𝑘−1⁻¹,𝜓_1,𝑘𝑥₀,𝑥₁,𝑦₀,…,𝑦_𝑘=𝑥₁−𝑥₀𝜓_0,𝑘𝑥₁,𝑦₀,…,𝑦_𝑘−𝜓_0,𝑘𝑥₀,𝑦₀,…,𝑦_𝑘,𝜓_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦_𝑘=𝑥_𝑙−𝑥_𝑙−1×𝜓_𝑙−1,0𝑥₀,…,𝑥_𝑙−2,𝑥_𝑙,𝑦_𝑘−𝜓_𝑙−1,0𝑥₀,…,𝑥_𝑙−1,𝑦_𝑘⁻¹,𝜓_𝑙,𝑘𝑥₀,𝑥₁,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘=𝑥_𝑙−𝑥_𝑙−1×𝜓_{𝑙−1,𝑘}𝑥₀,…,𝑥_𝑙−2,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−𝜓_{𝑙−1,𝑘}𝑥₀,…,𝑥_𝑙−1,𝑦₀,…,𝑦_𝑘⁻¹. (13)
From (13) we can get the antithetical formula for two-variable function 𝑅_𝑚,𝑛𝑥,𝑦=𝑉₀(𝑥)+𝑉₁(𝑥)𝑦−𝑦₀+⋯+𝑉_𝑛(𝑥)𝑦−𝑦₀𝑦−𝑦₁⋯𝑦−𝑦_𝑛−1, (14) where for 𝑘=0,1,…,𝑛𝑉_𝑘(𝑥)=𝜓_0,𝑘𝑥₀,𝑦₀,…,𝑦_𝑘+𝑥−𝑥₀𝜓_1,𝑘𝑥₀,𝑥₁,𝑦₀,…,𝑦_𝑘+⋯+𝑥−𝑥_𝑚−1∣∣𝜓_𝑚,𝑘𝑥₀,…,𝑥_𝑚,𝑦₀,…,𝑦_𝑘. (15)

Now we will give two theorems about 𝑅_𝑚,𝑛(𝑥,𝑦) as in (11) and (12) and 𝑅_𝑚,𝑛(𝑥,𝑦) as in (14) and (15). First of all, some notations need to be defined.

In (12), for 𝑙=0,1,…,𝑚, let 𝑈_𝑙^(𝑠)𝑦=𝜑_𝑙,𝑠𝑥₀,𝑥₁,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑠+𝑦−𝑦_𝑠𝑈_𝑙^(𝑠+1)𝑦,𝑠=0,1,…,𝑛−1, (16) where 𝑈_𝑙^(𝑛)𝑦=𝜑_𝑙,𝑛𝑥₀,𝑥₁,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑛,𝑈_𝑙⁽⁰⁾𝑦=𝑈_𝑙𝑦. (17) Similarly, in (15), for 𝑘=0,1,…,𝑛, let 𝑉_𝑘^(𝑠)(𝑥)=𝜓_𝑠,𝑘𝑥₀,𝑥₁,…,𝑥_𝑠,𝑦₀,…,𝑦_𝑘+𝑥−𝑥_𝑠𝑉_𝑘^(𝑠+1)(𝑥),𝑠=0,1,…,𝑚−1, (18) where 𝑉_𝑘^(𝑚)(𝑥)=𝜓_𝑚,𝑘𝑥₀,𝑥₁,…,𝑥_𝑚,𝑦₀,…,𝑦_𝑘,𝑉_𝑘⁽⁰⁾(𝑥)=𝑉_𝑘(𝑥). (19)

Theorem 5. Let(i)𝜑_𝑙,𝑘(𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘),𝑙=0,1,…,𝑚;𝑘=0,1,…,𝑛 exist and be nonzero (except for  𝜑_𝑙,0(𝑥₀,…,𝑥_𝑙,𝑦₀)),(ii)𝑈_𝑙^(𝑠)(𝑦)=𝜑_𝑙,𝑠(𝑥₀,𝑥₁,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑠)+((𝑦−𝑦_𝑠)/𝑈_𝑙^(𝑠+1)(𝑦)) satisfy 𝑈_𝑙^(𝑠+1)(𝑦_𝑠)≠0,  𝑠=0,1,…,𝑛−1, then 𝑅_𝑚,𝑛(𝑥,𝑦) as in (11) and (12) exists such that 𝑅_𝑚,𝑛(𝑥_𝑙,𝑦_𝑘)=𝑌_𝑙𝑘,(𝑥_𝑙,𝑦_𝑘)∈Λ_𝑚,𝑛.

Proof. For simpleness, let 𝜑_𝑙,𝑘𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘=𝜑_𝑙,𝑘 (20) if the conditions hold, thus for 𝑙=0,1,…,𝑚;𝑘=0,1,…,𝑛𝑈_𝑙𝑦_𝑘=𝑈_𝑙⁽⁰⁾𝑦_𝑘=𝜑_𝑙,0+𝑦_𝑘−𝑦₀𝜑_𝑙,1𝑦+⋯+^𝑘−𝑦_𝑘−1∣∣𝜑_𝑙,𝑛+𝑦_𝑘−𝑦_𝑘∣∣𝑈_𝑙^(𝑘+1)𝑦_𝑘 (21)(𝑦_𝑘−𝑦_𝑘)/𝑈_𝑙^(𝑘+1)(𝑦_𝑘)=0, since 𝑈_𝑙^(𝑘+1)(𝑦_𝑘)≠0, then we get𝑈_𝑙𝑦_𝑘=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦_𝑘−𝑦₀𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁𝑦+⋯+^𝑘−𝑦_𝑘−2∣∣𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−1+𝑦_𝑘−𝑦_𝑘−1∣∣𝜑_𝑙,𝑘𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦_𝑘−𝑦₀𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁+𝑦+⋯_𝑘−𝑦_𝑘−2∣∣𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−1+𝑦_𝑘−𝑦_𝑘−1∣∣𝑦_𝑘−𝑦_𝑘−1×𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘−𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−1⁻¹=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦_𝑘−𝑦₀𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁𝑦+⋯+^𝑘−𝑦_𝑘−2∣∣𝜑_{𝑙,𝑘−1}𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘−2,𝑦_𝑘=⋯=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦_𝑘−𝑦₀𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦_𝑘=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦_𝑘. (22)We can use (9), (11), and (22) to find that 𝑅_𝑚,𝑛𝑥_𝑙,𝑦_𝑘=𝑈₀𝑦_𝑘+𝑈₁𝑦_𝑘𝑥_𝑙−𝑥₀+⋯+𝑈_𝑙𝑦_𝑘𝑥_𝑙−𝑥₀𝑥_𝑙−𝑥₁⋯𝑥_𝑙−𝑥_𝑙−1=𝑌𝑥_𝑙,𝑦_𝑘. (23)

We can similarly prove the following result.

Theorem 6. Let(i)𝜓_𝑙,𝑘(𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘),𝑙=0,1,…,𝑚;𝑘=0,1,…,𝑛 exist and be nonzero (except for  𝜓_0,𝑘(𝑥₀,𝑦₀,…,𝑦_𝑘)),(ii)𝑉_𝑘^(𝑠)(𝑥)=𝜓_𝑠,𝑘(𝑥₀,𝑥₁,…,𝑥_𝑠,𝑦₀,…,𝑦_𝑘)+((𝑥−𝑥_𝑠)/𝑉_𝑘^(𝑠+1)(𝑥)) satisfy 𝑉_𝑘^(𝑠+1)(𝑥_𝑠)≠0,  𝑠=0,1,…,𝑚−1, then 𝑅_𝑚,𝑛(𝑥,𝑦) as in (14) and (15) exists such that 𝑅_𝑚,𝑛(𝑥_𝑙,𝑦_𝑘)=𝑌_𝑙𝑘,(𝑥_𝑙,𝑦_𝑘)∈Λ_𝑚,𝑛.

Now we give the Newton-Thiele Matrix Interpolation Algorithm (NTMIA).

Algorithm 7 (NTMIA). Input:  {(𝑥_𝑖,𝑦_𝑗),𝑌_𝑖,𝑗}(𝑖=0,1,…,𝑚;𝑗=0,1,…,𝑛)

Output:  𝑅_𝑚(𝑥,𝑦).

Step  1. For all (𝑥_𝑖,𝑦_𝑗)∈Λ_𝑚,𝑛, let 𝜑(𝑥_𝑖,𝑦_𝑗)=𝑌_𝑖,𝑗.

Step  2. For 𝑗=0,1,…,𝑛;𝑝=1,2,…,𝑚;𝑖=𝑝,𝑝+1,…,𝑚, compute 𝜑_𝑝,0𝑥₀,…,𝑥_𝑝−1,𝑥_𝑖,𝑦_𝑗=𝜑_𝑝−1,0𝑥₀,…,𝑥_𝑝−2,𝑥_𝑖,𝑦_𝑗−𝜑_𝑝−1,0𝑥₀,…,𝑥_𝑝−1,𝑦_𝑗𝑥_𝑖−𝑥_𝑝−1. (24)

Step  3. For 𝑖=0,1,…,𝑚;𝑞=1,2,…,𝑛;𝑗=𝑞,𝑞+1,…,𝑛, compute 𝜑_𝑖,𝑞𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−1,𝑦_𝑗=𝑦_𝑗−𝑦_𝑞−1×𝜑_{𝑖,𝑞−1}𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−2,𝑦_𝑗−𝜑_{𝑖,𝑞−1}𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−1⁻¹=𝑦_𝑗−𝑦_𝑞−1𝜑_{𝑖,𝑞−1}𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−2,𝑦_𝑗−𝜑_{𝑖,𝑞−1}𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−1×𝜑𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−2,𝑦_𝑗−𝜑𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑞−1²⁻¹. (25)

Step  4. For 𝑖=0,1,…,𝑚;𝑗=0,1,…,𝑛, compute 𝑎_𝑖,𝑗=𝜑_𝑖,𝑗(𝑥₀,…,𝑥_𝑖,𝑦₀,…,𝑦_𝑗).

Step  5. For 𝑖=0,1,…,𝑚;𝑗=1,…,𝑛 judge if 𝑎_𝑖,𝑗≠0; if yes go to Step 6, otherwise exist and show “the algorithm is not valid.”

Step  6. For 𝑖=0,1,…,𝑚;𝑗=𝑛−2,𝑛−3,…,0𝑄_𝑖,𝑛=1,𝑃_𝑖,𝑛=𝑎_𝑖,𝑛,𝑄_{𝑖,𝑛−1}=𝑎_𝑖,𝑛²,𝑃_{𝑖,𝑛−1}=𝑎_{𝑖,𝑛−1}𝑄_{𝑖,𝑛−1}+𝑦−𝑦_𝑛−1𝑃_𝑖,𝑛,𝑄_𝑖,𝑗=𝑎_𝑖,𝑗+1²𝑄_𝑖,𝑗+1+2𝑦−𝑦_𝑗+1tr𝑎_𝑖,𝑗+1𝑃^∗_𝑖,𝑗+2+𝑦−𝑦_𝑗+1²𝑄_𝑖,𝑗+2,𝑃_𝑖,𝑗=𝑎_𝑖,𝑗𝑄_𝑖,𝑗+𝑦−𝑦_𝑗𝑃_𝑖,𝑗+1. (26)

Step  7. For 𝑖=0,1,…,𝑚;𝑗=𝑛,𝑛−1,…,1, let 𝐵_𝑖,𝑗(𝑦)=𝑃_𝑖,𝑗(𝑦)/𝑄_𝑖,𝑗(𝑦); if 𝐵_𝑖,𝑗(𝑦_𝑗−1)≠0 then go to Step 8, otherwise exist and show “the algorithm is not valid.”

Step  8. For 𝑖=0,1,…,𝑚, let 𝐵_𝑖(𝑦)=𝐵_𝑖,0(𝑦)=𝑃_𝑖,0(𝑦)/𝑄_𝑖,0(𝑦).

Step  9. For 𝑘=1,…,𝑚 let 𝑅₀(𝑥,𝑦)=𝐵₀(𝑦), then compute 𝑅_𝑘(𝑥,𝑦)=𝑅_𝑘−1(𝑥,𝑦)+⋯+(𝑥−𝑥₀)⋯(𝑥−𝑥_𝑘−1)𝐵_𝑘(𝑦).

3. Some Properties

Lemma 8 (see [5]). Define 𝑈₀(𝑦)=𝜑_0,0(𝑥₀,𝑦₀)+(𝑦−𝑦₀)/𝜑_0,1(𝑥₀,𝑦₀,𝑦₁)+⋯+(𝑦−𝑦_𝑛−1)/𝜑_0,𝑛(𝑥₀,𝑦₀,…,𝑦_𝑛); if we use generalized inverse by a tail-to-head rationalization, then a polynomial matrix 𝐸₀(𝑦) and a real polynomial 𝐹₀(𝑦) exist such that(i)𝑈₀(𝑦)=𝐸₀(𝑦)/𝐹₀(𝑦),𝐹₀(𝑦)≥0, (ii)𝐹₀(𝑦)∣‖𝐸₀(𝑦)‖²,“∣” means the sign of divisibility.

Definition 9. 𝑅(𝑥,𝑦)=𝑁(𝑥,𝑦)/𝐷(𝑥,𝑦) is defined to be of type [𝑙/𝑡] if deg{𝑁_𝑖𝑗}≤𝑙 for 1≤𝑖≤𝑢,1≤𝑗≤𝑣,deg{𝑁_𝑖𝑗}=𝑙 for some (𝑖,𝑗) and deg{𝐷}=𝑡, where 𝑁(𝑥,𝑦)=(𝑁_𝑖𝑗(𝑥,𝑦))∈𝐶^𝑢×𝑣.

Lemma 10 (see [5]). Let 𝑈₀(𝑦) be defined as in Lemma 8, and 𝑈₀(𝑦)=𝐸₀(𝑦)/𝐹₀(𝑦); if 𝑛 is even, 𝑈₀(𝑦) is of type [𝑛/𝑛]; if 𝑛 is odd, 𝑈₀(𝑦) is of type [𝑛/𝑛−1].

Theorem 11. Let 𝜑_𝑙,𝑘∈𝐶^𝑢×𝑣, (𝑥_𝑙,𝑦_𝑘)∈Λ_𝑚,𝑛, and 𝑥,𝑦∈𝑅. Define 𝑈_𝑙(𝑦)𝑜𝑓𝑅_𝑚𝑛(𝑥,𝑦) as in (7) by a tail-to-head rationalization using generalized inverse and suppose every intermediate denominator be nonzero in the operation, then a square polynomial matrix 𝑁(𝑥,𝑦) and a real polynomial 𝐷(𝑥,𝑦) exist such that(i)𝑅_𝑚,𝑛(𝑥,𝑦)=𝑁(𝑥,𝑦)/𝐷(𝑥,𝑦), (ii)𝐷(𝑥,𝑦)≥0, (iii)𝐷(𝑥,𝑦)∣‖𝑁(𝑥,𝑦)‖².

Proof. Consider the following algorithm for the construction of 𝑁(𝑥,𝑦), 𝐷(𝑥,𝑦), and 𝑅_𝑚𝑛(𝑥,𝑦).
Initialization: let 𝐷₀(𝑥,𝑦)=1, and 𝑁₀𝑥,𝑦=𝑈₀𝑦=𝜑_0,0+𝑦−𝑦₀𝜑_0,1+⋯+𝑦−𝑦_𝑛−1∣∣𝜑_0,𝑛. (27) By Lemma 8, a polynomial matrix 𝐸₀(𝑦) and a real polynomial 𝐹₀(𝑦) exist such that(a1)𝑁₀(𝑥,𝑦)=𝐸₀(𝑦)/𝐹₀(𝑦), (a2)𝐹₀(𝑦)≥0, (a3)𝐹₀(𝑦)∣‖𝐸₀(𝑦)‖².
Recursion: for 𝑗=1,2,…,𝑚 let 𝑆^(𝑗)𝑥,𝑦=𝑆^(𝑗−1)𝑥,𝑦+𝑈_𝑗𝑦𝑥−𝑥₀⋯𝑥−𝑥_𝑗−1 (28) with 𝑈_𝑗(𝑦)=𝐸_𝑗(𝑦)/𝐹_𝑗(𝑦),𝐹_𝑗(𝑦)∣‖𝐸_𝑗(𝑦)‖² having the representation(b1)𝑆^(𝑗)(𝑥,𝑦)=𝑁_𝑗(𝑥,𝑦)/𝐷_𝑗(𝑥,𝑦), (b2)𝐷_𝑗(𝑥,𝑦)≥0, (b3)𝐷_𝑗(𝑥,𝑦)∣‖𝑁_𝑗(𝑥,𝑦)‖², where 𝐷_𝑗(𝑥,𝑦) is a real polynomial.
Then we can get 𝑆^(𝑗+1)𝑥,𝑦=𝑆^(𝑗)𝑥,𝑦+𝑈_𝑗+1𝑥−𝑥₀⋯𝑥−𝑥_𝑗=𝑁_𝑗𝑥,𝑦𝐷_𝑗+𝐸𝑥,𝑦^𝑗+1𝑦𝐹_𝑗+1𝑦𝑥−𝑥₀⋯𝑥−𝑥_𝑗=𝑁_𝑗𝑥,𝑦𝐹_𝑗+1𝑦+𝐷_𝑗𝑥,𝑦𝐸_𝑗+1𝑦×𝑥−𝑥₀⋯𝑥−𝑥_𝑗×𝐷_𝑗𝑥,𝑦𝐹_𝑗+1𝑦⁻¹=𝑁_𝑗+1𝑥,𝑦𝐷_𝑗+1,𝑥,𝑦 (29) where 𝑁_𝑗+1𝑥,𝑦=𝑁_𝑗𝑥,𝑦𝐹_𝑗+1𝑦+𝐷_𝑗𝑥,𝑦𝐸_𝑗+1𝑦𝑥−𝑥₀⋯𝑥−𝑥_𝑗,𝐷_𝑗+1𝑥,𝑦=𝐷_𝑗𝑥,𝑦𝐹_𝑗+1𝑦. (30) Obviously, 𝑁_𝑗+1²=𝑁_𝑗(𝑥,𝑦)²𝐹²_(𝑗+1)𝑦+𝐷²_𝑗×𝐸𝑥,𝑦_(𝑗+1)(𝑦)²𝑥−𝑥₀²⋯𝑥−𝑥_𝑗²+𝑁_𝑗𝑥,𝑦⋅𝐸_𝑗+1𝑦𝐹_𝑗+1𝑦𝐷_𝑗×𝑥,𝑦𝑥−𝑥₀⋯𝑥−𝑥_𝑗+𝑁_𝑗𝑥,𝑦⋅𝐸_𝑗+1𝑦𝐹_𝑗+1𝑦𝐷_𝑗×𝑥,𝑦𝑥−𝑥₀⋯𝑥−𝑥_𝑗 (31) since 𝐷_𝑗𝑥,𝑦∣𝑁_𝑗𝑥,𝑦²,𝐹_𝑗+1𝑦∣𝐸_𝑗+1𝑦²,𝐷_𝑗+1𝑥,𝑦∣𝑁_𝑗+1𝑥,𝑦². (32) Termination: 𝑆^(𝑚)(𝑥,𝑦)=𝑅_𝑚,𝑛(𝑥,𝑦).

Theorem 12 (Characterization). Let 𝑅_𝑚𝑛(𝑥,𝑦)=𝑁(𝑥,𝑦)/𝐷(𝑥,𝑦) be expressed as 𝑅=𝑅_𝑚,𝑛𝑥,𝑦=𝑈₀𝑦+𝑈₁𝑦𝑥−𝑥₀+⋯+𝑈_𝑚𝑦𝑥−𝑥₀⋯𝑥−𝑥_𝑚−1, (33) where 𝑈_𝑙(𝑦)=𝐸_𝑖/𝐹_𝑖, as in (12) we assume 𝐹_𝑖,𝐹_𝑗 have no common factor, 𝑖,𝑗=0,1,…,𝑚,𝑖≠𝑗, then(i)if n is even, 𝑅 is of type [𝑚𝑛+𝑚+𝑛/𝑚𝑛+𝑛];(ii)if n is odd, 𝑅 is of type [𝑚𝑛+𝑛/𝑚𝑛+𝑛−𝑚−1].

Proof. The proof is recursive. Let 𝑛 be even.
For 𝑚=0,𝑅_0,𝑛(𝑥,𝑦)=𝑈₀(𝑦)=𝐸₀/𝐹₀, by Lemma 10, we find that it is of type [𝑛/𝑛].
For 𝑚=1, because deg{𝐸₀}=deg{𝐹₀}=𝑛,deg{𝐸₁}=deg{𝐹₁}=𝑛𝑅_1,𝑛=𝑁𝑥,𝑦¹𝐷₁=𝑈₀𝑦+𝑈₁𝑦𝑥−𝑥₀=𝐸₀𝐹₀+𝐸₁𝐹₁𝑥−𝑥₀=𝐸₀𝐹₁+𝐸₁𝐹₀𝑥−𝑥₀𝐹₀𝐹₁ (34) is of type [2𝑛+1/2𝑛].
For 𝑚=2, from the formula above we know that deg{𝑁₁}=2𝑛+1,deg{𝐷₁}=2𝑛,deg{𝐸₂}=𝑛,anddeg{𝐹₂}=𝑛. one has 𝑅_2,𝑛=𝑁𝑥,𝑦²𝐷₂=𝑁₁𝐷₁+𝐸₂𝐹₂𝑥−𝑥₀𝑥−𝑥₁=𝑁₁𝐹₂+𝐸₂𝐷₁𝑥−𝑥₀𝑥−𝑥₁𝐷₁𝐹₂. (35) It is easy to find that 𝑅_2,𝑛(𝑥,𝑦) is of type [3𝑛+2/3𝑛].
For 𝑚=𝑘, let 𝑅_𝑘,𝑛(𝑥,𝑦)=𝑁_𝑘/𝐷_𝑘 be of type [(𝑘+1)𝑛+𝑘/(𝑘+1)𝑛].
For 𝑚=𝑘+1, we consider that 𝑅_𝑘+1,𝑛=𝑁𝑥,𝑦^𝑘𝐷_𝑘+𝐸_𝑘+1𝐹_𝑘+1𝑥−𝑥₀⋯𝑥−𝑥_𝑘=𝑁_𝑘𝐹_𝑘+1+𝐸_𝑘+1𝐷_𝑘𝑥−𝑥₀⋯𝑥−𝑥_𝑘𝐷_𝑘𝐹_𝑘=𝑁_𝑘+1𝐷_𝑘+1, (36) where deg𝑁_𝑘+1=()=()(),𝐷𝑘+1)(𝑛+1+𝑛𝑘+2𝑛+𝑘+1deg_𝑘+1=(𝑘+1)𝑛+𝑛=(𝑘+2)𝑛. (37) Therefore, when 𝑛 is even, 𝑅 is of type [𝑚𝑛+𝑚+𝑛/𝑚𝑛+𝑛].
Similarly, if 𝑛 is odd, for 𝑚=0,𝑅_0,𝑛(𝑥,𝑦)=𝑈₀(𝑦)=𝐸₀/𝐹₀, by Lemma 10, we find that it is of type [𝑛/𝑛−1].
For 𝑚=1, because deg{𝐸₀}=deg{𝐸₁}=𝑛,deg{𝐹₀}=deg{𝐹₁}=𝑛−1𝑅_1,𝑛=𝑁𝑥,𝑦¹𝐷₁=𝑈₀𝑦+𝑈₁𝑦𝑥−𝑥₀=𝐸₀𝐹₀+𝐸₁𝐹₁𝑥−𝑥₀=𝐸₀𝐹₁+𝐸₁𝐹₀𝑥−𝑥₀𝐹₀𝐹₁ (38) is of type [2𝑛/2𝑛−2].
For 𝑚=2, from the formula above we know that deg{𝑁₁}=2𝑛,deg{𝐷₁}=2𝑛−2,deg{𝐸₂}=𝑛,anddeg{𝐹₂}=𝑛−1. one has 𝑅_2,𝑛=𝑁𝑥,𝑦²𝐷₂=𝑁₁𝐷₁+𝐸₂𝐹₂𝑥−𝑥₀𝑥−𝑥₁=𝑁₁𝐹₂+𝐸₂𝐷₁𝑥−𝑥₀𝑥−𝑥₁𝐷₁𝐹₂. (39) We get that 𝑅_2,𝑛(𝑥,𝑦) is of type [3𝑛/3𝑛−3].
For 𝑚=𝑘, let 𝑅_𝑘,𝑛=𝑁𝑥,𝑦^𝑘𝐷_𝑘 (40) be of type [(𝑘+1)𝑛/(𝑘+1)(𝑛−1)].
For 𝑚=𝑘+1, we consider that 𝑅_𝑘+1,𝑛=𝑁𝑥,𝑦^𝑘𝐷_𝑘+𝐸_𝑘+1𝐹_𝑘+1𝑥−𝑥₀⋯𝑥−𝑥_𝑘=𝑁_𝑘𝐹_𝑘+1+𝐸_𝑘+1𝐷_𝑘𝑥−𝑥₀⋯𝑥−𝑥_𝑘𝐷_𝑘𝐹_𝑘=𝑁_𝑘+1𝐷_𝑘+1, (41) where deg𝑁_𝑘+1=()()𝐷𝑘+1)(𝑛−1+𝑘+1=𝑘+2𝑛,deg_𝑘+1=(𝑘+1)(𝑛−1)+(𝑛−1)=(𝑘+2)(𝑛−1). (42)
Therefore, when 𝑛 is odd, 𝑅 is of type [(𝑚+1)𝑛/(𝑚+1)(𝑛−1)].

Definition 13. A matrix-valued Newton-Thiele type rational fraction 𝑅_𝑚,𝑛(𝑥,𝑦)=𝑁(𝑥,𝑦)/𝐷(𝑥,𝑦) is defined to be a bivariate generalized inverse and rational interpolant (BGIRI_NT) on the rectangular grid Λ_𝑚,𝑛 if(i)𝑅_𝑚𝑛(𝑥_𝑙,𝑦_𝑘)=𝑌_𝑙𝑘,(𝑥_𝑙,𝑦_𝑘)∈Λ_𝑚,𝑛, (ii)𝐷(𝑥_𝑙,𝑦_𝑘)≠0,(𝑥_𝑙,𝑦_𝑘)∈Λ_𝑚,𝑛, (iii)(a) if 𝑛 is even, deg𝑁𝑥,𝑦=𝑚𝑛+𝑚+𝑛,deg𝐷𝑥,𝑦=𝑚𝑛+𝑛, (43) (b) if n is odd, deg𝑁𝑥,𝑦=𝑚𝑛+𝑛,deg𝐷𝑥,𝑦=𝑚𝑛+𝑛−𝑚−1, (44)(iv)𝐷(𝑥,𝑦)∣‖𝑁(𝑥,𝑦)‖², (v)𝐷(𝑥,𝑦) is real, and 𝐷(𝑥,𝑦)≥0.
Now let us turn to the error estimate of the BGIRI_NT. Firstly, we give Lemma 14.

Lemma 14. Let 𝑅^∗_𝑚,𝑛𝑥,𝑦=𝑈^∗₀𝑦+𝑈^∗₁𝑦𝑥−𝑥₀+⋯+𝑈^∗_𝑚𝑦𝑥−𝑥₀𝑥−𝑥₁⋯𝑥−𝑥_𝑚−1, (45) where for 𝑙=0,1,…,𝑚𝑈^∗_𝑙𝑦=𝜑_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦−𝑦₀𝜑_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁+⋯+𝑦−𝑦_𝑛−1∣∣𝜑_𝑙,𝑛𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑛+𝑦−𝑦_𝑛∣∣𝜑_𝑙,𝑛+1𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑛,,𝑦 (46) then 𝑅^∗_𝑚,𝑛(𝑥,𝑦) satisfies 𝑅^∗_𝑚,𝑛𝑥_𝑖,𝑦=𝑓𝑥_𝑖𝑅,𝑦,𝑖=0,1,…,𝑚,_∗𝑚,𝑛𝑥,𝑦_𝑗=𝑅_𝑚,𝑛𝑥,𝑦_𝑗,𝑗=0,1,…,𝑛. (47)

Remark 15. We delete the proof of the above lemma since it can be easily generalized from [3].

Now the error estimate Theorem 16 will be put which is also motivated by [3]. We let 𝑅_𝑚,𝑛=𝑁𝑥,𝑦𝑥,𝑦𝐷𝑥,𝑦,𝑅^∗_𝑚,𝑛=𝑁𝑥,𝑦^∗𝑥,𝑦𝐷^∗.𝑥,𝑦 (48)

Theorem 16. Let a matrix 𝑓(𝑥,𝑦) be max(𝑚+1,𝑛+1) times continuously differentiable on a set 𝐷={𝑎≤𝑥≤𝑏,𝑐≤𝑦≤𝑑} containing the points Λ_𝑚,𝑛={(𝑥_𝑙,𝑦_𝑘)∶𝑙=0,1,…,𝑚,𝑘=0,1,…,𝑛,(𝑥_𝑙,𝑦_𝑘)∈ℝ²}, then, for any point (𝑥,𝑦)∈𝐷, there exists a point (𝜉,𝜂)∈𝐷 such that 𝑓𝑥,𝑦−𝑅_𝑚,𝑛=1𝑥,𝑦𝐷𝑥,𝑦𝐷^∗×𝜔𝑥,𝑦_𝑚+1(𝑥)()!𝜕𝑚+1𝜕𝑥^𝑚+1𝑚+1𝐸₁+𝜛𝜉,𝑦^𝑛+1𝑦()!𝜕𝑛+1𝜕𝑦^𝑛+1𝑛+1𝐸₂𝑥,𝜂, (49) where 𝜔_𝑚+1(𝑥)=𝑥−𝑥₀𝑥−𝑥₁⋯𝑥−𝑥_𝑚,𝜛_𝑛+1𝑦=𝑦−𝑦₀𝑦−𝑦₁⋯𝑦−𝑦_𝑛,𝐸₁𝑥,𝑦=𝐷𝑥,𝑦𝐷^∗𝑥,𝑦𝑓𝑥,𝑦−𝑅^∗_𝑚,𝑛,𝐸𝑥,𝑦₂𝑥,𝑦=𝐷𝑥,𝑦𝐷^∗𝑥,𝑦𝑅^∗_𝑚,𝑛𝑥,𝑦−𝑅_𝑚,𝑛𝑥,𝑦. (50)

Proof. Let 𝐸₁𝑥,𝑦=𝐷𝑥,𝑦𝐷^∗𝑥,𝑦𝑓𝑥,𝑦−𝑅^∗_𝑚,𝑛,𝐸𝑥,𝑦₂𝑥,𝑦=𝐷𝑥,𝑦𝐷^∗𝑥,𝑦𝑅^∗_𝑚,𝑛𝑥,𝑦−𝑅_𝑚,𝑛,𝐸𝑥,𝑦𝑥,𝑦=𝐸₁𝑥,𝑦+𝐸₂𝑥,𝑦. (51) From Lemma 14, we know 𝐸₁𝑥_𝑖,𝑦=0,𝑖=0,1,…,𝑚, (52) which results in 𝐸₁𝜔(𝑥,𝑦)=^𝑚+1(𝑥)𝜕(𝑚+1)!𝜕𝑥^𝑚+1𝑚+1𝐸₁(𝑥,𝑦)_𝑥=𝜉, (53) where 𝜔_𝑚+1(𝑥)=(𝑥−𝑥₀)(𝑥−𝑥₁)⋯(𝑥−𝑥_𝑚), and 𝜉 is a number contained in the interval (𝑎,𝑏) which may depend on 𝑦. Similarly, from 𝐸₂𝑥,𝑦_𝑗=0,𝑗=0,1,…,𝑛, (54) we can get that 𝐸₂𝜛(𝑥,𝑦)=^𝑛+1𝑦()!𝜕𝑛+1𝜕𝑦^𝑛+1𝑛+1𝐸₂𝑥,𝑦_𝑦=𝜂, (55) where 𝜛_𝑛+1(𝑦)=(𝑦−𝑦₀)(𝑦−𝑦₁)⋯(𝑦−𝑦_𝑛), and 𝜂 is a number contained in the interval (𝑐,𝑑) which may depend on 𝑥.
Therefore, 𝐸𝑥,𝑦=𝐸₁𝑥,𝑦+𝐸₂=𝜔𝑥,𝑦_𝑚+1(𝑥)()!𝜕𝑚+1𝜕𝑥^𝑚+1𝑚+1𝐸₁𝑥,𝑦_𝑥=𝜉+𝜛_𝑛+1𝑦()!𝜕𝑛+1𝜕𝑦^𝑛+1𝑛+1𝐸₂𝑥,𝑦_𝑦=𝜂. (56) The proof is thus completed.

4. Numerical Examples

Example 17. Let (𝑥_𝑖,𝑦_𝑗) and 𝑌_𝑖,𝑗,𝑖,𝑗=0,1,2 be given in Table 1: 𝑅_2,2𝑥,𝑦=𝑈₀𝑦+𝑈₁𝑦𝑥+𝑈₂𝑦𝑥(𝑥−1). (57) Using Algorithm 7, we can get 𝑈₀𝑦=+10𝑦00⁰⁰⁰¹+𝑦−1∣∣⁰⁰_−2/5−1/5,𝑈₁𝑦=+0−1𝑦00^01/2−1/20+𝑦−1∣∣⁰¹₁₀,𝑈₂𝑦=10120+_0−3/𝑦_51/50+𝑦−1∣∣^0−15/14_{−5/14−5/7}. (58) Based on generalized inverse (4), we obtain 𝑅_2,2=𝑎𝑥,𝑦₁₁𝑎₁₂𝑎₂₁𝑎₂₂×70𝑦⁶−368𝑦⁵+892𝑦⁴−1240𝑦³+1040𝑦²−512𝑦+128⁻¹=𝑁𝑥,𝑦𝐷,𝑥,𝑦 (59) where 𝑎₁₁=70𝑦⁶−368𝑦⁵+892𝑦⁴−1240𝑦³+1040𝑦²𝑎−512𝑦+128,₁₂=178𝑦⁵𝑥−858𝑦⁴𝑥+1688𝑦³𝑥−1776𝑥𝑦²+1024𝑥𝑦−256𝑥−38𝑦⁵𝑥²+262𝑦⁴𝑥²−640𝑦³𝑥²−512𝑥²𝑦−5𝑦⁶𝑥²+5𝑦⁶𝑥+800𝑥²𝑦²+128𝑥²,𝑎₂₁=108𝑦⁵−200𝑦⁴+200𝑦³−112𝑦²+32𝑦−28𝑦⁶+60𝑦⁶𝑥−344𝑦⁵𝑥+736𝑦⁴𝑥−784𝑦³𝑥+408𝑥𝑦²−32𝑥𝑦−64𝑥−24𝑦⁵𝑥²+16𝑦⁴𝑥²−96𝑥²𝑦+10𝑦⁶𝑥²+40𝑥²𝑦²+64𝑥²,𝑎₂₂=68𝑦⁵−140𝑦⁴+160𝑦³−96𝑦²+32𝑦−14𝑦⁶−50𝑦⁶𝑥²+270𝑦⁵𝑥²−640𝑦⁴𝑥²+820𝑦³𝑥²−560𝑥²𝑦²+160𝑥²𝑦+50𝑦⁶𝑥−270𝑦⁵𝑥+640𝑦⁴𝑥−820𝑦³𝑥+560𝑥𝑦²−160𝑥𝑦. (60)
From Definition 13, we know that 𝑅_2,2(𝑥,𝑦) is a BGIRI_NT since 𝑅_2,2(𝑥_𝑖,𝑦_𝑗)=𝑌_𝑖,𝑗 and deg𝑁(𝑥,𝑦)=8,deg𝐷(𝑥,𝑦)=6, and 𝐷(𝑥,𝑦)∣‖𝑁(𝑥,𝑦)‖².
In paper [9], a bivariate Thieletype matrix-valued rational interpolant 𝑅_𝑇(𝑥,𝑦)=𝐺_0,0(𝑦)+(𝑥−𝑥₀)/𝐺_1,0(𝑦)+⋯+(𝑥−𝑥_𝑚−1)/𝐺_𝑛,0(𝑦) (page 73), where for 𝑙=0,1,…,𝑛, 𝐺_𝑙,0𝑦=𝐵_𝑙,0𝑥₀,…,𝑥_𝑙,𝑦₀+𝑦−𝑦₀𝐵_𝑙,1𝑥₀,…,𝑥_𝑙,𝑦₀,𝑦₁+⋯+𝑦−𝑦_𝑚−1∣∣𝐵_𝑙,𝑚𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑚. (61) Here 𝐵_𝑙,𝑘(𝑥₀,…,𝑥_𝑙,𝑦₀,…,𝑦_𝑘),𝑙=0,1,…,𝑛,𝑘=0,1,…,𝑚 is defined different from that of this paper (see [9]). We now give another numerical example to compare the two algorithms, which shows that the method of this paper is better than the one of [9].

Table 1 

Example 18. Let 𝑌=𝑒𝑥,𝑦cos𝑥+𝑦sin𝑥+𝑦^𝑦𝑥𝑒^𝑦𝑥+𝑦 (62) and we suppose {𝑥₀,𝑥₁,𝑥₂,𝑥₃,𝑥₄,𝑥₅}={𝑦₀,𝑦₁,𝑦₂,𝑦₃,𝑦₄,𝑦₅}={1.0,1.2,1.4,1.6,1.8,2.0}. The numerical results 𝑅_NT(𝑥,𝑦) in Step 9 of Algorithm 7 and 𝑅_𝑇(𝑥,𝑦) in [9] are given in Table 2.

Table 2 

Remark 19. In Table 2, from the F-norm of 𝑅_NT−𝑌 and 𝑅_𝑇−𝑌, we can see that the error using Newton-Thiele type formula is much less than the one using Thieletype formula.

5. Conclusion

In this paper, the bivariate generalized inverse Newton-Thiele type matrix interpolation to approximate a matrix function 𝑓(𝑥,𝑦) is given, and we also give a recursive algorithm accompanied with some other important conclusions such as divisibility, characterization and some numerical examples. From the second example, it is easy to find that the approximant using Newton-Thiele type formula is much better than the one using Thieletype formula.

Acknowledgments

The work is supported by Shanghai Natural Science Foundation (10ZR1410900) and by Key Disciplines of Shanghai Municipality (S30104).