# PIA and PPIA for Interpolating Points and Derivatives at Endpoints by Bézier Curves

Authors: Yeqing Yi; Lijuan Hu; Chengzhi Liu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this study, we are concerned with the interpolation problem that interpolates not only a given set of points but also derivatives at endpoints by using Bézier curves. The progressive iterative approximation (PIA) is proposed to interpolate the given points and derivatives. To speed up the convergence rate of PIA, we exploit the preconditioned PIA (PPIA), in which the diagonally compensated reduction is used to construct the preconditioner. The convergence of PIA and PPIA is also analyzed in this study. The proposed PPIA is applied to approximate higher order or/and rational Bézier curves. Numerical examples are given to illustrate the effectiveness of the proposed methods.


---

## Body

## 1. Introduction

We begin with the description of the data interpolation problem with constrains at endpoints. Given a set of organized pointspii=0n∈ℝ2 or ℝ3, rth derivatives d0rr=1,2,…,u∈ℝ2 or ℝ3 at the endpoint p0 and sth derivatives dnss=1,2,…,v∈ℝ2 or ℝ3 at the endpoint pn. For i=0,1,…,n, the ith point pi assigned an ordered parameter ti, i.e., t0<t1<⋯<tn. We want to find a m degree Bézier curve,(1)Ct=∑j=0mqjBjmt,that interpolates these points as well as the derivatives at the endpoints, i.e.,(2)Cti=pi,i=0,1,…,n;drCtdtr|t=t0=d0r,r=1,2,…,u;dsCtdts|t=tn=dns,s=1,2,…,v.The termsBjmt in (1) are the Bernstein polynomials of degree m, i.e., Bjmt=tj1−tm−j. According to the existence and uniqueness of polynomial interpolant, since there are n+u+v+1 different interpolation conditions in (2), we have m=n+u+v.Data interpolation and interpolation with tangents, curvatures, and other derivatives appear frequently in the fields of science and engineering [1, 2]. It is difficult to represent curves with complex shapes by a single curve, and a composite curve composed of several pieces of curves is necessary [3]. Hence, we need to control the smoothness of the composite curve when piecing curves together. It would be great if the derivatives at the endpoints are known. Very often, this knowledge is not given in practice, and we have to determine the derivatives at the endpoints. A careful choice of end conditions is important because it determines the shape of the interpolating curve near the endpoints. There are several strategies, e.g., natural end condition, Bessel end condition, quadratic end condition [4, 5]. In this study, we are especially interested in exploiting iterative methods for solving the interpolation curve (1) and do not consider the determination of constraints at the endpoints.It is known to all that the interpolation Bézier curve (1) can be obtained by solving a system of linear equations directly when the coefficient matrix is well-conditioned. However, the condition number of the linear system will be very large for large m; therefore, it is required to introduce an efficient algorithm to solve the linear system. In recent years, a geometric iterative method, named PIA, plays an important role in data interpolation. Due to its clear geometric meaning, stable convergence, and simple iterative scheme, PIA has intrigued the researchers for decades. For more details about PIA, readers can refer to a recent survey study [6], in which the authors summarized PIAs and their successful applications.In this study, we exploit the PIA format that interpolates not only a given set of points but also derivatives at endpoints by using Bézier curves. The study is organized as follows. In Section2, we propose the PIA format and then analyze its convergence. In order to accelerate the convergence rate of PIA, we exploit the preconditioning technique for PIA in Section 3. In Section 4, the proposed method is applied to approximate higher order or rational Bézier curves with constraints. Numerical examples are given to illustrate the effectiveness of our proposed methods in Section 5, and some conclusion remarks are given in the last section.


## 2. Progressive Iterative Approximation

### 2.1. Interpolating Derivatives Conditions at Endpoints

We note that thekth derivative of (1) is given by(3)dkCtdtk=m!m−k!∑j=0m−kΔkqjBjm−kt,where Δkqj is defined according to the recurrence relative formula:(4)Δkqj=Δk−1qj+1−Δk−1qj,withΔ0qj=qj.By direct calculation, from (2) and (3), we have(5)C0=q0=p0;C1=qm=pn;Δrq0=m−r!m!d0r,r=1,2,…,u;Δsqm−s=m−s!m!dns,s=1,2,…,v.Hence, the control pointsqjj=0,1,…,u,m−v,…,m in (1) can be determined according to (5). As an example, we outline the control points with respect to the case of u=v=3.(6)q0=p0,qm=pn,q1=q0+1md01,qm−1=qm−1mdn1,q2=m−2!m!d02+2q1−q0,qm−2=m−2!m!dn2+2qm−1−qm,q3=m−3!m!d03+3q2−3q1+q0,qm−3=−m−3!m!dn3+3qm−2−3qm−1+qm.Therefore, the Bézier curve (1) can be expressed by(7)Ct=∑j=0uqjBjmt+∑j=u+1m−v−1qjBjmt+∑j=m−vmqjBjmt,where the control points qjj=0,1,…,u,m−v,…,m can be obtained according to (5), and the control points qjj=u+1,u+2,…,m−v−1 need to be determined. We will give a computing method for solving these points in the following subsection.


### 2.2. Approximate Interpolation Algorithm

First, the given pointspii=1,2,…,n−1 and qjj=0,1,…,u,m−v,…,m are interpreted as the control points of a Bézier curve, and therefore, we can generate an initial approximate interpolation curve:(8)C0t=∑j=0uqjBjmt+∑j=u+1m−v−1qj0Bjmt+∑j=m−vmqjBjmt,where qj0=pj−uj=u+1,…,m−v−1.Second, forj=u+1,…,m−v−1, let δj0=pj−u−C0tj−u be the adjusting vector of the jth control point. Then, we can update the jth point qj1=qj0+δj0 and generate the first approximate interpolation curve:(9)C1t=∑j=0uqjBjmt+∑j=u+1m−v−1qj1Bjmt+∑j=m−vmqjBjmt.We repeat these procedures; thus, we obtain a sequence of approximate interpolation curves:(10)Ckt=∑j=0uqjBjmt+∑j=u+1m−v−1qjkBjmt+∑j=m−vmqjBjmt,k=1,2,…,where(11)qjk=qjk−1+δjk−1,δjk−1=pj−u−Ck−1tj−u,j=u+1,…,m−v−1.Thus, we obtain a sequence of approximate interpolation curvesCkt,k=0,1,…,. The initial Bézier curve is said to have the property of PIA with constraints at endpoints if(12)limk⟶∞Ckti=pi,i=0,1,…,n;limk⟶∞drCktdtr|t=t0=d0r,r=1,2,…,u;limk⟶∞dsCktdts|t=tn=dns,s=1,2,…,v.Remark 1.
By simple calculation, it is easy to verify that these approximate interpolation curvesCkt interpolate the endpoints p0andpn as well as the derivatives d0rr=1,2,…,u and dnss=1,2,…,v.
Let(13)Qk=qu+1k,qu+2k,…,qm−v−1kT,P¯=p1−∑j=0uqjBjmt1−∑j=m−vmqjBjmt1,…,pm−u−v−1−∑j=0uqjBjmtm−u−v−1−∑j=m−vmqjBjmtm−u−v−1T.
Then, the iterative process (11) can be written in the matrix form:(14)Qk=I−B¯Qk−1+P¯,k=1,2,…,where I is the n−1×n−1 identity matrix and B¯ is said to be the collocation matrix resulting from the Bernstein basis at parameters tii=1n−1, in detail(15)B¯=Bu+1mt1Bu+2mt1Bu+3mt1⋯Bm−v−1mt1Bu+1mt2Bu+2mt2Bu+3mt2⋯Bm−v−1mt2Bu+1mt3Bu+2mt3Bu+3mt3⋯Bm−v−1mt3⋮⋮⋮⋱⋮Bu+1mtn−1Br+2mtn−1Br+3mtn−1⋯Bm−v−1mtn−1.
The iterative process (14) is equivalent to the Richardson iteration for solving the system(16)B¯Q=P¯.


