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

Authors: Rongrong Cui; Chuanqing Gu
Journal: Journal of Applied Mathematics
 (2013)

---

## 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).


---

## Body

## 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. LetCu×v consists of all u×v matrices with their elements in the complex plan C, and let A=(aij),B=(bij),andA,B∈Cu×v.Definition 1 (see [8]).
The scalar product of matricesA and B is defined by
(1)A·B=∑i=1u‍∑j=1v‍aijbij=tr(AB*),
where B* denotes the transpose of B and the Euclidean norm of A is given by
(2)∥A∥=(∑i=1u‍∑j=1v‍|aij|2)1/2.
It follows from Definition 1 and (2) that
(3)A·A-=∑i=1u‍∑j=1v‍|aij|2=∥A∥2=tr(AA-*),
where A- denotes the complex conjugate matrix of A and A-* denotes the complex conjugate transpose matrix of A. On the basis of (2) and (3), the generalized inverse of the matrix A is defined as
(4)Ar-1=1A=A-∥A∥2,A≠0,A∈Cu×v,
in particular, for A∈Ru×v,Ar-1=1/A=A/∥A∥2,A≠0.By means of generalized inverseAr-1 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Λm,n={(xl,yk):l=0,1,…,m,k=0,1,…,n,(xl,yk)∈ℝ2}.A matrix-valued set(5)Πm,n={Yl,k:Yl,k=Y(xl,yk)∈Cu×v,(xl,yk)∈Λm,n}.We need to find a bivariate matrix rational function(6)Rm,n(x,y)=N(x,y)D(x,y),
whose numerator N(x,y) is a complex or real polynomial matrix and denominator D(x,y) is a real polynomial, and
(7)Rm,n(xl,yk)=N(xl,yk)D(xl,yk)=Yl,k,(xl,yk)∈Λm,n.First, some notations need to be given as follows:(8)φ0,0(xl,yk)=Yl,k,∀(xl,yk)∈Λm,n,φ1,0(xi,xj,yk)=φ0,0(xj,yk)-φ0,0(xi,yk)xj-xi,φl,0(xp,…,xq,xi,xj,yk)=(φl-1,0(xp,…,xq,xj,yk)as-φl-1,0(xp,…,xq,xi,yk))×(xj-xi)-1,φl,1(xp,…,xq,xi,xj,yr,ys)=(ys-yr)×(φl,0(xp,…,xq,xi,xj,ys)(ys-ssssyr)×-φl,0(xp,…,xq,xi,xj,yr))-1,φl,k(xp,…,xq,xi,xj,ya,…,yb,yr,ys)=(ys-yr)×(φl,k-1(xp,…,xq,xi,xj,ya,…,yb,ys)-φl,k-1(xp,…,xq,xi,xj,ya,…,yb,yr))-1,
where the first subscript l of φ means the number of nodes xp,…,xq,xi,xj minus 1, and the second subscript k of φ means the number of nodes ya,…,yb,yr,ys minus 1.For simplicity, we let the nodesxp,…,xq,xi,xj be replaced by x0,…,xl, and the nodes ya,…,yb,yr,ys be replaced by y0,…,yk, 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:
(9)φ0,0(xl,yk)=Yl,k,φ1,0(x0,x1,yk)=φ0,0(x1,yk)-φ0,0(x0,yk)x1-x0,φl,0(x0,…,xl,yk)=(φl-1,0(x0,…,xl-2,xl,yk)-φl-1,0(x0,…,xl-2,xl-1,yk))×(xl-xl-1)-1,φl,1(x0,…,xl,y0,y1)=(y1-y0)φl,0(x0,…,xl,y1)-φl,0(x0,…,xl,y0),φ0,k(xl,y0,…,yk)=(yk-yk-1)×(φ0,k-1(xl,y0,…,yk-2,yk)-φ0,k-1(xl,y0,…,yk-2,yk-1))-1,φl,k(x0,x1,…,xl,y0,…,yk)=(yk-yk-1)×(φl,k-1(x0,…,xl,y0,…,yk-2,yk)-φl,k-1(x0,…,xl,y0,…,yk-2,yk-1))-1.
We assume that for all l,xl≠xl-1 and
(10)φl,k-1(x0,…,xl,y0,…,yk-2,yk)≠φl,k-1(x0,…,xl,y0,…,yk-2,yk-1),∀l,k.
From (9) matrix-valued Newton-Thiele type continued fractions for two-variable function can be constructed as follows:
(11)Rm,n(x,y)=U0(y)+U1(y)(x-x0)+⋯+Um(y)(x-x0)×(x-x1)⋯(x-xm-1),
where for l=0,1,…,m(12)Ul(y)=φl,0(x0,…,xl,y0)+y-y0φl,1(x0,…,xl,y0,y1)+⋯+y-yn-1∣∣φl,n(x0,…,xl,y0,…,yn).Remark 3.
Because of the construction,D(x,y) in Rm,n(x,y) mentioned as (6) is actually D(y) and each Ul(y) depends on n.We can also construct the antithetical form of bivariate matrix continued fraction in light of Definition4.Definition 4.
By means of generalized inverse (4), we define bivariate antithetical Newton-Thiele type matrix blending differences
(13)ψ0,0(xl,yk)=Yl,k,ψ0,1(xl,y0,y1)=(ψ0,0(xl,y1)-ψ0,0(xl,y0))(y1-y0)ψ0,k(xl,y0,…,yk)=(ψ0,k-1(xl,y0,…,yk-2,yk)-ψ0,k-1(xl,y0,…,yk-2,yk-1))×(yk-yk-1)-1,ψ1,k(x0,x1,y0,…,yk)=(x1-x0)ψ0,k(x1,y0,…,yk)-ψ0,k(x0,y0,…,yk),ψl,0(x0,…,xl,yk)=(xl-xl-1)×(ψl-1,0(x0,…,xl-2,xl,yk)-ψl-1,0(x0,…,xl-1,yk))-1,ψl,k(x0,x1,…,xl,y0,…,yk)=(xl-xl-1)×(ψl-1,k(x0,…,xl-2,xl,y0,…,yk)-ψl-1,k(x0,…,xl-1,y0,…,yk))-1.
From (13) we can get the antithetical formula for two-variable function
(14)R^m,n(x,y)=V0(x)+V1(x)(y-y0)+⋯+Vn(x)(y-y0)(y-y1)⋯(y-yn-1),
where for k=0,1,…,n(15)Vk(x)=ψ0,k(x0,y0,…,yk)+x-x0ψ1,k(x0,x1,y0,…,yk)+⋯+x-xm-1∣∣ψm,k(x0,…,xm,y0,…,yk).Now we will give two theorems aboutRm,n(x,y) as in (11) and (12) and R^m,n(x,y) as in (14) and (15). First of all, some notations need to be defined.In (12), for l=0,1,…,m, let
(16)Ul(s)(y)=φl,s(x0,x1,…,xl,y0,…,ys)+y-ysUl(s+1)(y),s=0,1,…,n-1,
where
(17)Ul(n)(y)=φl,n(x0,x1,…,xl,y0,…,yn),Ul(0)(y)=Ul(y).
Similarly, in (15), for k=0,1,…,n, let
(18)Vk(s)(x)=ψs,k(x0,x1,…,xs,y0,…,yk)+x-xsVk(s+1)(x),s=0,1,…,m-1,
where
(19)Vk(m)(x)=ψm,k(x0,x1,…,xm,y0,…,yk),Vk(0)(x)=Vk(x).Theorem 5.
Let(i)
φl,k(x0,…,xl,y0,…,yk),l=0,1,…,m;k=0,1,…,n exist and be nonzero (except for  φl,0(x0,…,xl,y0)),(ii)
Ul(s)(y)=φl,s(x0,x1,…,xl,y0,…,ys)+((y-ys)/Ul(s+1)(y)) satisfy Ul(s+1)(ys)≠0,  s=0,1,…,n-1, then Rm,n(x,y) as in (11) and (12) exists such that Rm,n(xl,yk)=Ylk,(xl,yk)∈Λm,n.Proof.
For simpleness, let(20)φl,k(x0,…,xl,y0,…,yk)=φl,k
if the conditions hold, thus for l=0,1,…,m;k=0,1,…,n(21)Ul(yk)=Ul(0)(yk)=φl,0+yk-y0φl,1+⋯+yk-yk-1∣∣φl,n+yk-yk∣∣Ul(k+1)(yk)(yk-yk)/Ul(k+1)(yk)=0, since Ul(k+1)(yk)≠0, then we get(22)Ul(yk)=φl,0(x0,…,xl,y0)+yk-y0φl,1(x0,…,xl,y0,y1)+⋯+yk-yk-2∣∣φl,k-1(x0,…,xl,y0,…,yk-1)+yk-yk-1∣∣φl,k(x0,…,xl,y0,…,yk)=φl,0(x0,…,xl,y0)+yk-y0φl,1(x0,…,xl,y0,y1)+⋯+yk-yk-2∣∣φl,k-1(x0,…,xl,y0,…,yk-1)+yk-yk-1∣∣(yk-yk-1)×(φl,k-1(x0,…,xl,y0,…,yk-2,yk)-φl,k-1(x0,…,xl,y0,…,yk-1))-1=φl,0(x0,…,xl,y0)+yk-y0φl,1(x0,…,xl,y0,y1)+⋯+yk-yk-2∣∣φl,k-1(x0,…,xl,y0,…,yk-2,yk)=⋯=φl,0(x0,…,xl,y0)+yk-y0φl,1(x0,…,xl,y0,yk)=φl,0(x0,…,xl,yk).We can use (9), (11), and (22) to find that
(23)Rm,n(xl,yk)=U0(yk)+U1(yk)(xl-x0)+⋯+Ul(yk)(xl-x0)(xl-x1)⋯(xl-xl-1)=Y(xl,yk).We can similarly prove the following result.Theorem 6.
Let(i)
ψl,k(x0,…,xl,y0,…,yk),l=0,1,…,m;k=0,1,…,n exist and be nonzero (except for  ψ0,k(x0,y0,…,yk)),(ii)
Vk(s)(x)=ψs,k(x0,x1,…,xs,y0,…,yk)+((x-xs)/Vk(s+1)(x)) satisfy Vk(s+1)(xs)≠0,  s=0,1,…,m-1, then R^m,n(x,y) as in (14) and (15) exists such that R^m,n(xl,yk)=Ylk,(xl,yk)∈Λm,n.Now we give the Newton-Thiele Matrix Interpolation Algorithm (NTMIA).Algorithm 7 (NTMIA).
Input:  {(xi,yj),Yi,j}(i=0,1,…,m;j=0,1,…,n)Output:  R^m(x,y).Step  1. For all (xi,yj)∈Λm,n, let φ(xi,yj)=Yi,j.Step  2. For j=0,1,…,n;p=1,2,…,m;i=p,p+1,…,m, compute
(24)φp,0(x0,…,xp-1,xi,yj)a=φp-1,0(x0,…,xp-2,xi,yj)-φp-1,0(x0,…,xp-1,yj)xi-xp-1.Step  3. For i=0,1,…,m;q=1,2,…,n;j=q,q+1,…,n, compute
(25)φi,q(x0,…,xi,y0,…,yq-1,yj)=(yj-yq-1)×(φi,q-1(x0,…,xi,y0,…,yq-2,yj)-φi,q-1(x0,…,xi,y0,…,yq-1))-1=((yj-yq-1)(φi,q-1(x0,…,xi,y0,…,yq-2,yj)--φi,q-1(x0,…,xi,y0,…,yq-1))-)×(∥φ(x0,…,xi,y0,…,yq-2,yj)-φ(x0,…,xi,y0,…,yq-1)∥2)-1.Step  4. For i=0,1,…,m;j=0,1,…,n, compute ai,j=φi,j(x0,…,xi,y0,…,yj).Step  5. For i=0,1,…,m;j=1,…,n judge if ai,j≠0; if yes go to Step 6, otherwise exist and show “the algorithm is not valid.”Step  6. For i=0,1,…,m;j=n-2,n-3,…,0(26)Qi,n=1,Pi,n=ai,n,Qi,n-1=∥ai,n∥2,Pi,n-1=ai,n-1Qi,n-1+(y-yn-1)Pi,n¯,Qi,j=∥ai,j+1∥2Qi,j+1+2(y-yj+1)tr(ai,j+1Pi,j+2¯*)+(y-yj+1)2Qi,j+2,Pi,j=ai,jQi,j+(y-yj)Pi,j+1¯.Step  7. For i=0,1,…,m;j=n,n-1,…,1, let Bi,j(y)=Pi,j(y)/Qi,j(y); if Bi,j(yj-1)≠0 then go to Step 8, otherwise exist and show “the algorithm is not valid.”Step  8. For i=0,1,…,m, let B^i(y)=Bi,0(y)=Pi,0(y)/Qi,0(y).Step  9. For k=1,…,m let R^0(x,y)=B^0(y), then compute R^k(x,y)=R^k-1(x,y)+⋯+(x-x0)⋯(x-xk-1)B^k(y).


