# Continued Fraction Interpolation of Preserving Horizontal Asymptote

Authors: Yushan Zhao; Kaiwen Wu; Jieqing Tan
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The classical Thiele-type continued fraction interpolation is an important method of rational interpolation. However, the rational interpolation based on the classical Thiele-type continued fractions cannot maintain the horizontal asymptote when the interpolated function is of a horizontal asymptote. By means of the relationship between the leading coefficients of the numerator and the denominator and the reciprocal differences of the continued fraction interpolation, a novel algorithm for the continued fraction interpolation is constructed in an effort to preserve the horizontal asymptote while approximating the given function with a horizontal asymptote. The uniqueness of the interpolation problem is proved, an error estimation is given, and numerical examples are provided to verify the effectiveness of the presented algorithm.


---

## Body

## 1. Introduction

The classical Thiele-type continued fraction interpolation [1–9] is an important method of rational interpolation. Suppose y=fx is the interpolated function and yj=fxj,j=0,1,…,n where x0,x1,…,xn are n+1 different interpolating nodes. The classical Thiele-type continued fraction has the form as follows:(1)Rnx=b0+x−x0b1+x−x1b2+⋯+x−xn−1bn,where(2)bj=φx0,x1,x2,…,xj,j=0,1,…,n,is the jth inverse difference of the function fx with respect to x0,x1,…,xj, which can be calculated recursively as follows:(3)φxj=fxj,j=0,1,…,n,φxp,xq=xq−xpφxq−φxp,φxu,…,xv,xr,xs=xs−xrφxu,…,xv,xs−φxu,…,xv,xr.It is not difficult to show thatRnx is a rational function whose numerator and denominator are polynomials of degrees not exceeding n+1/2 and n/2, respectively, where u denotes the largest integer not exceeding u, and Rnx satisfies(4)Rnxj=yj,j=0,1,2,…,n.Thenth reciprocal difference ρx0,x1,…,xn of the function fx with respect to x0,x1,…,xn is defined recursively as follows:(5)ρxj=fxj,j=0,1,…,n,ρxp,xq=xq−xpρxq−ρxp,ρxu,…,xv,xr,xs=xs−xrρxu,…,xv,xs−ρxu,…,xv,xr+ρxu,…,xv.The inverse differences can be calculated via the reciprocal differences as follows [10, 11]:(6)φxu,…,xv,xr,xs=ρxu,…,xv,xr,xs−ρxu,…,xv.Let(7)Rnx=b0+x−x0b1+x−x1b2+⋯+x−xn−1bn=PnxQnx.Denote byLPnx the leading coefficient of the polynomial Pnx, then when n is odd, the reciprocal differences and the leading coefficients of the numerator polynomial and denominator polynomial of continued fraction interpolation have the following identity relationship [12]:(8)LPnx=1,LQnx=ρx0,x1,…,xn,when n is even,(9)LPnx=ρx0,x1,…,xn,(10)LQnx=1.The classical Thiele-type continued fraction interpolation may not necessarily maintain the original horizontal asymptote of the interpolated function when the interpolated function has a horizontal asymptote. This paper presents an algorithm to construct the continued fraction interpolation preserving the horizontal asymptote that the interpolated function possesses. The uniqueness of solution of the numerical problem is proved, an error estimation is worked out, and numerical examples are provided to show the effectiveness of the new algorithm.


## 2. The Algorithm for Continued Fraction Interpolation of Preserving Horizontal Asymptote

The problem for continued fraction interpolation of preserving horizontal asymptote: Lety=fx be defined in I and x0,x1,…,x2m−1⊂I be 2m distinct interpolation nodes such that yj=fxj,j=0,1,…,2m−1. Suppose y=fx has a horizontal asymptote y=A, i.e., limx⟶∞fx=A, where A is a constant. Our purpose is to seek for a rational function of the following form:(11)R2mx=b0+x−x0b1+x−x1b2+⋯+x−x2m−1b2m=P2mxQ2mx.such that(12)R2mxj=yj,i=0,1,…,2m−1,limx⟶∞R2mx=A,where(13)bl=φx0,x1,…,xl,l=0,1,…,2m.Sincex2m is unknown, the formula of inverse differences cannot be used to calculate b2m=φx0,x1,…,x2m directly.It is not difficult to show(14)degP2mx≤2m+12=m,degQ2mx≤2m2=m.R2mx can be written as in the following form:(15)R2mx=b0+x−x0b1+x−x1b2+⋯+x−x2m−1b2m=P2mxQ2mx=cmxm+cm−1xm−1+⋯+c1x+c0dmxm+dm−1xm−1+⋯+d1x+d0.It follows from (9) and (10),(16)LP2mxLQ2mx=ρx0,x1,…,x2m1=ρx0,x1,…,x2m=cmdm=limx⟶∞R2mx.Using the relationship between the inverse differences and the reciprocal differences gives(17)ρx0,x1,…,x2m=b0+b2+⋯+b2m−2+b2m.With (12), (16), and (17) in mind, we have(18)limx⟶∞R2mx=A=ρx0,x1,…,x2m=b0+b2+⋯+b2m−2+b2m,i.e.,(19)b2m=A−b0−b2−⋯−b2m−4−b2m−2.As a result, the continued fraction interpolation with the preserved horizontal asymptote is given by(20)R2mx=b0+x−x0b1+x−x1b2+⋯+x−x2m−2b2m−1+x−x2m−1A−b0−b2−⋯−b2m−4−b2m−2.