### 2.3. Convergence Analysis of PIA

Before analyzing the convergence, we review some definitions and conclusions.Definition 1 (see [7, 8]).
A matrix is called totally positive if all its minors are positive and positive stable (semistable) if all its eigenvalues have positive (nonnegative) real parts.
It is easy to verify that a positive stable matrix is nonsingular. We note that all the eigenvalues of a totally positive matrix are positive ([9]); thus, a totally positive matrix is also nonsingular.Definition 2 (see [7]).
A basisbiti=0n is called normalized if bit≥0i=0,…,n and ∑i=0nbit=1 and totally positive if its collocation matrix at any increasing sequence is a totally positive matrix.Lemma 1 (see [8]).
LetA be a given matrix. Then, A is positive stable if and only if for each nonzero vector x, there exists a positive definite matrix K, such that Rex∗KAx>0.Lemma 2 (see [10]).
LetA be an n×n matrix with entries aij. For i∈1,…,n, let Λi=∑j≠iaij be the sum of the absolute values of the nondiagonal entries in the ith row. Let Daii,Λi⊆ℂ be a closed disc centered at aii with radius Λi. Then, every eigenvalue of A lies within at least one of the discs Daii,Λi.Theorem 1.
The Bézier curve has the property of PIA with constraints at endpoints.Proof.
As stated in Remark1, the approximate interpolation curves Ckt interpolate endpoints p0andp1 and derivatives d0rr=1,2,…,u and dnss=1,2,…,v.
On the other hand, we insertu+v different parameters into the ordered parameters set tjj=0n; hence, we obtain a sequence of increasing parameters:(17)t0<t˜1<t˜2<⋯<t˜u<t1<t2<⋯<tn−1<t¯1<t¯2<⋯<t¯v<tn.
Note that the Bernstein basis is normalized and totally positive [7]. Let B be the collocation matrix resulting from the Bernstein basis at the parameter sequence (17). Then, it follows from Definition 2 that B is totally positive. Since the matrix B¯ defined in (15) is a submatrix of B, then B¯ is also totally positive according to Definition 1. Hence, all the eigenvalues of B¯ are positive.
Again, since the Bernstein basis is normalized, we haveB∞=1; thus, the infinite norm of the submatrix B¯ is less than 1, i.e., B¯∞<B∞=1a. Therefore, 0<λjB¯≤maxjλjB¯≤B¯∞<1,j=1,2,…,n−1. So, 0<λjI−B¯=1−λjB¯<1,j=1,2,…,n−1. This implies that the spectral radius of the iteration matrix I−B¯ is less than 1, i.e., ρI−B¯<1. Therefore, the iterative format (14) converges. The proof is thus complete.
Theorem1 tells us that the iteration process (14) always converges and we can obtain a satisfying approximate interpolant without solving a system of linear equations. However, as stated in [11], the collocation matrix resulting from the Bernstein basis is usually ill-conditioned, and the convergence rate of PIA is very slow. Therefore, it is necessary to introduce a preconditioner to speed up the convergence rate of PIA. In the following section, we propose a preconditioning technique for PIA by exploiting the properties of collocation matrix B¯.


## 2.1. Interpolating Derivatives Conditions at Endpoints

We note that thekth derivative of (1) is given by(3)dkCtdtk=m!m−k!∑j=0m−kΔkqjBjm−kt,where Δkqj is defined according to the recurrence relative formula:(4)Δkqj=Δk−1qj+1−Δk−1qj,withΔ0qj=qj.By direct calculation, from (2) and (3), we have(5)C0=q0=p0;C1=qm=pn;Δrq0=m−r!m!d0r,r=1,2,…,u;Δsqm−s=m−s!m!dns,s=1,2,…,v.Hence, the control pointsqjj=0,1,…,u,m−v,…,m in (1) can be determined according to (5). As an example, we outline the control points with respect to the case of u=v=3.(6)q0=p0,qm=pn,q1=q0+1md01,qm−1=qm−1mdn1,q2=m−2!m!d02+2q1−q0,qm−2=m−2!m!dn2+2qm−1−qm,q3=m−3!m!d03+3q2−3q1+q0,qm−3=−m−3!m!dn3+3qm−2−3qm−1+qm.Therefore, the Bézier curve (1) can be expressed by(7)Ct=∑j=0uqjBjmt+∑j=u+1m−v−1qjBjmt+∑j=m−vmqjBjmt,where the control points qjj=0,1,…,u,m−v,…,m can be obtained according to (5), and the control points qjj=u+1,u+2,…,m−v−1 need to be determined. We will give a computing method for solving these points in the following subsection.


## 2.2. Approximate Interpolation Algorithm

First, the given pointspii=1,2,…,n−1 and qjj=0,1,…,u,m−v,…,m are interpreted as the control points of a Bézier curve, and therefore, we can generate an initial approximate interpolation curve:(8)C0t=∑j=0uqjBjmt+∑j=u+1m−v−1qj0Bjmt+∑j=m−vmqjBjmt,where qj0=pj−uj=u+1,…,m−v−1.Second, forj=u+1,…,m−v−1, let δj0=pj−u−C0tj−u be the adjusting vector of the jth control point. Then, we can update the jth point qj1=qj0+δj0 and generate the first approximate interpolation curve:(9)C1t=∑j=0uqjBjmt+∑j=u+1m−v−1qj1Bjmt+∑j=m−vmqjBjmt.We repeat these procedures; thus, we obtain a sequence of approximate interpolation curves:(10)Ckt=∑j=0uqjBjmt+∑j=u+1m−v−1qjkBjmt+∑j=m−vmqjBjmt,k=1,2,…,where(11)qjk=qjk−1+δjk−1,δjk−1=pj−u−Ck−1tj−u,j=u+1,…,m−v−1.Thus, we obtain a sequence of approximate interpolation curvesCkt,k=0,1,…,. The initial Bézier curve is said to have the property of PIA with constraints at endpoints if(12)limk⟶∞Ckti=pi,i=0,1,…,n;limk⟶∞drCktdtr|t=t0=d0r,r=1,2,…,u;limk⟶∞dsCktdts|t=tn=dns,s=1,2,…,v.Remark 1.
By simple calculation, it is easy to verify that these approximate interpolation curvesCkt interpolate the endpoints p0andpn as well as the derivatives d0rr=1,2,…,u and dnss=1,2,…,v.
Let(13)Qk=qu+1k,qu+2k,…,qm−v−1kT,P¯=p1−∑j=0uqjBjmt1−∑j=m−vmqjBjmt1,…,pm−u−v−1−∑j=0uqjBjmtm−u−v−1−∑j=m−vmqjBjmtm−u−v−1T.
Then, the iterative process (11) can be written in the matrix form:(14)Qk=I−B¯Qk−1+P¯,k=1,2,…,where I is the n−1×n−1 identity matrix and B¯ is said to be the collocation matrix resulting from the Bernstein basis at parameters tii=1n−1, in detail(15)B¯=Bu+1mt1Bu+2mt1Bu+3mt1⋯Bm−v−1mt1Bu+1mt2Bu+2mt2Bu+3mt2⋯Bm−v−1mt2Bu+1mt3Bu+2mt3Bu+3mt3⋯Bm−v−1mt3⋮⋮⋮⋱⋮Bu+1mtn−1Br+2mtn−1Br+3mtn−1⋯Bm−v−1mtn−1.
The iterative process (14) is equivalent to the Richardson iteration for solving the system(16)B¯Q=P¯.


## 2.3. Convergence Analysis of PIA