## 3. Some Properties

Lemma 8 (see [5]).
DefineU0(y)=φ0,0(x0,y0)+(y-y0)/φ0,1(x0,y0,y1)+⋯+(y-yn-1)/φ0,n(x0,y0,…,yn); if we use generalized inverse by a tail-to-head rationalization, then a polynomial matrix E0(y) and a real polynomial F0(y) exist such that(i)
U0(y)=E0(y)/F0(y),F0(y)≥0,(ii)
F0(y)∣∥E0(y)∥2,“∣” means the sign of divisibility.Definition 9.
R
(
x
,
y
)
=
N
(
x
,
y
)
/
D
(
x
,
y
) is defined to be of type [l/t] if deg{Nij}≤l for 1≤i≤u,1≤j≤v,deg{Nij}=l for some (i,j) and deg{D}=t, where N(x,y)=(Nij(x,y))∈Cu×v.Lemma 10 (see [5]).
LetU0(y) be defined as in Lemma 8, and U0(y)=E0(y)/F0(y); if n is even, U0(y) is of type [n/n]; if n is odd, U0(y) is of type [n/n-1].Theorem 11.
Letφl,k∈Cu×v, (xl,yk)∈Λm,n, and x,y∈R. Define Ul(y)ofRmn(x,y) 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 N(x,y) and a real polynomial D(x,y) exist such that(i)
Rm,n(x,y)=N(x,y)/D(x,y),(ii)
D(x,y)≥0,(iii)
D(x,y)∣∥N(x,y)∥2.Proof.
Consider the following algorithm for the construction ofN(x,y), D(x,y), and Rmn(x,y).
Initialization: letD0(x,y)=1, and
(27)N0(x,y)=U0(y)=φ0,0+y-y0φ0,1+⋯+y-yn-1∣∣φ0,n.
By Lemma 8, a polynomial matrix E0(y) and a real polynomial F0(y) exist such that(a1)
N0(x,y)=E0(y)/F0(y),(a2)
F0(y)≥0,(a3)
F0(y)∣∥E0(y)∥2.
Recursion: forj=1,2,…,m let
(28)S(j)(x,y)=S(j-1)(x,y)+Uj(y)(x-x0)⋯(x-xj-1)
with Uj(y)=Ej(y)/Fj(y),Fj(y)∣∥Ej(y)∥2 having the representation(b1)
S(j)(x,y)=Nj(x,y)/Dj(x,y),(b2)
Dj(x,y)≥0,(b3)
Dj(x,y)∣∥Nj(x,y)∥2,where Dj(x,y) is a real polynomial.
Then we can get(29)S(j+1)(x,y)=S(j)(x,y)+Uj+1(x-x0)⋯(x-xj)=Nj(x,y)Dj(x,y)+Ej+1(y)Fj+1(y)(x-x0)⋯(x-xj)=(Nj(x,y)Fj+1(y)+Dj(x,y)Ej+1(y)×(x-x0)⋯(x-xj))×(Dj(x,y)Fj+1(y))-1=Nj+1(x,y)Dj+1(x,y),
where
(30)Nj+1(x,y)=Nj(x,y)Fj+1(y)+Dj(x,y)Ej+1(y)(x-x0)⋯(x-xj),Dj+1(x,y)=Dj(x,y)Fj+1(y).
Obviously,
(31)∥Nj+1∥2=∥Nj(x,y)∥2F(j+1)2(y)+Dj2(x,y)×∥E(j+1)(y)∥2(x-x0)2⋯(x-xj)2+Nj(x,y)·E-j+1(y)Fj+1(y)Dj(x,y)×(x-x0)⋯(x-xj)+N-j(x,y)·Ej+1(y)Fj+1(y)Dj(x,y)×(x-x0)⋯(x-xj)
since
(32)Dj(x,y)∣∥Nj(x,y)∥2,Fj+1(y)∣∥Ej+1(y)∥2,Dj+1(x,y)∣∥Nj+1(x,y)∥2.
Termination: S(m)(x,y)=Rm,n(x,y).Theorem 12 (Characterization).
LetRmn(x,y)=N(x,y)/D(x,y) be expressed as
(33)R=Rm,n(x,y)=U0(y)+U1(y)(x-x0)+⋯+Um(y)(x-x0)⋯(x-xm-1),
where Ul(y)=Ei/Fi, as in (12) we assume Fi,Fj have no common factor, i,j=0,1,…,m,i≠j, then(i)
if n is even,R is of type [mn+m+n/mn+n];(ii)
if n is odd,R is of type [mn+n/mn+n-m-1].Proof.
The proof is recursive. Letn be even.
Form=0,R0,n(x,y)=U0(y)=E0/F0, by Lemma 10, we find that it is of type [n/n].
Form=1, because deg{E0}=deg{F0}=n,deg{E1}=deg{F1}=n(34)R1,n(x,y)=N1D1=U0(y)+U1(y)(x-x0)=E0F0+E1F1(x-x0)=E0F1+E1F0(x-x0)F0F1
is of type [2n+1/2n].
Form=2, from the formula above we know that deg{N1}=2n+1,deg{D1}=2n,deg{E2}=n,anddeg{F2}=n. one has
(35)R2,n(x,y)=N2D2=N1D1+E2F2(x-x0)(x-x1)=N1F2+E2D1(x-x0)(x-x1)D1F2.
It is easy to find that R2,n(x,y) is of type [3n+2/3n].
Form=k, let Rk,n(x,y)=Nk/Dk be of type [(k+1)n+k/(k+1)n].
Form=k+1, we consider that
(36)Rk+1,n(x,y)=NkDk+Ek+1Fk+1(x-x0)⋯(x-xk)=NkFk+1+Ek+1Dk(x-x0)⋯(x-xk)DkFk=Nk+1Dk+1,
where
(37)deg{Nk+1}=(k+1)(n+1)+n=(k+2)n+(k+1),deg{Dk+1}=(k+1)n+n=(k+2)n.
Therefore, when n is even, R is of type [mn+m+n/mn+n].
Similarly, ifn is odd, for m=0,R0,n(x,y)=U0(y)=E0/F0, by Lemma 10, we find that it is of type [n/n-1].
Form=1, because deg{E0}=deg{E1}=n,deg{F0}=deg{F1}=n-1(38)R1,n(x,y)=N1D1=U0(y)+U1(y)(x-x0)=E0F0+E1F1(x-x0)=E0F1+E1F0(x-x0)F0F1
is of type [2n/2n-2].
Form=2, from the formula above we know that deg{N1}=2n,deg{D1}=2n-2,deg{E2}=n,anddeg{F2}=n-1. one has
(39)R2,n(x,y)=N2D2=N1D1+E2F2(x-x0)(x-x1)=N1F2+E2D1(x-x0)(x-x1)D1F2.
We get that R2,n(x,y) is of type [3n/3n-3].
Form=k, let
(40)Rk,n(x,y)=NkDk
be of type [(k+1)n/(k+1)(n-1)].
Form=k+1, we consider that
(41)Rk+1,n(x,y)=NkDk+Ek+1Fk+1(x-x0)⋯(x-xk)=NkFk+1+Ek+1Dk(x-x0)⋯(x-xk)DkFk=Nk+1Dk+1,
where
(42)deg{Nk+1}=(k+1)(n-1)+k+1=(k+2)n,deg{Dk+1}=(k+1)(n-1)+(n-1)=(k+2)(n-1).
Therefore, whenn is odd, R is of type [(m+1)n/(m+1)(n-1)].Definition 13.
A matrix-valued Newton-Thiele type rational fractionRm,n(x,y)=N(x,y)/D(x,y) is defined to be a bivariate generalized inverse and rational interpolant (BGIRINT) on the rectangular grid Λm,n if(i)
Rmn(xl,yk)=Ylk,(xl,yk)∈Λm,n,(ii)
D(xl,yk)≠0,(xl,yk)∈Λm,n,(iii)
(a) ifn is even,
(43)deg{N(x,y)}=mn+m+n,deg{D(x,y)}=mn+n,
(b) if n is odd,
(44)deg{N(x,y)}=mn+n,deg{D(x,y)}=mn+n-m-1,(iv)
D(x,y)∣∥N(x,y)∥2,(v)
D(x,y) is real, and D(x,y)≥0.
Now let us turn to the error estimate of theBGIRINT. Firstly, we give Lemma 14.Lemma 14.
Let(45)Rm,n*(x,y)=U0*(y)+U1*(y)(x-x0)+⋯+Um*(y)(x-x0)(x-x1)⋯(x-xm-1),
where for l=0,1,…,m(46)Ul*(y)=φl,0(x0,…,xl,y0)+y-y0φl,1(x0,…,xl,y0,y1)+⋯+y-yn-1∣∣φl,n(x0,…,xl,y0,…,yn)+y-yn∣∣φl,n+1(x0,…,xl,y0,…,yn,y),
then Rm,n*(x,y) satisfies
(47)Rm,n*(xi,y)=f(xi,y),i=0,1,…,m,Rm,n*(x,yj)=Rm,n(x,yj),j=0,1,…,n.Remark 15.
We delete the proof of the above lemma since it can be easily generalized from [3].Now the error estimate Theorem16 will be put which is also motivated by [3]. We let
(48)Rm,n(x,y)=N(x,y)D(x,y),Rm,n*(x,y)=N*(x,y)D*(x,y).Theorem 16.
Let a matrixf(x,y) be max(m+1,n+1) times continuously differentiable on a set D-={a≤x≤b,c≤y≤d} containing the points Λm,n={(xl,yk):l=0,1,…,m,k=0,1,…,n,(xl,yk)∈ℝ2}, then, for any point (x,y)∈D-, there exists a point (ξ,η)∈D- such that
(49)f(x,y)-Rm,n(x,y)=1D(x,y)D*(x,y)×{ωm+1(x)(m+1)!∂m+1∂xm+1E1(ξ,y)+ϖn+1(y)(n+1)!∂n+1∂yn+1E2(x,η)},
where
(50)ωm+1(x)=(x-x0)(x-x1)⋯(x-xm),ϖn+1(y)=(y-y0)(y-y1)⋯(y-yn),E1(x,y)=D(x,y)D*(x,y)[f(x,y)-Rm,n*(x,y)],E2(x,y)=D(x,y)D*(x,y)[Rm,n*(x,y)-Rm,n(x,y)].Proof.
Let(51)E1(x,y)=D(x,y)D*(x,y)[f(x,y)-Rm,n*(x,y)],E2(x,y)=D(x,y)D*(x,y)[Rm,n*(x,y)-Rm,n(x,y)],E(x,y)=E1(x,y)+E2(x,y).
From Lemma 14, we know
(52)E1(xi,y)=0,i=0,1,…,m,
which results in
(53)E1(x,y)=ωm+1(x)(m+1)!∂m+1∂xm+1E1(x,y)|x=ξ,
where ωm+1(x)=(x-x0)(x-x1)⋯(x-xm), and ξ is a number contained in the interval (a,b) which may depend on y. Similarly, from
(54)E2(x,yj)=0,j=0,1,…,n,
we can get that
(55)E2(x,y)=ϖn+1(y)(n+1)!∂n+1∂yn+1E2(x,y)|y=η,
where ϖn+1(y)=(y-y0)(y-y1)⋯(y-yn), and η is a number contained in the interval (c,d) which may depend on x.
Therefore,(56)E(x,y)=E1(x,y)+E2(x,y)=ωm+1(x)(m+1)!∂m+1∂xm+1E1(x,y)|x=ξ+ϖn+1(y)(n+1)!∂n+1∂yn+1E2(x,y)|y=η.
The proof is thus completed.