## 3. The Uniqueness of Interpolant

Theorem 1.
Let(21)R2mx=b0+x−x0b1+x−x1b2+⋯+x−x2m−2b2m−1+x−x2m−1A−b0−b2−⋯−b2m−4−b2m−2.If all theith inverse differences bj=φx0,x1,x2,…,xj,j=0,1,2,…,2m−1 exist and A−b0−b2−⋯−b2m−4−b2m−2≠0, then set (see [13, 14])(22)t2mx=A−b0−b2−⋯−b2m−4−b2m−2,(23)tkx=bk+x−xktk+1x,k=2m−1,…,1,0.If(24)tk+1xk≠0,k=2m−1,…,1,0,then, we have(25)R2mxj=yj,j=0,1,…,2m−1,limx⟶∞R2mx=A.Proof.
∀i∈0,1,2,…,2m−1, using (20), (22)–(24) gives(26)R2mxj=b0+xj−x0b1+xj−x1b2+⋯+xj−xj−1bj=b0+xj−x0b1+xj−x1b2+⋯+xj−xj−2bj−1+=b0+xj−x0b1+xj−x1b2+⋯+xj−xj−2φx0,x1,…,xj−2,xj=⋯=b0+xj−x0φx0,xj=φxj=fxj.
The following can be obtained with the formulas (16), (17), and (19):(27)limx⟶∞R2mx=ρx0,x1,…,x2m=b0+b2+⋯+b2m−2+b2m=A.
The proof is completed.Theorem 2.
If a rational interpolation function of typem,m with the preserved horizontal asymptote exists, it must be the unique one.Proof.
Suppose the two rational functions of typem,m.(28)R2mx=P2mxQ2mx=cmxm+cm−1xm−1+⋯+c1x+c0dmxm+dm−1xm−1+⋯+d1x+d0,R2m∗x=P2m∗xQ2m∗x=cm∗xm+cm−1∗xm−1+⋯+c1∗x+c0∗dm∗xm+dm−1∗xm−1+⋯+d1∗x+d0∗.
Both meet the interpolation conditions in formula (12), we have(29)P2mxjQ2mxj=P2m∗xjQ2m∗xj,j=0,1,…,2m−1,limx⟶∞P2mxQ2mx=cmdm=limx⟶∞P2m∗xQ2m∗x=cm∗dm∗=A.
That is,(30)P2mxjQ2m∗xj−P2m∗xjQ2mxj=0,j=0,1,…,2m−1,cmdm∗−cm∗dm=0.
Since the leading term ofP2mxQ2m∗x−P2m∗xQ2mx is cmdm∗−cm∗dmx2m=0, P2mxQ2m∗x−P2m∗xQ2mx turns out to be a polynomial of degree not exceeding 2m−1, which has 2m distinct zeros. Therefore,(31)P2mxQ2m∗x≡P2m∗xQ2mx,namely,(32)R2mx≡R2m∗x.


## 4. The Error Estimation

Theorem 3.
Supposec,d is the smallest interval containing X2m−1=x0,x1,…,x2m−1 and fx is 2m times differentiable in c,d. Let(33)R2mx=b0+x−x0b1+x−x1b2+⋯+x−x2m−2b2m−1+x−x2m−1A−b0−b2−⋯−b2m−4−b2m−2=P2mxQ2mx,satisfy R2mxj=yj,j=0,1,…,2m−1,andlimx⟶∞R2mx=A. Then, for each x∈c,d, there exists a point ξ∈c,d such that(34)fx−R2mx=ω2mxQ2mx⋅fxQ2mxx=ξ2m2m!,and(35)limx⟶∞fx−R2mx=0,where ω2mx=x−x0x−x1⋯x−x2m−1.Proof.
LetEx=fxQ2mx−P2mx. Then, from R2mxj=yjj=0,1,…,2m−1, it follows,(36)Exj=0,j=0,1,…,2m−1.
Using the Lagrange interpolation formula with remainder term yields (see [15]),(37)Ex=ω2mx⋅fxQ2mx−P2mxx=ξ2m2m!=ω2mx⋅fxQ2mxx=ξ2m2m!.
Therefore,(38)fx−R2mx=ω2mxQ2mx⋅fxQ2mxx=ξ2m2m!,and(39)limx⟶∞fx−R2mx=limx⟶∞fx−limx⟶∞R2mx=A−A=0.