Before analyzing the convergence, we review some definitions and conclusions.Definition 1 (see [7, 8]).
A matrix is called totally positive if all its minors are positive and positive stable (semistable) if all its eigenvalues have positive (nonnegative) real parts.
It is easy to verify that a positive stable matrix is nonsingular. We note that all the eigenvalues of a totally positive matrix are positive ([9]); thus, a totally positive matrix is also nonsingular.Definition 2 (see [7]).
A basisbiti=0n is called normalized if bit≥0i=0,…,n and ∑i=0nbit=1 and totally positive if its collocation matrix at any increasing sequence is a totally positive matrix.Lemma 1 (see [8]).
LetA be a given matrix. Then, A is positive stable if and only if for each nonzero vector x, there exists a positive definite matrix K, such that Rex∗KAx>0.Lemma 2 (see [10]).
LetA be an n×n matrix with entries aij. For i∈1,…,n, let Λi=∑j≠iaij be the sum of the absolute values of the nondiagonal entries in the ith row. Let Daii,Λi⊆ℂ be a closed disc centered at aii with radius Λi. Then, every eigenvalue of A lies within at least one of the discs Daii,Λi.Theorem 1.
The Bézier curve has the property of PIA with constraints at endpoints.Proof.
As stated in Remark1, the approximate interpolation curves Ckt interpolate endpoints p0andp1 and derivatives d0rr=1,2,…,u and dnss=1,2,…,v.
On the other hand, we insertu+v different parameters into the ordered parameters set tjj=0n; hence, we obtain a sequence of increasing parameters:(17)t0<t˜1<t˜2<⋯<t˜u<t1<t2<⋯<tn−1<t¯1<t¯2<⋯<t¯v<tn.
Note that the Bernstein basis is normalized and totally positive [7]. Let B be the collocation matrix resulting from the Bernstein basis at the parameter sequence (17). Then, it follows from Definition 2 that B is totally positive. Since the matrix B¯ defined in (15) is a submatrix of B, then B¯ is also totally positive according to Definition 1. Hence, all the eigenvalues of B¯ are positive.
Again, since the Bernstein basis is normalized, we haveB∞=1; thus, the infinite norm of the submatrix B¯ is less than 1, i.e., B¯∞<B∞=1a. Therefore, 0<λjB¯≤maxjλjB¯≤B¯∞<1,j=1,2,…,n−1. So, 0<λjI−B¯=1−λjB¯<1,j=1,2,…,n−1. This implies that the spectral radius of the iteration matrix I−B¯ is less than 1, i.e., ρI−B¯<1. Therefore, the iterative format (14) converges. The proof is thus complete.
Theorem1 tells us that the iteration process (14) always converges and we can obtain a satisfying approximate interpolant without solving a system of linear equations. However, as stated in [11], the collocation matrix resulting from the Bernstein basis is usually ill-conditioned, and the convergence rate of PIA is very slow. Therefore, it is necessary to introduce a preconditioner to speed up the convergence rate of PIA. In the following section, we propose a preconditioning technique for PIA by exploiting the properties of collocation matrix B¯.


## 3. Preconditioning Technique

### 3.1. Construction of Preconditioner

To derive a preconditioner for PIA, we first splitB¯ as B¯=B¯q+R, where(18)R=0⋯0Bu+q+1mt1⋯Bm−v−1mt1⋮⋱⋱⋱⋱⋮0⋱⋱⋱⋱Bm−v−1mtm−q−v−2Bu+1mtu+q+1⋱⋱⋱⋱0⋮⋱⋱⋱⋱⋮Bu+1mtn−1⋯Bm−q−2r−2mtn−10⋯0,and B¯q=B¯−R is a band matrix with bandwidth 2q+10≤q≤n−2.Then, we define a diagonal matrix(19)D=diagd1,d2,…,dn−1,such that De=Re with e=1,1,…,1T.Finally, let(20)Mq=B¯q+Dbe the diagonally compensated reduction matrix of B¯. Then, the matrix R is called the reduced matrix and D is called the compensation matrix for R [10, 12].IfMq is invertible, it can serve as a preconditioner for (16) and yield the preconditioned system(21)Mq−1B¯X=Mq−1P¯.


### 3.2. Preconditioned PIA

By applying Richardson method to the preconditioned system (21), we obtain(22)Qk=I−Mq−1B¯Qk−1+Mq−1P¯,=Qk−1+Mq−1P¯−B¯Qk−1,=Qk−1+Mq−1Δk−1,k=1,2,…,We refer the iterative process (22) as PPIA.


### 3.3. Convergence Analysis of PPIA

Theorem 2.
The preconditionerMq defined in (20) is invertible.Proof.
In the proof of Theorem1, we show that B¯ is positive stable. By Lemma 1, for any nontrivial y, there exists a positive definite matrix K, such that Rey∗KB¯y>0. According to the definitions of D and R, both the diagonal entries of D−R and the sum of the absolute values of the nondiagonal entries in the jth row of D−R are djj∈1,…,n−1. Let λ=a+b⋅i∈ℂ be any eigenvalue of D−R and x be the eigenvector corresponding to λ. Then, it follows from Lemma 2 that λ−dj=a−dj2+b2≤dj. Therefore, 0≤a≤2dj. This means that all the eigenvalues of D−R have nonnegative real parts and D−R is positive semistable. Hence, for the positive definite matrix K and the nontrivial x, we have Rex∗KD−Rx=x∗KxReλ≥0.
From (20), we have Mq=B¯+D−R, then Rex∗KMqx=Rex∗KB¯+D−Rx=Rex∗KB¯x+x∗KD−Rx=Rex∗KB¯x+Rex∗KD−Rx=Rex∗KB¯x+Rex∗Kλx>0.
It follows from Lemma1 that Mq is positive stable and there is no eigenvalue equal to zero. This completes the proof.
As it is stated in [13], the error matrix E¯=Mq−B¯ or the residual matrix R¯=I−Mq−1B¯ is often used to measure the efficiency of preconditioner Mq.
SinceMq=B¯+D−R, then(23)E¯∞=Mq−B¯∞=D−R∞=2R∞.
From (18) and (19), we have(24)D∞=R∞=max1≤j≤n−1∑u+1≤i≤m−v−1i−j>qBimtj.
As is stated in [12], Bj+q+1mtj>⋯>Bnmtj=0,Bj−q−1mtj>⋯>B0mtj=0. Thus, the sum of the entries of each column outside of the q diagonals decays to 0. This means that E¯∞ approaches to 0 as q increases. Combined with (23), we conclude that the matrix Mq is a good approximation to B¯ for large q. Thus, the spectrum of Mq−1B¯ is clustered around 1 for large q.Remark 2.
As it is known to all, the iterative process (22) is convergent if the spectral radius of the iteration matrix is less than 1, that is, ρI−Mq−1B¯<1. Note that the eigenvalues of I−Mq−1B¯ are 1−λiMq−1B¯. In particular, when the eigenvalues of Mq−1B¯ are distributed in 0,1, we have ρI−Mq−1B¯=1−λminMq−1B¯. Therefore, the convergence of (22) mainly depend on the value of λminMq−1B¯. Moreover, the larger the value of λminMq−1B¯ is, the faster the PPIA converges. It would be great if λMq−1B¯∈0,1. Unfortunately, we cannot give a rigid theoretical proof that λMq−1B¯∈0,1. Numerous experiments indicate that the eigenvalues of Mq−1B¯ are separated in 0,1, and most of them are clustered around 0. Based on the above analysis, we make a conjecture that the eigenvalues of Mq−1B¯ are separated in 0,1.
Recalling thatMq=B¯+D−R, we have B¯−1Mq=B¯−1B¯+D−R=I+B¯−1D−R; hence,(25)ρB¯−1Mq=ρI+B¯−1D−R≤1+ρB¯−1D−R≤1+B¯−1D−R∞.
Note that whenλiMq−1B¯>0, we have(26)λminMq−1B¯=1ρB¯−1Mq.
According to (25) and the conjecture in Remark 2, we have(27)ρI−Mq−1B¯=1−λminMq−1B¯=1−1ρB¯−1Mq≤1−11+B¯−1D−R∞<1.
Thus, PPIA converges.Remark 3.
In (27), we give a upper bound for ρI−Mq−1B¯. As mentioned earlier, D−R∞ tends to 0 as q increases; therefore, the upper bound 1−1/1+B¯−1D−R∞ also tends to 0. In particular, ρI−Mq−1B¯ equals to 0 when q=n−2. Hence, the spectral radius ρI−Mq−1B¯ may decrease as q increases, i.e., the bigger the value of q is, the faster the PPIA format converges.
On the other hand, it should be pointed out that at thekk=0,1,…,th iteration of PPIA, we have to calculate Mq−1Δk or solve the equation MqXk=Δk equivalently. It is costly, especially for large q. Therefore, we need to seek a tradeoff between the convergence rate and computational complexity. Experimentally, we found that q≈n/2 is a suitable choice.