## 4. Numerical Examples

Example 17.
Let(xi,yj) and Yi,j,i,j=0,1,2 be given in Table 1:
(57)R2,2(x,y)=U0(y)+U1(y)x+U2(y)x(x-1).
Using Algorithm 7, we can get
(58)U0(y)=(1000)+y(0001)+y-1∣∣(00-2/5-1/5),U1(y)=(0-100)+y(01/2-1/20)+y-1∣∣(0110),U2(y)=(01120)+y(0-3/51/50)+y-1∣∣(0-15/14-5/14-5/7).
Based on generalized inverse (4), we obtain
(59)R2,2(x,y)=(a11a12a21a22)×(70y6-368y5+892y4-1240y3+1040y2-512y+128)-1=N(x,y)D(x,y),
where
(60)a11=70y6-368y5+892y4-1240y3+1040y2-512y+128,a12=178y5x-858y4x+1688y3x-1776xy2+1024xy-256x-38y5x2+262y4x2-640y3x2-512x2y-5y6x2+5y6x+800x2y2+128x2,a21=108y5-200y4+200y3-112y2+32y-28y6+60y6x-344y5x+736y4x-784y3x+408xy2-32xy-64x-24y5x2+16y4x2-96x2y+10y6x2+40x2y2+64x2,a22=68y5-140y4+160y3-96y2+32y-14y6-50y6x2+270y5x2-640y4x2+820y3x2-560x2y2+160x2y+50y6x-270y5x+640y4x-820y3x+560xy2-160xy.
From Definition13, we know that R2,2(x,y) is a BGIRINT since R2,2(xi,yj)=Yi,j and degN(x,y)=8,degD(x,y)=6, and D(x,y)∣∥N(x,y)∥2.
In paper [9], a bivariate Thieletype matrix-valued rational interpolant RT(x,y)=G0,0(y)+(x-x0)/G1,0(y)+⋯+(x-xm-1)/Gn,0(y) (page 73), where for l=0,1,…,n,
(61)Gl,0(y)=Bl,0(x0,…,xl,y0)+y-y0Bl,1(x0,…,xl,y0,y1)+⋯+y-ym-1∣∣Bl,m(x0,…,xl,y0,…,ym).
Here Bl,k(x0,…,xl,y0,…,yk),l=0,1,…,n,k=0,1,…,m 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
y
j
/
x
i
x
0
=
0
x
1
=
1
x
2
=
2
y
0
=
0
(
1
0
0
0
)
(
1
-
1
0
0
)
(
1
0
1
0
)
y
1
=
1
(
1
0
0
1
)
(
1
0
-
1
1
)
(
1
-
1
0
1
)
y
2
=
2
(
1
0
-
1
0
)
(
1
1
-
1
0
)
(
1
1
0
-
1
)Example 18.
Let(62)Y(x,y)=(cos(x+y)sin(x+y)eyxeyx+y)
and we suppose {x0,x1,x2,x3,x4,x5}={y0,y1,y2,y3,y4,y5}={1.0,1.2,1.4,1.6,1.8,2.0}. The numerical results RNT(x,y) in Step 9 of Algorithm 7 and RT(x,y) in [9] are given in Table 2.Table 2
Y
(
1.3,1.3
)
R
NT
(
1.3,1.3
)
∥
R
NT
-
Y
∥
F
(
-
0.856888
0.5155013
3.669296
1.30000
3.66929
2.600
)
(
-
0.856888
0.5155010
3.669297
1.30000
3.66929
2.600
)
1.09
e
-
006
R
T
(
1.3,1.3
)
∥
R
T
-
Y
∥
F
(
-
0.98999
0.14111
5.47394
1.29999
5.47394
2.999
)
2.613697
Y
(
1.5,1.5
)
R
NT
(
1.5,1.5
)
∥
R
NT
-
Y
∥
F
(
-
0.989999
0.1411200
4.481689
1.50000
4.481689
3.00
)
(
-
0.989992
0.1411203
4.481688
1.50000
4.481688
3.00
)
9.15
e
-
007
R
T
(
1.5,1.5
)
∥
R
T
-
Y
∥
F
(
-
0.99829372
-
0.05837
5.47394
1.49999
5.47394
3.1999
)
1.43144Remark 19.
In Table2, from the F-norm of RNT-Y and RT-Y, 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 functionf(x,y) 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.



---

*Source: 642818-2013-03-10.xml*