## 5. Numerical Examples

Example 1.
Given six interpolation nodesx0=1,x1=4,x2=8,x3=12,x4=16,x5=20. Suppose fx=arctanx, then, fx0=0.78539816,fx1=1.32581766,fx2=1.44644133,fx3=1.48765509,fx4=1.50837752,fx5=1.52083793 and limx⟶∞arctanx=1.57079633. We want to construct the continued fraction interpolant R6x such that it meets the interpolation conditions and limx⟶∞R6x=1.57079633 (keep eight decimal places).
According to what is known, the involved inverse differences can be calculated as shown in Table1.
By equation (19), we have(40)b6=A−b0−b2−b4=1.57079633−0.78539816−0.79395317−−0.00856673=0.00001173.
Substitutingb6 into R6x gives(41)R6x=0.78539816+x−15.55124306+x−40.79395317+x−8−1400.48475216+x−12−0.00856673+x−161707699.07647493+x−200.00001173,which can be simplified as(42)P6x=1.57079633x3−0.415046576808782x2+0.641166706097807x−0.112286748066322,Q6x=x3+0.372393353974841x2+0.645238679864088x+0.127305141976055,so,(43)R6x=1.57079633x3−0.415046576808782x2+0.641166706097807x−0.112286748066322x3+0.372393353974841x2+0.645238679864088x+0.127305141976055.
It is obvious that(44)limx⟶∞R6x=1.57079633.
Using the classical Thiele-type continued fraction interpolation, one can get(45)R5x=0.78539816+x−15.55124306+x−40.79395317+x−8−1400.48475216+x−12−0.00856673+x−161707699.07647493,which can be simplified as(46)R5x=0.5855832647x3+1569477.287x2−416900.8348x+225451.515999183.1503x2+370155.227x+385222.0773.
A comparison is made between the curvesy=fx and y=R6x as shown in Figure 1. The values of fx−R6x and fx−R5x at certain points are calculated as shown in Table 2. The errors fx−R6x and fx−R5x are illustrated in Figure 2.Table 1
Table of inverse differences.
jxjfxjφx0,xjφx0,x1,xjφx0,x1,xjφx0,x1,x3,xjφx0,x1,x3,x4,x5010.78539816141.325817665.55124306281.4464413310.589323560.793953173121.4876550915.663782740.79109702−1400.484752164161.5083775220.747480450.78966916−1867.40726226−0.008566735201.5208379325.834882520.78881307−2334.58517597−0.008564391707699.07647493Figure 1
Function plots offx and R6x.Table 2
Numerical experiments about the error.
xfx−R6xfx−R5x100.00000000520.0000000750500.00000000430.00000677341000.00000000190.000030074610000.00000000250.0005507282100000.00000000310.0058245238Figure 2
Function plots offx−R6x and fx−R5x.Example 2.
Given six interpolation nodesxj=j+2,j=0,1,…,5. Suppose fx=1+x33/x, then, fx0=1.04004191,fx1=1.01219632,fx2=1.00518144,fx3=1.00265959,fx4=1.00154083,fx5=1.00097087, and limx⟶∞1+x33/x=1. We want to construct the continued fraction interpolant R6x such that it meets the interpolation conditions and limx⟶∞R6x=1.
According to what is known, the involved inverse differences can be calculated as shown in Table3.
From equation (19), it follows:(47)b6=A−b0−b2−b4=1−1.04004191−−0.04659998−0.00713454=−0.00057647.
Substitutingb6 into R6x gives(48)R6x=1.04004191+x−2−35.91233252+x−3−0.04659998+x−4669.58185365+x−50.00713454+x−6−15861.70879642+x−7−0.00057647,which can be simplified as(49)P6x=x3−0.0270442257372862x2+0.0708647913780895x+0.587171073596883,Q6x=x3−0.0270376783697333x2+0.0706956613059999x+0.255582822560929.
So,(50)R6x=x3−0.0270442257372862x2+0.0708647913780895x+0.587171073596883x3−0.0270376783697333x2+0.0706956613059999x+0.255582822560929.
Clearly,(51)limx⟶∞R6x=1.
Using the classical Thiele-type continued fraction interpolation, one can get(52)R5x=1.04004191+x−2−35.91233252+x−3−0.04659998+x−4669.58185365+x−50.00713454+x−6−15861.70879642,which can be simplified as(53)R5x=−0.00006304490976x3+0.9617508703x2−1.53746219x+1.1433669680.96005036232x2−1.518908091x+1.03883064.
A comparison is made between the curvesy=fx and y=R6x, as shown in Figure 3. The values of fx−R6x and fx−R5x at certain points are calculated, as shown in Table 4. The errors fx−R6x and fx−R5x are shown in Figure 4.Table 3
Table of inverse differences.
jxjfxjφx0,xjφx0,x1,xjφx0,x1,xjφx0,x1,x3,xjφx0,x1,x3,x4,x5021.04004191131.01219632−35.91233252241.00518144−57.37156992−0.04659998351.00265959−80.25183125−0.04510651669.58185365461.00154083−103.89319748−0.04413007809.745032690.00713454571.00097087−127.97203165−0.04345007952.407421060.00707150−15861.70879642Figure 3
Function plots offx and R6x.Table 4
Numerical experiments about the error.
xfx−R6xfx−R5x100.000000001730.0000155716500.00000007590.00191754511000.00000005020.005057737010000.00000000640.0640174977100000.00000000010.6550176460Figure 4
Function plots offx−R6x and fx−R5x.Example 3.
Given six interpolation nodesxj=j/2+2,j=0,1,…,5. Suppose fx=1/2πe−x2/2+1, then, fx0=1.05399097,fx1=1.01752830,fx2=1.00443185,fx3=1.00087268,fx4=1.00013383,fx5=1.00001598 and limx⟶∞1/2πe−x2/2+1=1. We try to construct a continued fraction interpolant R6x such that it meets the interpolation conditions and limx⟶∞R6x=1.
According to what is known, the involved inverse differences can be calculated and listed in the following Table5. From (19), it follows:(54)b6=A−b0−b2−b4=1−1.05399097−−0.07733634−0.03057007=−0.00722470.
Substitutingb6 into R6x gives(55)R6x=1.05399097+x−2−13.71265611+x−2.5−0.07733634+x−358.85639112+x−3.50.03057007+x−4−259.12153694+x−4.5−0.00722470,which can be simplified as(56)P6x=x3−6.51105513479376x2+14.6566808471970x−10.9463472624099,Q6x=x3−6.51354079716349x2+14.6793603725717x−10.9982988663582.
So,(57)R6x=x3−6.51105513479376x2+14.6566808471970x−10.9463472624099x3−6.51354079716349x2+14.6793603725717x−10.9982988663582.
Obviously,(58)limx⟶∞R6x=1.
Using the classical Thiele-type continued fraction interpolation, we have(59)R5x=1.05399097+x−2−13.71265611+x−2.5−0.07733634+x−358.85639112+x−3.50.03057007+x−4−259.12153694,which can be simplified as(60)R5x=−0.003859192918x3+0.8758988689x2−2.912238657x+3.3502795540.8257816176x2−2.694383623x+3.03301528.
A comparison is conducted between the curvesy=fx and y=R6x as shown in Figure 5. The values of fx−R6x and fx−R5x at certain points are calculated and listed in Table 6. The errors fx−R6x and fx−R5x are shown in Figure 6.Table 5
Table of inverse differences.
jxjfxjφx0,xjφx0,x1,xjφx0,x1,xjφx0,x1,x3,xjφx0,x1,x3,x4,x5021.0539909712.51.01752830−13.71265611231.00443185−20.17792161−0.0773363433.51.00087268−28.23886414−0.0688410958.85639112441.00013383−37.13528304−0.0640406475.212256360.0305700754.51.00001598−46.31775448−0.0613401093.772010860.02864048−259.12153694Figure 5
Function plots offx and R6x.Table 6
Numerical experiments about the error.
xfx−R6xfx−R5x100.00015218070.0120808928500.00004689490.19023571231000.00002418040.422890115010000.00000247914.6280382082100000.000000248546.688387239Figure 6
Function plots offx−R6x and fx−R5x.


## 6. Conclusion

As classical approximation tool, continued fractions have been playing an important role in numerical rational approximation. However, continued fractions are rarely involved in shape-preserving design which is an interesting research topic in geometric modeling. In this paper, we construct an interpolating rational function based on the continued fractions, which serves to approximate the functions with the horizontal asymptotes. An algorithm is presented for the interpolating rational function to preserve the horizontal asymptote, the uniqueness of the interpolating rational function is proved and the error is analyzed. Numerical examples are given to verify the effectiveness of the new method.



---

*Source: 5662542-2022-06-01.xml*