## 3.1. Construction of Preconditioner

To derive a preconditioner for PIA, we first splitB¯ as B¯=B¯q+R, where(18)R=0⋯0Bu+q+1mt1⋯Bm−v−1mt1⋮⋱⋱⋱⋱⋮0⋱⋱⋱⋱Bm−v−1mtm−q−v−2Bu+1mtu+q+1⋱⋱⋱⋱0⋮⋱⋱⋱⋱⋮Bu+1mtn−1⋯Bm−q−2r−2mtn−10⋯0,and B¯q=B¯−R is a band matrix with bandwidth 2q+10≤q≤n−2.Then, we define a diagonal matrix(19)D=diagd1,d2,…,dn−1,such that De=Re with e=1,1,…,1T.Finally, let(20)Mq=B¯q+Dbe the diagonally compensated reduction matrix of B¯. Then, the matrix R is called the reduced matrix and D is called the compensation matrix for R [10, 12].IfMq is invertible, it can serve as a preconditioner for (16) and yield the preconditioned system(21)Mq−1B¯X=Mq−1P¯.


## 3.2. Preconditioned PIA

By applying Richardson method to the preconditioned system (21), we obtain(22)Qk=I−Mq−1B¯Qk−1+Mq−1P¯,=Qk−1+Mq−1P¯−B¯Qk−1,=Qk−1+Mq−1Δk−1,k=1,2,…,We refer the iterative process (22) as PPIA.


## 3.3. Convergence Analysis of PPIA

Theorem 2.
The preconditionerMq defined in (20) is invertible.Proof.
In the proof of Theorem1, we show that B¯ is positive stable. By Lemma 1, for any nontrivial y, there exists a positive definite matrix K, such that Rey∗KB¯y>0. According to the definitions of D and R, both the diagonal entries of D−R and the sum of the absolute values of the nondiagonal entries in the jth row of D−R are djj∈1,…,n−1. Let λ=a+b⋅i∈ℂ be any eigenvalue of D−R and x be the eigenvector corresponding to λ. Then, it follows from Lemma 2 that λ−dj=a−dj2+b2≤dj. Therefore, 0≤a≤2dj. This means that all the eigenvalues of D−R have nonnegative real parts and D−R is positive semistable. Hence, for the positive definite matrix K and the nontrivial x, we have Rex∗KD−Rx=x∗KxReλ≥0.
From (20), we have Mq=B¯+D−R, then Rex∗KMqx=Rex∗KB¯+D−Rx=Rex∗KB¯x+x∗KD−Rx=Rex∗KB¯x+Rex∗KD−Rx=Rex∗KB¯x+Rex∗Kλx>0.
It follows from Lemma1 that Mq is positive stable and there is no eigenvalue equal to zero. This completes the proof.
As it is stated in [13], the error matrix E¯=Mq−B¯ or the residual matrix R¯=I−Mq−1B¯ is often used to measure the efficiency of preconditioner Mq.
SinceMq=B¯+D−R, then(23)E¯∞=Mq−B¯∞=D−R∞=2R∞.
From (18) and (19), we have(24)D∞=R∞=max1≤j≤n−1∑u+1≤i≤m−v−1i−j>qBimtj.
As is stated in [12], Bj+q+1mtj>⋯>Bnmtj=0,Bj−q−1mtj>⋯>B0mtj=0. Thus, the sum of the entries of each column outside of the q diagonals decays to 0. This means that E¯∞ approaches to 0 as q increases. Combined with (23), we conclude that the matrix Mq is a good approximation to B¯ for large q. Thus, the spectrum of Mq−1B¯ is clustered around 1 for large q.Remark 2.
As it is known to all, the iterative process (22) is convergent if the spectral radius of the iteration matrix is less than 1, that is, ρI−Mq−1B¯<1. Note that the eigenvalues of I−Mq−1B¯ are 1−λiMq−1B¯. In particular, when the eigenvalues of Mq−1B¯ are distributed in 0,1, we have ρI−Mq−1B¯=1−λminMq−1B¯. Therefore, the convergence of (22) mainly depend on the value of λminMq−1B¯. Moreover, the larger the value of λminMq−1B¯ is, the faster the PPIA converges. It would be great if λMq−1B¯∈0,1. Unfortunately, we cannot give a rigid theoretical proof that λMq−1B¯∈0,1. Numerous experiments indicate that the eigenvalues of Mq−1B¯ are separated in 0,1, and most of them are clustered around 0. Based on the above analysis, we make a conjecture that the eigenvalues of Mq−1B¯ are separated in 0,1.
Recalling thatMq=B¯+D−R, we have B¯−1Mq=B¯−1B¯+D−R=I+B¯−1D−R; hence,(25)ρB¯−1Mq=ρI+B¯−1D−R≤1+ρB¯−1D−R≤1+B¯−1D−R∞.
Note that whenλiMq−1B¯>0, we have(26)λminMq−1B¯=1ρB¯−1Mq.
According to (25) and the conjecture in Remark 2, we have(27)ρI−Mq−1B¯=1−λminMq−1B¯=1−1ρB¯−1Mq≤1−11+B¯−1D−R∞<1.
Thus, PPIA converges.Remark 3.
In (27), we give a upper bound for ρI−Mq−1B¯. As mentioned earlier, D−R∞ tends to 0 as q increases; therefore, the upper bound 1−1/1+B¯−1D−R∞ also tends to 0. In particular, ρI−Mq−1B¯ equals to 0 when q=n−2. Hence, the spectral radius ρI−Mq−1B¯ may decrease as q increases, i.e., the bigger the value of q is, the faster the PPIA format converges.
On the other hand, it should be pointed out that at thekk=0,1,…,th iteration of PPIA, we have to calculate Mq−1Δk or solve the equation MqXk=Δk equivalently. It is costly, especially for large q. Therefore, we need to seek a tradeoff between the convergence rate and computational complexity. Experimentally, we found that q≈n/2 is a suitable choice.


## 4. Applications of PPIA with Constraints at Endpoints

Polynomial approximation for higher order or rational Bézier curves has been a continuous hotspot in computer-aided geometric design (CAGD). A number of approximation methods have been proposed in the literature, for example, polynomial approximation for rational Bézier curves [14–20] and polynomial approximation for higher Bézier curves (degree reduction) [3, 21–28].As a sample-based polynomial approximate method, Lu employed the weighted PIA (WPIA) to iteratively approximate the points sampled from the rational Bézier curves and achieved good approximations [17]. As mentioned in [19], Lu’s method is slow due to the slow convergence of WPIA and cannot deal with polynomial approximation with higher order continuity conditions at the endpoints. In fact, polynomial approximations with constraints are more practical than those without constraints. Due to the effectiveness in interpolation, these problems resulted from Lu’s method can be solved efficiently by employing the presented PPIA with constraints.


### 4.1. Polynomial Iterative Approximation for Higher Order and Rational Bézier Curves with Constraints

Given a rational Bézier curve of degreeN(28)Rt=∑j=0NωjvjBjNt∑j=0mωjBjNt,where vj and ωj∈ℝ+ are the control points and the associated weights, respectively.First, we samplen+1 points Rtjj=0n at the parameter values tjj=0,1,…,n and the derivatives Rrt0r=1u and Rst1s=1v at endpoints, where 0=t0<t1<⋯<tn=1. Then, PPIA for Bézier curves with constraints is employed to interpolate the points Rtjj=0n and the derivatives Rrt0r=1u and Rst1s=1v. Therefore, we can generate a sequence of m degree Bézier curves Ckt, which are the polynomial approximations of Rt.Remark 4.
We remark here that if all the weightsωj in (18) equal to 1, the rational Bézier curve (28) will degenerate into the Bézier curve of degree N. At this time, if m<N, it is the problem of approximating a given curve by a lower degree Bézier curve, which is an important approximation problem in CAGD, i.e., degree reduction of Bézier curves.
As stated in [17], the iterative method guarantees the convergence of data interpolation but does not necessary converge to the given curve Rt. Therefore, it is not necessary for PPIA to iterate infinitely. To ensure the computing efficiency, the iteration can be terminated if(29)Lpk+1≥θLpk,θ∈0,1,is satisfied, where Lpk is the approximation error at the kth iteration given by(30)Lpk=∫01Ckt−Rtpdt1/p,1≤p<∞.
Since there may exist rational functions in the integrand of (30), we need an effective and stable numerical method to calculate the integral (30). We note in [29] that the Gaussian rule or/and the adaptive quadrature method would be a better choice for rational cases. Therefore, we employ the τ-order Gauss–Legendre rule to calculate the integral (30).
By solving the roots of the Legendre polynomial of degreeτ, we can obtain Gauss points xi,i=1,2,…, and then calculate τ positive weights ωi. Finally, the integral in (30) can be evaluated approximately by(31)Lpk≈∑1≤i≤τwiCkxi−Rxip1/p.
We refer the reader to read [30] for more details about the τ-order Gauss–Legendre quadrature. In our tests, we used 15-order Gauss–Legendre quadrature to estimate the approximation error L2k, i.e., τ=15 and p=2.
Finally, we summarize the approximation algorithm in the following Algorithm1.Algorithm 1: Polynomial approximation for higher order and rational Bézier curves with constraints.
Input: original curveRt, largest admissible number kmax, and outer stopping tolerance θ in (19)Output: control points of the approximate curve, number of iterationsk, and approximation error Lpk+1(1)
Sample the pointsRtjj=0n and the derivatives Rrt0r=1u,Rst1s=1v(2)
Compute the initial errorL20 according to (31)(3)
Fork=0,1,…,kmax(i)
Employ PPIA for Bézier curves with constraints to generate thek+1th approximation curve(ii)
ComputeL2k+1 according to (31)(iii)
IfL2k+1≥θL2k, break forEnd for


## 4.1. Polynomial Iterative Approximation for Higher Order and Rational Bézier Curves with Constraints

Given a rational Bézier curve of degreeN(28)Rt=∑j=0NωjvjBjNt∑j=0mωjBjNt,where vj and ωj∈ℝ+ are the control points and the associated weights, respectively.First, we samplen+1 points Rtjj=0n at the parameter values tjj=0,1,…,n and the derivatives Rrt0r=1u and Rst1s=1v at endpoints, where 0=t0<t1<⋯<tn=1. Then, PPIA for Bézier curves with constraints is employed to interpolate the points Rtjj=0n and the derivatives Rrt0r=1u and Rst1s=1v. Therefore, we can generate a sequence of m degree Bézier curves Ckt, which are the polynomial approximations of Rt.Remark 4.
We remark here that if all the weightsωj in (18) equal to 1, the rational Bézier curve (28) will degenerate into the Bézier curve of degree N. At this time, if m<N, it is the problem of approximating a given curve by a lower degree Bézier curve, which is an important approximation problem in CAGD, i.e., degree reduction of Bézier curves.
As stated in [17], the iterative method guarantees the convergence of data interpolation but does not necessary converge to the given curve Rt. Therefore, it is not necessary for PPIA to iterate infinitely. To ensure the computing efficiency, the iteration can be terminated if(29)Lpk+1≥θLpk,θ∈0,1,is satisfied, where Lpk is the approximation error at the kth iteration given by(30)Lpk=∫01Ckt−Rtpdt1/p,1≤p<∞.
Since there may exist rational functions in the integrand of (30), we need an effective and stable numerical method to calculate the integral (30). We note in [29] that the Gaussian rule or/and the adaptive quadrature method would be a better choice for rational cases. Therefore, we employ the τ-order Gauss–Legendre rule to calculate the integral (30).
By solving the roots of the Legendre polynomial of degreeτ, we can obtain Gauss points xi,i=1,2,…, and then calculate τ positive weights ωi. Finally, the integral in (30) can be evaluated approximately by(31)Lpk≈∑1≤i≤τwiCkxi−Rxip1/p.
We refer the reader to read [30] for more details about the τ-order Gauss–Legendre quadrature. In our tests, we used 15-order Gauss–Legendre quadrature to estimate the approximation error L2k, i.e., τ=15 and p=2.
Finally, we summarize the approximation algorithm in the following Algorithm1.Algorithm 1: Polynomial approximation for higher order and rational Bézier curves with constraints.
Input: original curveRt, largest admissible number kmax, and outer stopping tolerance θ in (19)Output: control points of the approximate curve, number of iterationsk, and approximation error Lpk+1(1)
Sample the pointsRtjj=0n and the derivatives Rrt0r=1u,Rst1s=1v(2)
Compute the initial errorL20 according to (31)(3)
Fork=0,1,…,kmax(i)
Employ PPIA for Bézier curves with constraints to generate thek+1th approximation curve(ii)
ComputeL2k+1 according to (31)(iii)
IfL2k+1≥θL2k, break forEnd for


## 5. Numerical Examples

In this section, several numerical experiments are given to illustrate the effectiveness of the proposed methods. All the numerical experiments were carried out on a computer with Intel(R) Core(TM) i5-5200U CPU @2.20 GHz by Matlab R2012b.


### 5.1. Tests of PIA and PPIA with Constraints

Since we have shown in Theorem1 that the approximate interpolation curves Ckt interpolate the endpoints and the derivatives at the endpoints exactly, we need to measure the interpolation error at the points pi,i=1,2,…,n−1. Therefore, the interpolation error of the kth approximate interpolation curve Ckt can be expressed by(32)εk=max1≤i≤n−1Ckti−pi,where the norm is the Euclidean norm.Example 1.
We sample 9 pointspii=0,1,…,8, rr=1,2th derivatives d0r, and ss=1,2th derivatives dns from the semicircle of radius 5 in the following way:(33)pi=5cosπ8i,5sinπ8i,i=0,1,…,8;d01=0,5,d81=0,−5;d02=−5,0,d82=5,0.Example 2.
We sample 19 pointspii=0,1,…,18, rr=1,2,3th derivatives d0r, and ss=1,2,3th derivatives dns from the helix of radius 5 in the following way:(34)pi=5cosπ3i,5sinπ3i,π3i,i=0,1,…,18;d01=d181=0,5,1;d02=d182=−5,0,0.
PIA and PPIA formats for constraints at the endpoints are employed to test Example1 and Example 2. When we test PPIA, we take the half bandwidth q=4 in Example 1 and q=10 in Example 2. In Table 1, we list the spectral radii of the iteration matrices of PIA and PPIA with different constraints at the endpoints. For one thing, the spectral radii of the iteration matrices of PIA and PPIA are less than 1; this means that both PIA and PPIA formats converge. For another, the spectral radii of PPIA are much less that those of PIA; hence, we can expect that PPIA would have a better convergence behavior than PIA.
We list in Table2 the interpolation errors of PIA and PPIA when we test Example 1, and we list in Table 3 the interpolation errors of PIA and PPIA when we test Example 2. As we can see from Tables 2 and 3, the interpolation errors of PIA and WPIA decrease as the number of iterations increases, which indicates that both PIA and PPIA converge. With the same iterations, the interpolate errors of PPIA are much less than those of PIA. This means that the preconditioning technique can accelerate the convergence rate of PIA significantly.
In Figures1 and 2, we display the interpolation Bézier curves (dashed lines) when approximating the points with different derivatives given in Example 1. We also show in Figures 1 and 2 the semicircle (solid lines) for comparison. It should be pointed out that the approximate interpolation Bézier curves generated by PPIA almost coincide with the semicircle, no matter which endpoint conditions are used. In Figures 3 and 4, we display the interpolation Bézier curves when approximating the points combined with different derivatives given in Example 2. All the approximate interpolation Bézier curves in Figures 1–4 are obtained by preforming 5 PIA or PPIA iterations. Clearly, the IPPIA can yield better approximations than PIA.Table 1
Spectral radii of iteration matrices of PIA and PPIA in Example1 and Example 2.
d0u and dnvExample1Example2PIAPPIAPIAPPIAu=v=10.999182980.238468970.999999850.14929065u=1,v=20.997854110.059944280.999999850.32779082u=2,v=10.997854110.059944280.999999850.32779082u=v=20.999647440.578873710.999999850.24570266Table 2
Interpolation errors of PIA and PPIA in Example1 with different constraints at endpoints.
ku=v=1u=1,v=2u=2,v=1u=v=2PIAPPIAPIAPPIAPIAPPIAPIAPPIA15.4395e − 013.9638e − 031.4104e004.3811e − 021.4104e004.3811e − 029.1860e − 011.8233e − 0223.8274e − 011.4605e − 049.1358e − 012.9530e − 039.1358e − 012.9530e − 037.7540e − 017.6179e − 0432.8078e − 015.3842e − 067.7453e − 011.7897e − 047.7453e − 011.7897e − 046.6697e − 013.1829e − 0542.1528e − 011.9849e − 076.6933e − 011.0740e − 056.6933e − 011.0740e − 055.8133e − 011.3299e − 0651.7055e − 017.3177e − 095.8848e − 016.4387e − 075.8848e − 016.4387e − 075.1179e − 015.5567e − 08101.1297e − 015.6300e − 143.6339e − 015.2694e − 133.6339e − 014.9105e − 133.4893e − e − 018.0428e − 15158.6368e − 025.5417e − e − 143.4564e − 018.9706e − 143.4564e − 012.4995e − 143.4688e − 018.0428e − 15206.4608e − 025.5417e − 143.2913e − 018.9706e − 143.2913e − 012.4995e − 143.0648e − 018.0428e − 15402.4501e − 025.5417e − 142.4740e − 018.9706e − 142.4740e − 012.4995e − 141.6225e − 018.0428e − 15502.2440e − 025.5417e − 142.1753e − 018.9706e − 142.1753e − 012.4995e − 141.2254e − 018.0428e − 15Table 3
Interpolation errors of PIA and PPIA in Example2 with different constraints at endpoints.
ku=v=1u=1,v=2u=2,v=1u=v=2PIAPPIAPIAPPIAPIAPPIAPIAPPIA14.6885e007.5971e − 051.1145e015.9422e − 031.1191e015.9448e − 039.4468e002.1520e − 0324.2701e002.6533e − 067.0081e002.3159e − 036.9937e002.3164e − 035.7317e001.9501e − 0433.7030e001.5379e − 075.8458e007.6479e − 045.8565e007.6496e − 044.9157e003.7394e − 0543.1942e001.9400e − 084.6120e002.5081e − 044.6235e002.5087e − 043.5536e009.0175e − 0652.7764e002.8613e − 093.5412e008.2217e − 053.5411e008.2235e − 052.4739e002.2132e − 06101.5062e002.1517e − 132.4607e003.1113e − 072.4604e003.1120e − 071.9172e001.9817e − 09206.7818e − 011.7061e − 141.6841e003.3866e − 121.6823e004.6786e − 121.3289e001.9304e − 13502.5157e − 011.3557e − 141.2789e009.5874e − 131.2786e008.1082e − 139.9004e − 012.2176e − 131009.7764e − 022.4470e − 147.6721e − 019.6635e − 137.6715e − 019.3884e − 136.2363e − 011.6799e − 13Figure 1
Bézier curves interpolatingd01 and d81 in Example 1. (a) PIA. (b) PPIA.
(a)(b)Figure 2
Bézier curves interpolatingd02 and d82 in Example 1. (a) PIA. (b) PPIA.
(a)(b)Figure 3
Bézier curves obtained by performing 5 PIA or PPIA iterations in Example2. (a) Interpolating d01 and d181. (b) Interpolating d01 and d182.
(a)(b)Figure 4
Bézier curves obtained by performing 5 PIA or PPIA iterations in Example2. (a) Interpolating d02 and d181. (b) Interpolating d02 and d182.
(a)(b)


### 5.2. Tests of Polynomial Approximation for Higher Order and Rational Bézier Curves

First, we employ Algorithm1 to test the degree reduction of higher order Bézier curves in Example 3 and then test the polynomial approximation for rational Bézier curve in Example 4 and Example 5. In our tests, we take the endpoint interpolation conditions u=v=1 and take the half bandwidth q=⌈n/2⌉ when we construct the preconditioner for PIA, where ⌈⋅⌉ is the round-up operation. When we test Example 3 and Example 4, we take θ=0.9. Again, when we test Example 5, we take θ=0.98.Example 3 (Degree reduction [17]).
A Bézier curve of degree 15 is defined by the control points0,0,1.5,−2,4.5,−1,9,0,4.5,1.5,2.5,3,0,5,−4,8.5,3,9.5,4.4,10.5,6,12,8,11,9,10,9.5,5,7,6,and5,7.
By employing Algorithm1, we use the Bézier curves of degree mm=6,7,…,14 to approximate the Bézier curve given in Example 3. In Table 4, we list the iteration number k, the approximation errors L2k, and the elapsed CPU time T (in seconds) of Algorithm 1; the numerical results run by Lu’s method ([17]) are also listed for comparison. The results for degree reduction with constraints, reported in Table 4, indicate that Algorithm 1 outperforms Lu’s method, in terms of the approximation error and the computing time.
In Figure5(a), we show the approximations of degree 8 obtained by Algorithm 1 (blue dashed line) and Lu’s method (red dashdot line). In Figure 5(b), we display the corresponding error distance curves of the approximations shown in Figure 5(a). It is evident from Figure 5 that Algorithm 1 can produce better approximations than Lu’s method. What is more, the approximations obtained by Algorithm 1 are C1-continuity at two endpoints, whereas the approximations obtained by Lu’s method only satisfy the simplest C0-continuity.Table 4
Numerical results of Algorithm1 and Lu’s method in Example 3.
mAlgorithm1Lu’s methodkL2kT(s)kL2kT(s)613.6163e − 010.5934.9761e − 012.14711.3234e − 010.6043.6118e − 012.59819.2304e − 020.6071.9136e − 013.58912.1803e − 020.6091.0012e − 013.991018.6645e − 030.6097.4032e − 024.141115.8744e − 030.6196.3833e − 024.071222.5835e − 030.94114.7558e − 024.671312.3246e − 040.61114.5005e − 024.791431.1873e − 041.23152.6818e − 026.09Figure 5
Reduction from degree 15 to 8 in Example3. (a) Comparison of polynomial approximations. (b) The corresponding error distance curves.
(a)(b)Example 4.
(Polynomial approximation for rational Bézier curves [19]). A rational Bézier curve of degree 8 is defined by the control points 14,1, 34,25, 40,38, −12,24, 5,21, 26,7, 18,41, −13,34, and 14,1 and the associated weights 1,2,1/3,2,2,1/3,2,and1.Example 5.
(Polynomial approximation for rational Bézier curves [19]). A rational Bézier curve of degree 9 is defined by the control points 17,12, 32,34, −23,24, 33,6, −23,15, 25,3, 30,−2, −5,−8, −5,15, 11,8, and the associated weights 1,2,3,6,4,5,3,4,2,1.
We employ Algorithm1 and Lu’s method to generate polynomial approximations for the rational Bézier curves given in Example 4 and Example 5. In Table 5, we list the iteration number k, the approximation error L2k, and the elapsed CPU time T (in seconds) of Algorithm 1 and Lu’s method. The results indicate that Algorithm 1 also performs better than Lu’s method for polynomial approximation for rational Bézier curves with constraints.
Themm=9,12,15 degree approximations for the rational Bézier curve in Example 4 and their corresponding error distance curves are shown in Figures 6 and 7, respectively. Again, mm=10,12,15 degree approximations for the rational Bézier curve in Example 5 and their corresponding error distance curves are shown in Figures 8 and 9, respectively. It is already clear from these figures the reliability of the proposed approximation algorithm.Table 5
Numerical results of Algorithm1 and Lu’s method in Examples 4 and 5.
mExample4mExample5Algorithm1Lu’s methodAlgorithm1Lu’s methodkL2kT(s)kL2kT(s)kL2kT(s)kL2kT(s)811.9356e + 000.5851.9556e + 001.76915.7514e − 010.5989.9194e − 012.64911.0215e + 000.5951.6452e + 001.751022.9092e − 010.87107.6682e − 013.301016.6233e − 010.5851.4297e + 001.801111.5379e − 010.59125.8658e − 013.961114.0403e − 010.6051.2715e + 001.801226.9603e − 020.92164.2321e − 015.231212.6760e − 010.5951.1512e + 001.781314.3342e − 020.63183.1705e − 015.751311.6397e − 010.5951.0575e + 001.791432.0350e − 021.19182.4794e − 015.861411.2175e − 010.5877.7411e − 012.391521.5264e − 020.89201.8593e − 016.47Figure 6
Polynomial approximations of the rational Bézier curve in Example4 with different m. (a) m = 9. (b) m = 12. (c) m = 14.
(a)(b)(c)Figure 7
The corresponding error distance curves of the approximations in Figure6. (a) m = 9. (b) m = 12. (c) m = 14.
(a)(b)(c)Figure 8
Polynomial approximations of the rational Bézier curve in Example5 with different m. (a) m = 10. (b) m = 12. (c) m = 15.
(a)(b)(c)Figure 9
The corresponding error distance curves of the approximations in Figure8. (a) m = 10. (b) m = 12. (c) m = 15.
(a)(b)(c)


## 5.1. Tests of PIA and PPIA with Constraints

Since we have shown in Theorem1 that the approximate interpolation curves Ckt interpolate the endpoints and the derivatives at the endpoints exactly, we need to measure the interpolation error at the points pi,i=1,2,…,n−1. Therefore, the interpolation error of the kth approximate interpolation curve Ckt can be expressed by(32)εk=max1≤i≤n−1Ckti−pi,where the norm is the Euclidean norm.Example 1.
We sample 9 pointspii=0,1,…,8, rr=1,2th derivatives d0r, and ss=1,2th derivatives dns from the semicircle of radius 5 in the following way:(33)pi=5cosπ8i,5sinπ8i,i=0,1,…,8;d01=0,5,d81=0,−5;d02=−5,0,d82=5,0.Example 2.
We sample 19 pointspii=0,1,…,18, rr=1,2,3th derivatives d0r, and ss=1,2,3th derivatives dns from the helix of radius 5 in the following way:(34)pi=5cosπ3i,5sinπ3i,π3i,i=0,1,…,18;d01=d181=0,5,1;d02=d182=−5,0,0.
PIA and PPIA formats for constraints at the endpoints are employed to test Example1 and Example 2. When we test PPIA, we take the half bandwidth q=4 in Example 1 and q=10 in Example 2. In Table 1, we list the spectral radii of the iteration matrices of PIA and PPIA with different constraints at the endpoints. For one thing, the spectral radii of the iteration matrices of PIA and PPIA are less than 1; this means that both PIA and PPIA formats converge. For another, the spectral radii of PPIA are much less that those of PIA; hence, we can expect that PPIA would have a better convergence behavior than PIA.
We list in Table2 the interpolation errors of PIA and PPIA when we test Example 1, and we list in Table 3 the interpolation errors of PIA and PPIA when we test Example 2. As we can see from Tables 2 and 3, the interpolation errors of PIA and WPIA decrease as the number of iterations increases, which indicates that both PIA and PPIA converge. With the same iterations, the interpolate errors of PPIA are much less than those of PIA. This means that the preconditioning technique can accelerate the convergence rate of PIA significantly.
In Figures1 and 2, we display the interpolation Bézier curves (dashed lines) when approximating the points with different derivatives given in Example 1. We also show in Figures 1 and 2 the semicircle (solid lines) for comparison. It should be pointed out that the approximate interpolation Bézier curves generated by PPIA almost coincide with the semicircle, no matter which endpoint conditions are used. In Figures 3 and 4, we display the interpolation Bézier curves when approximating the points combined with different derivatives given in Example 2. All the approximate interpolation Bézier curves in Figures 1–4 are obtained by preforming 5 PIA or PPIA iterations. Clearly, the IPPIA can yield better approximations than PIA.Table 1
Spectral radii of iteration matrices of PIA and PPIA in Example1 and Example 2.
d0u and dnvExample1Example2PIAPPIAPIAPPIAu=v=10.999182980.238468970.999999850.14929065u=1,v=20.997854110.059944280.999999850.32779082u=2,v=10.997854110.059944280.999999850.32779082u=v=20.999647440.578873710.999999850.24570266Table 2
Interpolation errors of PIA and PPIA in Example1 with different constraints at endpoints.
ku=v=1u=1,v=2u=2,v=1u=v=2PIAPPIAPIAPPIAPIAPPIAPIAPPIA15.4395e − 013.9638e − 031.4104e004.3811e − 021.4104e004.3811e − 029.1860e − 011.8233e − 0223.8274e − 011.4605e − 049.1358e − 012.9530e − 039.1358e − 012.9530e − 037.7540e − 017.6179e − 0432.8078e − 015.3842e − 067.7453e − 011.7897e − 047.7453e − 011.7897e − 046.6697e − 013.1829e − 0542.1528e − 011.9849e − 076.6933e − 011.0740e − 056.6933e − 011.0740e − 055.8133e − 011.3299e − 0651.7055e − 017.3177e − 095.8848e − 016.4387e − 075.8848e − 016.4387e − 075.1179e − 015.5567e − 08101.1297e − 015.6300e − 143.6339e − 015.2694e − 133.6339e − 014.9105e − 133.4893e − e − 018.0428e − 15158.6368e − 025.5417e − e − 143.4564e − 018.9706e − 143.4564e − 012.4995e − 143.4688e − 018.0428e − 15206.4608e − 025.5417e − 143.2913e − 018.9706e − 143.2913e − 012.4995e − 143.0648e − 018.0428e − 15402.4501e − 025.5417e − 142.4740e − 018.9706e − 142.4740e − 012.4995e − 141.6225e − 018.0428e − 15502.2440e − 025.5417e − 142.1753e − 018.9706e − 142.1753e − 012.4995e − 141.2254e − 018.0428e − 15Table 3
Interpolation errors of PIA and PPIA in Example2 with different constraints at endpoints.
ku=v=1u=1,v=2u=2,v=1u=v=2PIAPPIAPIAPPIAPIAPPIAPIAPPIA14.6885e007.5971e − 051.1145e015.9422e − 031.1191e015.9448e − 039.4468e002.1520e − 0324.2701e002.6533e − 067.0081e002.3159e − 036.9937e002.3164e − 035.7317e001.9501e − 0433.7030e001.5379e − 075.8458e007.6479e − 045.8565e007.6496e − 044.9157e003.7394e − 0543.1942e001.9400e − 084.6120e002.5081e − 044.6235e002.5087e − 043.5536e009.0175e − 0652.7764e002.8613e − 093.5412e008.2217e − 053.5411e008.2235e − 052.4739e002.2132e − 06101.5062e002.1517e − 132.4607e003.1113e − 072.4604e003.1120e − 071.9172e001.9817e − 09206.7818e − 011.7061e − 141.6841e003.3866e − 121.6823e004.6786e − 121.3289e001.9304e − 13502.5157e − 011.3557e − 141.2789e009.5874e − 131.2786e008.1082e − 139.9004e − 012.2176e − 131009.7764e − 022.4470e − 147.6721e − 019.6635e − 137.6715e − 019.3884e − 136.2363e − 011.6799e − 13Figure 1
Bézier curves interpolatingd01 and d81 in Example 1. (a) PIA. (b) PPIA.
(a)(b)Figure 2
Bézier curves interpolatingd02 and d82 in Example 1. (a) PIA. (b) PPIA.
(a)(b)Figure 3
Bézier curves obtained by performing 5 PIA or PPIA iterations in Example2. (a) Interpolating d01 and d181. (b) Interpolating d01 and d182.
(a)(b)Figure 4
Bézier curves obtained by performing 5 PIA or PPIA iterations in Example2. (a) Interpolating d02 and d181. (b) Interpolating d02 and d182.
(a)(b)


## 5.2. Tests of Polynomial Approximation for Higher Order and Rational Bézier Curves

First, we employ Algorithm1 to test the degree reduction of higher order Bézier curves in Example 3 and then test the polynomial approximation for rational Bézier curve in Example 4 and Example 5. In our tests, we take the endpoint interpolation conditions u=v=1 and take the half bandwidth q=⌈n/2⌉ when we construct the preconditioner for PIA, where ⌈⋅⌉ is the round-up operation. When we test Example 3 and Example 4, we take θ=0.9. Again, when we test Example 5, we take θ=0.98.Example 3 (Degree reduction [17]).
A Bézier curve of degree 15 is defined by the control points0,0,1.5,−2,4.5,−1,9,0,4.5,1.5,2.5,3,0,5,−4,8.5,3,9.5,4.4,10.5,6,12,8,11,9,10,9.5,5,7,6,and5,7.
By employing Algorithm1, we use the Bézier curves of degree mm=6,7,…,14 to approximate the Bézier curve given in Example 3. In Table 4, we list the iteration number k, the approximation errors L2k, and the elapsed CPU time T (in seconds) of Algorithm 1; the numerical results run by Lu’s method ([17]) are also listed for comparison. The results for degree reduction with constraints, reported in Table 4, indicate that Algorithm 1 outperforms Lu’s method, in terms of the approximation error and the computing time.
In Figure5(a), we show the approximations of degree 8 obtained by Algorithm 1 (blue dashed line) and Lu’s method (red dashdot line). In Figure 5(b), we display the corresponding error distance curves of the approximations shown in Figure 5(a). It is evident from Figure 5 that Algorithm 1 can produce better approximations than Lu’s method. What is more, the approximations obtained by Algorithm 1 are C1-continuity at two endpoints, whereas the approximations obtained by Lu’s method only satisfy the simplest C0-continuity.Table 4
Numerical results of Algorithm1 and Lu’s method in Example 3.
mAlgorithm1Lu’s methodkL2kT(s)kL2kT(s)613.6163e − 010.5934.9761e − 012.14711.3234e − 010.6043.6118e − 012.59819.2304e − 020.6071.9136e − 013.58912.1803e − 020.6091.0012e − 013.991018.6645e − 030.6097.4032e − 024.141115.8744e − 030.6196.3833e − 024.071222.5835e − 030.94114.7558e − 024.671312.3246e − 040.61114.5005e − 024.791431.1873e − 041.23152.6818e − 026.09Figure 5
Reduction from degree 15 to 8 in Example3. (a) Comparison of polynomial approximations. (b) The corresponding error distance curves.
(a)(b)Example 4.
(Polynomial approximation for rational Bézier curves [19]). A rational Bézier curve of degree 8 is defined by the control points 14,1, 34,25, 40,38, −12,24, 5,21, 26,7, 18,41, −13,34, and 14,1 and the associated weights 1,2,1/3,2,2,1/3,2,and1.Example 5.
(Polynomial approximation for rational Bézier curves [19]). A rational Bézier curve of degree 9 is defined by the control points 17,12, 32,34, −23,24, 33,6, −23,15, 25,3, 30,−2, −5,−8, −5,15, 11,8, and the associated weights 1,2,3,6,4,5,3,4,2,1.
We employ Algorithm1 and Lu’s method to generate polynomial approximations for the rational Bézier curves given in Example 4 and Example 5. In Table 5, we list the iteration number k, the approximation error L2k, and the elapsed CPU time T (in seconds) of Algorithm 1 and Lu’s method. The results indicate that Algorithm 1 also performs better than Lu’s method for polynomial approximation for rational Bézier curves with constraints.
Themm=9,12,15 degree approximations for the rational Bézier curve in Example 4 and their corresponding error distance curves are shown in Figures 6 and 7, respectively. Again, mm=10,12,15 degree approximations for the rational Bézier curve in Example 5 and their corresponding error distance curves are shown in Figures 8 and 9, respectively. It is already clear from these figures the reliability of the proposed approximation algorithm.Table 5
Numerical results of Algorithm1 and Lu’s method in Examples 4 and 5.
mExample4mExample5Algorithm1Lu’s methodAlgorithm1Lu’s methodkL2kT(s)kL2kT(s)kL2kT(s)kL2kT(s)811.9356e + 000.5851.9556e + 001.76915.7514e − 010.5989.9194e − 012.64911.0215e + 000.5951.6452e + 001.751022.9092e − 010.87107.6682e − 013.301016.6233e − 010.5851.4297e + 001.801111.5379e − 010.59125.8658e − 013.961114.0403e − 010.6051.2715e + 001.801226.9603e − 020.92164.2321e − 015.231212.6760e − 010.5951.1512e + 001.781314.3342e − 020.63183.1705e − 015.751311.6397e − 010.5951.0575e + 001.791432.0350e − 021.19182.4794e − 015.861411.2175e − 010.5877.7411e − 012.391521.5264e − 020.89201.8593e − 016.47Figure 6
Polynomial approximations of the rational Bézier curve in Example4 with different m. (a) m = 9. (b) m = 12. (c) m = 14.
(a)(b)(c)Figure 7
The corresponding error distance curves of the approximations in Figure6. (a) m = 9. (b) m = 12. (c) m = 14.
(a)(b)(c)Figure 8
Polynomial approximations of the rational Bézier curve in Example5 with different m. (a) m = 10. (b) m = 12. (c) m = 15.
(a)(b)(c)Figure 9
The corresponding error distance curves of the approximations in Figure8. (a) m = 10. (b) m = 12. (c) m = 15.
(a)(b)(c)


## 6. Conclusion

The presented PIA format for Bézier curves iteratively interpolates not only a set of points but also derivatives at endpoints. It can construct higher order continuous piecewise interpolation Bézier curves without solving a linear system directly. To accelerate the convergence rate of PIA, the preconditioning technique is introduced. Numerical examples demonstrate that the given methods are effective, and the convergence rate of PIA is accelerated significantly by the presented preconditioning technique.Furthermore, due to the efficient performance in data interpolation, the presented PPIA for Bézier curves with constraints is applied to approximate higher order and rational Bézier curves. Our numerical experiments show that the proposed method achieves higher efficiency and results in a better approximation than the similar sample-based polynomial approximation method. More importantly, they can satisfyCkk=1,2,…, continuity at the endpoints.



---

*Source: 9940610-2021-07-01.xml*