A Note on Rational Lagrange Polynomials for CAGD Applications and Isogeometric Analysis

1. Introduction

Since the mid-1970s, research in computer-aided geometric design (CAGD) has heavily relied on Bernstein-Bézier polynomials as a foundational tool for representing and manipulating curves and surfaces. This preference stems primarily from their strong geometric intuition: the shape of a Bézier curve or surface is governed by a finite set of control points, which usually lie outside the geometric entity itself and act as handles for interactive design and editing. A second important reason is the mathematical flexibility of the Bernstein basis, which supports essential operations such as degree elevation, subdivision, and knot insertion, enabling local refinement without altering the global geometry’s key features for adaptive modeling and analysis [1].

Furthermore, Bernstein polynomials are central to the construction of non-uniform rational B-splines (NURBS), which extend polynomial representations by incorporating weights. This allows for the exact representation of conic sections and other analytic shapes, a capability critical to engineering design and manufacturing applications. As a result, Bernstein-based representations remain a cornerstone of modern CAD/CAM systems and form the basis of isogeometric analysis, which unifies geometric modeling ([1,2,3]) and finite element analysis ([4]).

It is well established in the literature that the weighted, or rational, Bernstein-Bézier formulation enables the exact representation of conic sections and quadrics, provided that suitable weights are assigned to the appropriate control points (e.g., see [1,2,3], among hundreds of other sources). This capability is one of the key advantages of rational Bézier and NURBS representations in computer-aided geometric design and CAD systems. However, there have also been a limited number of investigations into using Lagrange polynomials for similar geometric modeling tasks, particularly for approximating circular arcs and ellipses, but for efficient analysis as well [5,6]. Focusing on circular arcs, these approaches can produce visually acceptable results but typically suffer from small but non-negligible numerical errors due to the inability of polynomial Lagrange interpolation to exactly represent circular shapes [7,8,9,10].

On the other hand, Lagrange interpolation in itself is quite robust and efficient. The mentioned instabilities are related to Runge oscillations, which only occur on equidistant nodal distributions (for example, see Ref. [11]). Within this context, Lagrange interpolation might break down in such a case for high polynomial orders. However, a simple remedy is to use non-equidistant nodal distributions that accumulate nodes near the interval boundaries. This has been done in the Spectral Element Method (SEM), where Gauss-Lobatto-Legendre (GLL) or Gauss-Lobatto-Chebyshev (GLC) nodes are employed to replace equidistant nodes [12,13]. Such non-equidistant nodes result in excellent interpolation properties and mitigate any instabilities. Employing GLL or GLC nodes, it has been noticed that the element matrices are well conditioned and not numerical problems arise [14].

From the above discussion, it follows that Lagrange polynomials offer high accuracy when applied in conjunction with non-equidistant nodal distributions (e.g., GLL or GLC points), whereas Bernstein polynomials possess the advantage of incorporating weights that enable the exact representation of conics and quadrics.

In this paper, we lay the groundwork for combining these two advantages into a unified framework. Specifically, we determine the weights that can be directly applied to both uniform and non-uniform Lagrange polynomials, thereby constructing an equivalent functional basis of rational Lagrange polynomials.

2. General Formulation

2.1. Degree Elevation Using Bernstein-Bézier Polynomials

Rational Bézier curves can be readily handled within the framework of computer-aided geometric design (CAGD). Beginning with a low degree, such as $n=2$, and employing an appropriate set of homogeneous (projected) control points ${\mathbf{P}}_0^w,\dots ,{\mathbf{P}}_n^w$, with weights ($w_i,i=0,\dots ,n$), the shape of the curve remains invariant under degree elevation, with the control points being updated accordingly at each step as follows:

$${\mathbf{Q}}_i^w=(1-a_i){\mathbf{P}}_i^w+a_i{\mathbf{P}}_{i-1}^w,$$

(1)

with

$$\begin{array}{c}\hfill a_i={\displaystyle \frac{i}{n+1}},i=1,\dots ,n.\end{array}$$

(2)

The updated homogeneous weights $w_i^{*}$ that accompany the elevated control points ${\mathbf{Q}}_i^w=w_i^{*}{[x_i^{*},y_i^{*},z_i^{*}]}^T$ (or the 2-D analogue) follow exactly the same convex-combination pattern as for the degree elevation formula for Bernstein (polynomial) Bézier curves (cf. Eq. (1)):

$$w_i^{*}={\displaystyle \frac{i}{n+1}}{w}_{i-1}+\left(1-{\displaystyle \frac{i}{n+1}}\right)w_i,i=1,\dots ,n.$$

(3)

with the boundary values preserved,

$$w_0^{*}=w_0,w_{n+1}^{*}=w_n.$$

(4)

Once the homogeneous points are updated, we obtain the Cartesian control points simply by dehomogenising:

$${\mathbf{Q}}_i={\displaystyle \frac{{\mathbf{Q}}_i^w}{w_i^{*}}},i=0,\dots ,n+1.$$

(5)

This degree-elevation step may be applied repeatedly; after k successive elevations (to degree $n+k$) the same recurrence is used at each stage, guaranteeing that the rational Bézier curve’s shape remains unchanged while its control-polygon flexibility increases.

2.2. Replacing Bernstein Polynomials by Lagrange Polynomials

While the procedure of degree elevation is straightforward in conjunction with Bernstein polynomials, it is not the case with Lagrange polynomials. Of course, since both sets span the same space (of powers ${\xi }^i,i=0,\dots$), these two sets are linearly interrelated.

In more detail, let us consider the reference unit length $0\le \xi \le 1$. For a certain degree n, i.e., for the nodal points $0={\xi }_0,\dots ,{\xi }_n=1$, let $\mathbf{L}$ be the row vector including all relevant univariate Lagrange polynomials of degree n:

$$\mathbf{L}=[L_0\left(\xi \right),\dots ,L_n\left(\xi \right)],$$

(6)

and $\mathbf{B}$ be the row vector including all univariate Bernstein polynomials of degree n:

$$\mathbf{B}=[B_0\left(\xi \right),\dots ,B_n\left(\xi \right)].$$

(7)

Then, we can write:

$$\begin{array}{cc}\hfill \mathbf{L}& ={\mathbf{T}}_{BL}\mathbf{B},\hfill \\ \hfill \mathrm{and}\end{array}$$

(8)

$$\begin{array}{cc}\hfill \mathbf{B}& ={\mathbf{T}}_{LB}\mathbf{L},\hfill \end{array}$$

(9)

where ${\mathbf{T}}_{BL}$ (from Bernstein to Lagrange) and ${\mathbf{T}}_{LB}$ (from Lagrange to Bernstein) are transformation matrices, which have been discussed in detail elsewhere [15]. For the purposes of this paper, the determination of transformation matrix ${\mathbf{T}}_{LB}$ that must be applied to Lagrange polynomials–according to Eq. (9)–plays a critical role.

Actually, since a linear transformation occurs between two functional sets (here the subscripts denote, B: Bernstein-Bézier, L: Lagrange), not only the different sets of polynomials (cf. Eq. (8) and Eq. (9)), but also the homogeneous coordinates (${\mathbf{Q}}_L^w,{\mathbf{P}}_B^w$) and the weights (${\mathbf{w}}_L,{\mathbf{w}}_B$) are interrelated according to the following theorem.

Theorem

The homogeneous (projected) cordinates of nodal points (${\mathbf{Q}}_L^w$) associated with the true boundary of a curve $C\left(\xi \right)$ in the Lagrange formulation, are related with the projected control points (${\mathbf{P}}_B^w$) of the same curve in the Bernstein formulation, as follows:

$$\begin{array}{|c|}\hline {\mathbf{Q}}_L^w={\mathbf{T}}_{LB}^T{\mathbf{P}}_B^w\\ \hline\end{array}.$$

(10)

Moreover, the weights ${\mathbf{w}}_L$ in the Lagrange formulation are related with the weights ${\mathbf{w}}_B$ in the initial Bernstein formulation, as follows:

$$\begin{array}{|c|}\hline {\mathbf{w}}_L={\mathbf{T}}_{LB}^T{\mathbf{w}}_B\\ \hline\end{array}.$$

(11)

In more detail, to obtain the nodal points ${\mathbf{Q}}_L$ on the true curve, it is necessary to divide the vector of homogeneous coordinates ${\mathbf{Q}}_L^w$ (calculated using Eq. (10)) by the vector of associated weights ${\mathbf{w}}_L$ (estimated by Eq. (11)).

Proof. 

The proof is based on the fact that the shape $\mathbf{x}\left(\xi \right)$ of the curve under consideration is preseved (in shape and parametrically [2]). Therefore, expressing $\mathbf{x}\left(\xi \right)$ first in terms of Lagrange polynomials (forming the vector $\mathbf{L}$) and second in terms of Bernstein polynomials (forming the vector $\mathbf{B}$), we have:

$$\begin{array}{cc}\hfill \mathbf{x}\left(\xi \right)& ={\displaystyle \frac{{\mathbf{L}}^T{\mathbf{Q}}_L^w}{{\mathbf{L}}^T{\mathbf{w}}_L}}={\displaystyle \frac{{\mathbf{B}}^T{\mathbf{P}}_B^w}{{\mathbf{B}}^T{\mathbf{w}}_B}}.\hfill \end{array}$$

(12)

Substituting $\mathbf{B}$ through Eq. (9) into the utmost right term of Eq. (12), we receive:

$$\begin{array}{cc}\hfill \mathbf{x}\left(\xi \right)& ={\displaystyle \frac{{\mathbf{L}}^T{\mathbf{Q}}_L^w}{{\mathbf{L}}^T{\mathbf{w}}_L}}={\displaystyle \frac{{\left({\mathbf{T}}_{LB}\mathbf{L}\right)}^T{\mathbf{P}}_B^w}{{\left({\mathbf{T}}_{LB}\mathbf{L}\right)}^T{\mathbf{w}}_B}}={\displaystyle \frac{{\mathbf{L}}^T\left({\mathbf{T}}_{LB}^T{\mathbf{P}}_B^w\right)}{{\mathbf{L}}^T\left({\mathbf{T}}_{LB}^T{\mathbf{w}}_B\right)}}.\hfill \end{array}$$

(13)

Obviously, the last equality of Eq. (13) suggests that Eq. (10) and Eq. (11) provide the projected nodal points and associated weights, respectively, and this completes the proof.    □

Remark-1: Weights. In a previous work (see, [15], pp. 80–83])) regarding a circular arc of 90 degrees modeled by quadratic Bernstein polynomials, Eq. (11) has been proven from scratch (after manipulation), but here the formula was generalized in matrix form.

Therefore, if we restrict our study, for example, to a circular arc with a central angle of ${90}^{\circ }$, we may begin with a quadratic polynomial ($n=2$) associated with the weights $w=[0,\frac{\sqrt{2}}{2},1]$ (see, Ref. [16]). Subsequently, Eq. (3) can be applied to compute the weights resulting from successive degree elevations until a certain degree $n>2$.

For each of the abovementioned degree elevations, we can determine the transformation matrix ${\mathbf{T}}_{LB}$, and then implement Eq. (11) to determine the weights associated with the rational Lagrange polynomials.

Remark-2: Nodal points. By definition, non-uniform rational B-splines (NURBS) is a non-uniform interpolation. This means that even if the interval $[0,1]$ is uniformly subdivided by the points $\xi =[0,{\displaystyle \frac{1}{n}},\dots ,{\displaystyle \frac{(n-1)}{n}},1]$, the corresponding nodal points do not generally subdivide the curve $C\left(\xi \right)$ in equal arc lengths. This issue is shown later in Section 3.1.2.

2.3. Determination of Transformation Matrix

In this section, we present two alternative procedures for constructing the transformation matrix.

2.3.1. Detailed Algebraic Procedure

Since Bernstein and Lagrange polynomials span the same space (of monomials ${\xi }^i,i=0,1,\dots ,n$), any univariate function $U\left(\xi \right)$ will be equally written in any of them. By equating the coefficients of the same powers ${\xi }^i$, we can derive a relationship between nodal values $U_i$ and generalized coefficients $a_i$. This relationship is set in the form:

$${\mathbf{U}}_{\mathrm{nodal}}={\mathbf{T}}_{LB}^T\mathbf{a}.$$

(14)

Example: Let we set the function $U\left(\xi \right)$ in terms of Lagrange polynomials polynomials (based on uniform nodal points at $\xi =0,\frac{1}{2},1$), as follows:

$$\begin{array}{cc}\hfill U\left(\xi \right)& =L_{0,2}\left(\xi \right)U_0+L_{1,2}\left(\xi \right)U_1+L_{2,2}\left(\xi \right)U_2\hfill \\ \hfill & =(2{\xi }^2-3\xi +1)U_0+(-4{\xi }^2+4\xi )U_1+(2{\xi }^2-\xi )U_2\hfill \\ \hfill & ={\xi }^2(2U_0-4U_1+2U_2)+\xi (-3U_0+4U_1-U_2)+U_0.\hfill \end{array}$$

(15)

Alternatively, let us we set the same function $U\left(\xi \right)$ in terms of Bernstein polynomials:

$$\begin{array}{cc}\hfill U\left(\xi \right)& =B_{0,2}\left(\xi \right)a_0+B_{1,2}\left(\xi \right)a_1+B_{2,2}\left(\xi \right)a_2\hfill \\ \hfill & ={(1-\xi )}^2{a}_0+2(1-\xi )\xi a_1+\left({\xi }^2\right)a_2\hfill \\ \hfill & =({\xi }^2-2\xi +1)a_0+(-2{\xi }^2+2\xi )a_1+\left({\xi }^2\right)a_2\hfill \\ \hfill & ={\xi }^2(a_0-2a_1+a_2)+\xi (-2a_0+2a_1)+a_0.\hfill \end{array}$$

(16)

By equating the coefficients of the same powers between the quadratic polynomials given by Eq. (15) and Eq. (16), we receive the following system:

$$\begin{array}{cc}\hfill 2U_0-4U_1+2U_2& =a_0-2a_1+a_2\hfill \\ \hfill -3U_0+4U_1-U_2& =-2a_0+2a_1\hfill \\ \hfill U_0=a_0.\end{array}$$

(17)

Equivalently, Eq. (17) can be written as follows:

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}2& -4& 2\\ -3& 4& -1\\ 1& 0& 0\end{array}\right]\left[\begin{array}{c}{U}_0\\ U_1\\ U_2\end{array}\right]=\left[\begin{array}{ccc}1& -2& 1\\ -2& 2& 0\\ 1& 0& 0\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right].\end{array}$$

(18)

Equation (18) can be solved for either of the two vectors. Choosing the second one, after matrix inversion we have:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ -{\displaystyle \frac{1}{2}}& 2& -{\displaystyle \frac{1}{2}}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{U}_0\\ U_1\\ U_2\end{array}\right].\end{array}$$

(19)

Inserting the three coefficients from Eq. (19) into Eq. (16), we have:

$$\begin{array}{cc}\hfill U\left(\xi \right)& =B_{0,2}\left(\xi \right)a_0+B_{1,2}\left(\xi \right)a_1+B_{2,2}\left(\xi \right)a_2\hfill \\ \hfill & =B_{0,2}\left(\xi \right)U_0+B_{1,2}\left(\xi \right)(-\frac{1}{2}{U}_0+2U_1-\frac{1}{2}{U}_2)+B_{2,2}\left(\xi \right)U_2\hfill \\ \hfill & =\left[B_{0,2}\left(\xi \right)-\frac{1}{2}{B}_{1,2}\left(\xi \right)\right]U_0+\left[2B_{1,2}\left(\xi \right)\right]U_1+\left[-\frac{1}{2}{B}_{1,2}\left(\xi \right)+B_{2,2}\left(\xi \right)\right]U_2.\hfill \end{array}$$

(20)

Comparing Eq. (20) with the first equality of Eq. (15), the unique representation –in terms of ($U_0,U_1,U_2$)– leads to the relationships:

$$\begin{array}{cc}\hfill L_{0,2}\left(\xi \right)& =B_{0,2}\left(\xi \right)-\frac{1}{2}{B}_{1,2}\left(\xi \right)\hfill \\ \hfill L_{1,2}\left(\xi \right)& =2B_{1,2}\left(\xi \right)\hfill \\ \hfill L_{2,2}\left(\xi \right)& =-\frac{1}{2}{B}_{1,2}\left(\xi \right)+B_{2,2}\left(\xi \right).\hfill \end{array}$$

(21)

In matrix form, Eq. (21) is rewritten as:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{L}_{0,2}\left(\xi \right)\\ L_{1,2}\left(\xi \right)\\ L_{2,2}\left(\xi \right)\end{array}\right]=\underset{{\mathbf{T}}_{BL}}{\underbrace{\left[\begin{array}{ccc}1& -\frac{1}{2}& 0\\ 0& 2& 0\\ 0& -\frac{1}{2}& 1\end{array}\right]}}\left[\begin{array}{c}{B}_{0,2}\left(\xi \right)\\ B_{1,2}\left(\xi \right)\\ B_{2,2}\left(\xi \right)\end{array}\right].\end{array}$$

(22)

By inverting Eq. (22), we receive:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{B}_{0,2}\left(\xi \right)\\ B_{1,2}\left(\xi \right)\\ B_{2,2}\left(\xi \right)\end{array}\right]=\underset{{\mathbf{T}}_{LB}}{\underbrace{\left[\begin{array}{ccc}1& \frac{1}{4}& 0\\ 0& \frac{1}{2}& 0\\ 0& \frac{1}{4}& 1\end{array}\right]}}\left[\begin{array}{c}{L}_{0,2}\left(\xi \right)\\ L_{1,2}\left(\xi \right)\\ L_{2,2}\left(\xi \right)\end{array}\right].\end{array}$$

(23)

Therefore, if we pack all triplets into the following relevant column vectors:

$$\mathbf{U}=\left[\begin{array}{c}{U}_0\\ U_1\\ U_2\end{array}\right],\mathbf{a}=\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right],\mathbf{L}=\left[\begin{array}{c}{L}_{0,2}\left(\xi \right)\\ L_{1,2}\left(\xi \right)\\ L_{2,2}\left(\xi \right)\end{array}\right],\left[\begin{array}{c}{B}_{0,2}\left(\xi \right)\\ B_{1,2}\left(\xi \right)\\ B_{2,2}\left(\xi \right)\end{array}\right],$$

(24)

the above relationships are written in matrix form as follows:

$$\mathbf{L}={\mathbf{T}}_{BL}\mathbf{B},\mathbf{B}={\mathbf{T}}_{LB}\mathbf{L},$$

(25)

and

$$\mathbf{a}={\mathbf{T}}_{BL}^T\mathbf{U},\mathbf{U}={\mathbf{T}}_{LB}^T\mathbf{a},$$

(26)

with

$${\mathbf{T}}_{BL}={\left({\mathbf{T}}_{LB}\right)}^{-1}.$$

(27)

2.3.2. A Shorter Numerical Procedure

Since the transformation matrices merely relate the nodal vector $\mathbf{U}$ with the coefficients vector $\mathbf{a}$, in a unique way (because of the same polynomial degree), we can start by implementing the well known formula:

$$U\left(\xi \right)=\sum _{i=0}^2{B}_{i,2}\left(\xi \right)a_i,$$

(28)

at the uniform points $\xi =0,\frac{1}{2},1$ (with which the nodal values $U_0,U_1,U_2$ are associated), and thus we shortly receive:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{U}_0\\ U_1\\ U_2\end{array}\right]& =\left[\begin{array}{ccc}{B}_{0,2}\left(0\right)& B_{1,2}\left(0\right)& B_{2,2}\left(0\right)\\ B_{0,2}\left(\frac{1}{2}\right)& B_{1,2}\left(\frac{1}{2}\right)& B_{2,2}\left(\frac{1}{2}\right)\\ B_{0,2}\left(1\right)& B_{1,2}\left(1\right)& B_{2,2}\left(1\right)\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]\hfill \\ \hfill & =\underset{{\mathbf{T}}_{LB}^T}{\underbrace{\left[\begin{array}{ccc}1& 0& 0\\ \frac{1}{4}& \frac{1}{2}& \frac{1}{4}\\ 0& 0& 1\end{array}\right]}}\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right].\hfill \end{array}$$

(29)

Therefore, by applying the Bernstein-based series expansion Eq. (28) at three uniform parameters $\xi$, we can shortly obtain the transpose of the transformation matrix ${\mathbf{T}}_{LB}$, which denotes the change from Lagrange (L) to Bernstein (B) basis.

2.4. From Approximation Theory and Numerical Analysis to CAGD

The conclusions of Section 2.3 pertain to the relationship between the nodal values $U_i$ and the generalized coefficients $a_i$. Within the context of computer-aided geometric design (CAGD), the nodal values $U_i$ correspond to the Cartesian coordinates ${\mathbf{x}}_i$, whereas the generalized coefficients $a_i$ are interpreted as the control points ${\mathbf{x}}_i^c$.

In this context, in non-rational formulation, the substitute of Eq. (28) is

$$\mathbf{x}\left(\xi \right)=\sum _{i=0}^2{B}_{i,2}\left(\xi \right){\mathbf{x}}_i^c.$$

(30)

Therefore, when the weights equal unity ($w_i=1$), we can apply Eq. (30) at the nodal points ${\mathbf{x}}_i$ on the true curve (at specific–so far uniform–parameters $\xi$) and then solve in ${\mathbf{x}}_i^c$. The matrix which relates the vectors $\mathbf{x}$ and ${\mathbf{x}}^c$ is the transformation matrix under demand, i.e., $\mathbf{x}={\mathbf{T}}_{LB}^T{\mathbf{x}}^c$.

In the general case of non-unity weights ($w_i\ne 1$), Eq. (30) converts to:

$$\mathbf{x}\left(\xi \right)={\displaystyle \frac{{\sum }_{i=0}^2{B}_{i,2}\left(\xi \right){\mathbf{x}}_{i,w}}{{\sum }_{i=0}^2{B}_{i,2}\left(\xi \right)w_i}}={\displaystyle \frac{\left[B_{0,2}\left(\xi \right),B_{1,2}\left(\xi \right),B_{2,2}\left(\xi \right)\right]\left[\begin{array}{c}{\mathbf{x}}_{0,w}\\ {\mathbf{x}}_{1,w}\\ {\mathbf{x}}_{2,w}\end{array}\right]}{{\sum }_{i=0}^2{B}_{i,2}\left(\xi \right)w_i}},$$

(31)

where the projected (homogeneous) control points are given by the well-known formula:

$${\mathbf{x}}_{\mathrm{w}}=\left[\begin{array}{c}{\mathbf{x}}_{0,w}\\ {\mathbf{x}}_{1,w}\\ {\mathbf{x}}_{2,w}\end{array}\right]\equiv \left[\begin{array}{c}{w}_0{\mathbf{x}}_0\\ w_1{\mathbf{x}}_1\\ w_2{\mathbf{x}}_2\end{array}\right].$$

(32)

Depending on the value of parameter $\xi =\{0,\frac{1}{2},1\}$, the row vector in the numerator of Eq. (31) is one out of the three rows of the matrix ${\mathbf{T}}_{LB}^T$ involved in Eq. (28). In more detail, considering Eq. (32), we have:

$$\begin{array}{cc}\hfill {\mathbf{x}}_0& ={\displaystyle \frac{\left[\begin{array}{ccc}{B}_{0,2}\left(0\right)& B_{1,2}\left(0\right)& B_{2,2}\left(0\right)\end{array}\right]{\mathbf{x}}_{\mathrm{w}}}{{\sum }_{i=0}^2{B}_{i,2}\left(0\right)w_i}}\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\mathbf{x}}_1& ={\displaystyle \frac{\left[\begin{array}{ccc}{B}_{0,2}\left(\frac{1}{2}\right)& B_{1,2}\left(\frac{1}{2}\right)& B_{2,2}\left(\frac{1}{2}\right)\end{array}\right]{\mathbf{x}}_{\mathrm{w}}}{{\sum }_{i=0}^2{B}_{i,2}\left(\frac{1}{2}\right)w_i}}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\mathbf{x}}_2& ={\displaystyle \frac{\left[\begin{array}{ccc}{B}_{0,2}\left(1\right)& B_{1,2}\left(1\right)& B_{2,2}\left(1\right)\end{array}\right]{\mathbf{x}}_{\mathrm{w}}}{{\sum }_{i=0}^2{B}_{i,2}\left(1\right)w_i}}.\hfill \end{array}$$

(35)

The totality of Eqs. (33) to (35) makes the column vector of nodal points on the curve, at the uniform values of parameter $\xi$, and thus we can write:

$${\mathbf{x}}_{\mathrm{nodal}}={\mathbf{T}}_{LB}^T{\mathbf{x}}_{\mathrm{w}},$$

(36)

where

$${\mathbf{x}}_{\mathrm{nodal}}=\left[\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{x}}_1\\ {\mathbf{x}}_2\end{array}\right].$$

(37)

Below we shall see how the nodal coordinates ${\mathbf{x}}_{\mathrm{nodal}}$ are introduced into the formulation. Actually, we return to Eq. (31) and substitute the row vector ${\mathbf{B}}^T$ of Bernstein polynomials by the row vector ${\mathbf{L}}^T$ of Lagrange polynomials using the second equality of Eq. (25). Furthermore, by virtue of Eq. (36), we have:

$$\mathbf{x}\left(\xi \right)={\displaystyle \frac{\left[\begin{array}{ccc}{L}_{0,2}\left(\frac{1}{2}\right)& L_{1,2}\left(\frac{1}{2}\right)& L_{2,2}\left(\frac{1}{2}\right)\end{array}\right]{\mathbf{T}}_{LB}^T{\mathbf{x}}_w}{{\sum }_{i=0}^2{B}_{i,2}\left(\xi \right)w_i}},$$

(38)

Regarding Eq. (38), one may observe that the weight in the numerator is incorporated into the projected (homogeneous) coordinate ${\mathbf{x}}_{i,w}$ of the i-th control point ${\mathbf{P}}_i$, whereas in the denominator appears as a factor of the Bernstein polynomial associated with current control point.

It is worth-mentioning that although the three nodal points (${\mathbf{x}}_0,{\mathbf{x}}_1,{\mathbf{x}}_2$) subdivide the circular arc into two equal parts, due to the weights $w_i$, in general, even if each part is further subdivided into uniform increments $\Delta \xi$ the mapping from the prameter space to the natural space is non-uniform. This means that the increments $\Delta \alpha$ of the corresponding central angles will be non-uniform, as wll be later discussed in Section 3.1.2.

The determination of the transformation matrix, ${\mathbf{T}}_{LB}^T$, is discussed in Section 3.

3. Weights

3.1. Circular Arc with Central Angle ${90}^{\circ }$

3.1.1. Uniform Distribution

We consider a circular arc of central angle ${90}^{\circ }$ and radius $R=1$, in the first quadrant. We start with degree $p=2$ in conjunction with the well-known control points ${\mathbf{P}}_0(1,0)$, ${\mathbf{P}}_1(1,1)$, and ${\mathbf{P}}_2(0,1)$, and the initial weights ($w_0=1,w_1=\frac{\sqrt{2}}{2},w_2=1$, according to [16]).

As a first step, we continuously elevate the polynomial degree applying Eq. (3), thus determining the updated weights for each degree. Furthermore, applying Eq. (1), we determine the Lagrange-based weights.

The second step is to determine the transpose of the transformation matrix ${\mathbf{T}}_{LB}$, for the desired degree p. Collectively, for degrees $p=2$ until $p=6$, we have:

$$\begin{array}{cc}\hfill {\left({\mathbf{T}}_{LB}^T\right)}_{p=2}& =\left[\begin{array}{ccc}1& 0& 0\\ 1/4& 1/2& 1/4\\ 0& 0& 1\end{array}\right]\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\left({\mathbf{T}}_{LB}^T\right)}_{p=3}& =\left[\begin{array}{cccc}1& 0& 0& 0\\ 8/27& 4/9& 2/9& 1/27\\ 1/27& 2/9& 4/9& 8/27\\ 0& 0& 0& 1\end{array}\right]\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\left({\mathbf{T}}_{LB}^T\right)}_{p=4}& =\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 81/256& 27/64& 27/128& 3/64& 1/256\\ 1/16& 1/4& 3/8& 1/4& 1/16\\ 1/256& 3/64& 27/128& 27/64& 81/256\\ 0& 0& 0& 0& 1\end{array}\right]\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\left({\mathbf{T}}_{LB}^T\right)}_{p=5}& =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 541/1651& 256/625& 128/625& 32/625& 4/625& 1/3125\\ 243/3125& 162/625& 216/625& 144/625& 48/625& 32/3125\\ 32/3125& 48/625& 144/625& 216/625& 162/625& 243/3125\\ 1/3125& 4/625& 32/625& 128/625& 256/625& 541/1651\\ 0& 0& 0& 0& 0& 1\end{array}\right]\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\left({\mathbf{T}}_{LB}^T\right)}_{p=6}& =\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 214/639& 428/1065& 214/1065& 234/4367& 125/15552& 5/7776& 1/46656\\ 64/729& 64/243& 80/243& 160/729& 20/243& 4/243& 1/729\\ 1/64& 3/32& 15/64& 5/16& 15/64& 3/32& 1/64\\ 1/729& 4/243& 20/243& 160/729& 80/243& 64/243& 64/729\\ 1/46656& 5/7776& 125/15552& 234/4367& 214/1065& 428/1065& 214/639\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right].\hfill \end{array}$$

(43)

Regarding the associated weights, these are shown in Table 1.

3.1.2. Non-Uniform Distribution

Now, we consider the case of non-uniform distribution of nodal points along the edges of a 27-node classical transfinite element (shown in Figure 1), where the nodes along the horizontal edge AB are located at the non-uniform positions:

$$\xi =\{0,1/6,2/6,3/6,3/4,1\}.$$

(44)

The formulation is the same, with the only difference regarding the ${\left({\mathbf{T}}_{LB}^T\right)}_{p=5}^{\mathrm{non}}-\mathrm{uniform}$, which now becomes:

$$\begin{array}{cc}\hfill {\left({\mathbf{T}}_{LB}^T\right)}_{p=5}^{\mathrm{non}}-\mathrm{uniform}& =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 428/1065& 428/1065& 625/3888& 125/3888& 25/7776& 1/7776\\ 32/243& 80/243& 80/243& 40/243& 10/243& 1/243\\ 1/32& 5/32& 5/16& 5/16& 5/32& 1/32\\ 1/1024& 15/1024& 45/512& 135/512& 405/1024& 243/1024\\ 0& 0& 0& 0& 0& 1\end{array}\right],\hfill \end{array}$$

(45)

whence the updated weight vector (according to Eq. (11)) becomes:

$${\mathbf{w}}_L=\left[\begin{array}{cccccc}1& {\displaystyle \frac{26+5\sqrt{2}}{36}}& {\displaystyle \frac{5+2\sqrt{2}}{9}}& {\displaystyle \frac{2+\sqrt{2}}{4}}& {\displaystyle \frac{10+3\sqrt{2}}{16}}& 1\end{array}\right].$$

(46)

One may obseve in Eq. (46), that the non-uniform node distribution (Eq. (44)) leads to the non-symmetric form of vector ${\mathbf{w}}_L$.

In the following, we shall demonstrate that both a straight segment and a circular arc can be represented accurately.

I. Straight segment

Taking the uniform set of control points associated with the Bernstein basis of degree $p=5$ (considering weights ${\mathbf{w}}_B=1$):

$${\mathbf{P}}_{\mathrm{B}}^w={\mathbf{P}}_{\mathrm{B}}={\left[0,{\displaystyle \frac{1}{5}},{\displaystyle \frac{2}{5}},{\displaystyle \frac{3}{5}},{\displaystyle \frac{4}{5}},1\right]}^T,$$

(47)

and left-multiplying by the non-uniform transformation matrix given by in Eq.(45) according to Eq. (10), we derive the true non-uniform coordinates:

$${\mathbf{U}}_{\mathrm{nodal}}={\left({\mathbf{T}}_{LB}^T\right)}_{p=5}^{\mathrm{non}}-\mathrm{uniform}·{\mathbf{P}}_{\mathrm{B}}^w={\left[0,{\displaystyle \frac{1}{6}},{\displaystyle \frac{2}{6}},{\displaystyle \frac{3}{6}},{\displaystyle \frac{3}{4}},1\right]}^T.$$

(48)

The result of Eq. (48) is validated by the non-uniform nodal points shown on the edge AB of Figure 1.

Remark: In contrast, if the uniform set of Eq. (47) had been multiplied by the uniform matrix of Eq. (42), the location of the nodal points on the straight segment would be uniform as well.

II. Circular Arc

Concerning the projected coordinates in Lagrange formulation, they can be calculated through Eq. (10), using the Cartesian coordinates of control points in Bernstein formulation (Table 2) in combination with the associated weights shown in Table 1. Afterwards, the projected coordinates in Lagrange formulation are divided by the corresponding entries of vector ${\mathbf{w}}_L$ given by Eq. (46).

The results of the above procedure is illustrated in Figure 2, where the nodal sequence {1,2,3,4,5,6} of Lagrangian formulation is in the counter-clockwise direction. The five central angles formed by the radii associated with the ends of each nodal segment (i.e., 1-2,…,5-6), are shown in Table 2. One may observe that nodal point 4 (at $\xi =0.50$) is exactly in the middle of the circular arc (i.e., $a_1+a_2+a_3={45}^{\circ }$). However, the other nodes do not accurately split the two parts in three (nodes 1 to 4) and two (nodes 4 to 6) segments, respectively, of equal magnitude. In other words, the increments of the parameter $\Delta {\xi }_i$ differ from the increments of central angle ${\alpha }_i$ of Lagrangian formulation, as also happens with the identical central angle ${\beta }_i$ in the Bernstein formulation (with control points $P_0$ to $P_5$ illustrated in magenta color).

Finally, the rational and non-rational shape functions associated with the nodal points 1-to-6 on the circular arc are shown in Figure 3, where the same color has been used for corresponding curves.

3.2. Central Angle Different than ${90}^{\circ }$

The control triangle of the arc must be isosceles, the end weights can be set to 1, and the middle weight is the cosine of the angle between the chord joining the endpoints and one of the legs [16].

Therefore, if the central angle of the circular arc equals to $\alpha$, to accurately represent the geometry, one choice of the weights is:

$$w_0=1,w_1=cos\left({\displaystyle \frac{\alpha }{2}}\right),w_2=1.$$

(49)

However, the end weights do not have to remain 1; the arc can be "reweighted" using the shape invariace ratio:

$$\lambda ={\displaystyle \frac{w_0{w}_2}{w_1^2}}.$$

(50)

For example, an equivalence choice to that of Eq. (49) may be:

$$w_0=2,w_1=2cos\left({\displaystyle \frac{\alpha }{2}}\right),w_2=2.$$

(51)

4. Examples

4.1. Example-1: The Case of Degree p=2

Let us consider a circular arc in the first quadrant, with a central angle of ${90}^{\circ }$. Considering the transformation matrix given by Eq. (39) as well as the well-known vector of weights ${\mathbf{w}}_B=[1,{\displaystyle \frac{\sqrt{2}}{2}},1]$ (Ref. [16]), the Lagrangian weights are given by the vector (see Eq. (11)):

$$\begin{array}{cc}\hfill {\mathbf{w}}_L& ={\mathbf{T}}_{LB}^T{\mathbf{w}}_B\hfill \\ \hfill & =\left[\begin{array}{ccc}1& 0& 0\\ 1/4& 1/2& 1/4\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}1\\ {\displaystyle \frac{\sqrt{2}}{2}}\\ 1\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}1\\ ({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\sqrt{2}}{4}})\\ 1\end{array}\right]\hfill \end{array}$$

(52)

Moreover, Eq. (10) depicts that the nodal points on the circular arc are given as:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_L^w& ={\mathbf{T}}_{LB}^T{\mathbf{P}}_B^w\hfill \\ \hfill & =\left[\begin{array}{ccc}1& 0& 0\\ 1/4& 1/2& 1/4\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cc}1·1& 0·1\\ 1·{\displaystyle \frac{\sqrt{2}}{2}}& 1·{\displaystyle \frac{\sqrt{2}}{2}}\\ 0·1& 1·1\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{cc}1& 0\\ {\displaystyle \frac{1+\sqrt{2}}{4}}& {\displaystyle \frac{1+\sqrt{2}}{4}}\\ 0& 1\end{array}\right]\hfill \end{array}$$

(53)

To obtain the Cartesian coordinates of the nodal points onto the circular arc, the entries of matrix ${\mathbf{Q}}_L^w$ (Eq. (53)) must be divided by the corresponding elements of vector ${\mathbf{w}}_L$ (Eq. (52)), thus obtaining:

$$\begin{array}{cc}\hfill {\mathbf{Q}}_L& =\left[\begin{array}{cc}1& 0\\ {\displaystyle \frac{1+\sqrt{2}}{4}}:({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\sqrt{2}}{4}})& {\displaystyle \frac{1+\sqrt{2}}{4}}:({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\sqrt{2}}{4}})\\ 0& 1\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{cc}1& 0\\ {\displaystyle \frac{\sqrt{2}}{2}}& {\displaystyle \frac{\sqrt{2}}{2}}\\ 0& 1\end{array}\right].\hfill \end{array}$$

(54)

Therefore, in Eq. (54) we eventually derived the two endpoints plus the middle point of the circular arc.

Remark: In the above case where $p=2$, the nodal points on which the equivalent rational Lagrange polynomials are based, have divided the circular arc in two equal parts, i.e., $[0,{\displaystyle \frac{1}{2}}]$ and $[{\displaystyle \frac{1}{2}},1]$.

Nevertheless, for higher degrees (i.e., $p>2$) the nodal points are not uniformly distributed. Since the parameterization does not change with degree elevation [2], the findings of this section for $p=2$ are sufficient to determine the position of the nodal points for higher degrees.

Within this context, for any polynomial degree $p>2$ which is related to parameter values $\xi =[0,{\displaystyle \frac{1}{p}},\dots ,{\displaystyle \frac{p-1}{p}},1]$, the nodal points are calculated through the quadratic model of this subsection as follows:

$$\begin{array}{cc}\hfill \mathbf{x}\left(\xi \right)& ={\displaystyle \frac{B_{0,2}\left(\xi \right)}{B_{0,2}\left(\xi \right)+\frac{\sqrt{2}}{2}{B}_{1,2}\left(\xi \right)+B_{2,2}\left(\xi \right)}}\left[\begin{array}{c}1\\ 0\end{array}\right]\hfill \\ \hfill & +{\displaystyle \frac{\frac{\sqrt{2}}{2}{B}_{1,2}\left(\xi \right)}{B_{0,2}\left(\xi \right)+\frac{\sqrt{2}}{2}{B}_{1,2}\left(\xi \right)+B_{2,2}\left(\xi \right)}}\left[\begin{array}{c}1\\ 1\end{array}\right]\hfill \\ \hfill & +{\displaystyle \frac{B_{2,2}\left(\xi \right)}{B_{0,2}\left(\xi \right)+\frac{\sqrt{2}}{2}{B}_{1,2}\left(\xi \right)+B_{2,2}\left(\xi \right)}}\left[\begin{array}{c}0\\ 1\end{array}\right].\hfill \end{array}$$

(55)

Sweeping the spectrum of polynomial degrees from $p=2$ to $p=6$, using Eq. (55) the results are tabulated in Table 3.

4.2. Example-2: The Case of Degree p=3

4.2.1. Straight Line

First, it is worthy to comment on the representation of a straight line through rational interpolation. Within this context, let us consider a segment AB, made of points A(0,0) and B(1,1), which is uniformly subdivided in three equal parts by points C and D, as shown in Figure 4. Assuming that the parametric space along AB corresponds to the interval $\xi \in [0,1]$, it is reasonable to consider that points C and D correspond to $\xi ={\displaystyle \frac{1}{3}}$ and $\xi ={\displaystyle \frac{2}{3}}$, respectively. Therefore, the associated Lagrange polynomials are given by:

$$\begin{array}{cc}\hfill L_0\left(\xi \right)& ={\displaystyle \frac{(\xi -1/3)(\xi -2/3)(\xi -1)}{(0-1/3)(0-2/3)(0-1)}}=-9(\xi -1/3)(\xi -2/3)(\xi -1)\hfill \\ \hfill L_1\left(\xi \right)& ={\displaystyle \frac{(\xi -0)(\xi -2/3)(\xi -1)}{(1/3-0)(1/3-2/3)(1/3-1)}}={\displaystyle \frac{27}{2}}\xi (\xi -2/3)(\xi -1)\hfill \\ \hfill L_2\left(\xi \right)& ={\displaystyle \frac{(\xi -0)(\xi -1/3)(\xi -1)}{(2/3-0)(2/3-1/3)(2/3-1)}}=-{\displaystyle \frac{27}{2}}\xi (\xi -1/3)(\xi -1)\hfill \\ \hfill L_3\left(\xi \right)& ={\displaystyle \frac{(\xi -0)(\xi -1/3)(\xi -2/3)}{(1-0)(1-1/3)(1-2/3)}}={\displaystyle \frac{9}{2}}\xi (\xi -1/3)(\xi -2/3).\hfill \end{array}$$

(56)

Given the weights

$$w_0=1,w_1,w_2,w_3=1,$$

(57)

associated with the points (A, C, D, B), the rational approximation is given as:

$$x\left(\xi \right)={\displaystyle \frac{w_0{L}_0{X}_0+w_1{L}_1{X}_1+w_2{L}_2{X}_2+w_3{L}_3{X}_3}{w_0{L}_0+w_1{L}_1+w_2{L}_2+w_3{L}_3}}$$

(58)

and

$$y\left(\xi \right)={\displaystyle \frac{w_0{L}_0{Y}_0+w_1{L}_1{Y}_1+w_2{L}_2{Y}_2+w_3{L}_3{Y}_3}{w_0{L}_0+w_1{L}_1+w_2{L}_2+w_3{L}_3}}$$

(59)

By definition, due to the cardinality of rational shape functions, the parameters ${\xi }_i=0,{\displaystyle \frac{1}{3}},{\displaystyle \frac{2}{3}},1$ correspond to the points A, C, D, and B, respectively. In all cases, it was found that the whole mapping is the straight segment AB itself. For example, if the segment lies on the bisector of the right angle in the first quadrant, the nodal points are equal each other, and thus the same rational basis functions lead to the same output ($x\left(\xi \right)=y\left(\xi \right)$).

Moreover:

• When $w_1=w_2$, the mapping is uniform, i.e., $x\left(\xi \right)=y\left(\xi \right)=\xi$.
• In contrast, when $w_1\ne w_2$, the mapping is non-uniform.

4.2.2. Circular Arc of ${90}^{\circ }$

Now, we focus on the representation of the circular arc, using rational Lagrange polynomials of degree $p=3$.

While Eq. (56) is still valid, in current case–due to symmetry with respect to the bisector of first quadrant–the weights are taken as

$$w_0=1,w_1=w_2,w_3=1,$$

(60)

The question is whether it is possible to determine a unique value of the weight $w_1$ such that the corresponding rational Lagrange polynomials, based on the uniform nodes (A, C, D, B) in the physical space, can exactly represent the circular arc. In this regard, by means of a counterexample, we will show that this is not generally possible.

Actually, substituting Eq. (58) and Eq. (59) into the desired identity:

$$x{\left(\xi \right)}^2+y{\left(\xi \right)}^2=1,$$

(61)

we derive a polynomial of sixth degree in $\xi$, as follows:

$$p\left(\xi \right)=\sum _{i=0}^6{c}_i{\xi }^i\equiv 0,$$

(62)

with

$$\begin{array}{cc}\hfill c_0& =0,\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill c_1& =\left(9\sqrt{3}-{\displaystyle \frac{27}{2}}\right)w_1-2,\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill c_2& ={\displaystyle \frac{(162-81\sqrt{3})}{2}}{w}_1^2+{\displaystyle \frac{(405-306\sqrt{3})}{4}}{w}_1+20,\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill c_3& ={\displaystyle \frac{(1053\sqrt{3}-2106)}{4}}{w}_1^2+{\displaystyle \frac{(513\sqrt{3}-594)}{2}}{w}_1-{\displaystyle \frac{153}{2}},\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill c_4& ={\displaystyle \frac{(4698-2349\sqrt{3})}{4}}{w}_1^2+{\displaystyle \frac{(1809-1728\sqrt{3})}{4}}{w}_1+{\displaystyle \frac{279}{2}},\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill c_5& ={\displaystyle \frac{243}{4}}\left[9(\sqrt{3}-2)w_1^2+18(\sqrt{3}-1)w_1-6\right],\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill c_6& ={\displaystyle \frac{1}{4}}\left[729(2-\sqrt{3})w_1^2+486(1-\sqrt{3})w_1+162\right].\hfill \end{array}$$

(69)

One may observe that the constant term $c_0$ a priori vanishes (Eq. (63)), while the root of Eq. (64) provides the value ${\left(w_1\right)}_1\approx 0.9576447$. Moreover, the closest root of the binomial given by Eq. (65) is slightly different as it equals to ${\left(w_1\right)}_2\approx 0.9599365$, whereas Eq. (66) vanishes at a higher value ${\left(w_1\right)}_3\approx 0.9715457$. Furthermore, Eq. (67) has two complex roots, and the rest two equations substantially differ from ${\left(w_1\right)}_1\approx 0.9576447$. As a result, there is no real value $w_1$ which makes all the coefficients $c_0$ to $c_6$ be zero.

In conclusion:

• It is not generally possible to accurately represent a circular arc using equi-distant nodal points, even if rational Lagrange polynomials are used. Based on optimization techniques, where an objective function such as $f\left(w_1\right)={(x^2+y^2-1)}^2$ becomes minimum, it is possible to determine a value $w_1$ which compromizes the conditions $c_i\approx 0,(i=1,\dots ,6)$. Nevertheless, small deviations from the ideal circle will be noticed (e.g., after the fourth decimal point).
• In contrast, if the nodal points i the physical space are taken at those non-equi-distant positions shown in Table 3, whereas the associated weights are taken according to Table 1, then the equation of the circle ($x^2+y^2=1$) is identically satisfied.

4.3. Example-3: The Case of Degree p=5

Let us deal with the uniform case of $p=5$. For $0\le \xi \le 1$, the Lagrange polynomials are as follows:

$$\begin{array}{cc}\hfill L_1\left(\xi \right)& =-\left(625\right(\xi -1\left)\right(\xi -1/5\left)\right(\xi -2/5\left)\right(\xi -3/5\left)\right(\xi -4/5\left)\right)/24\hfill \\ \hfill & =-\left(625{\xi }^5\right)/24+\left(625{\xi }^4\right)/8-\left(2125{\xi }^3\right)/24+\left(375{\xi }^2\right)/8-\left(137\xi \right)/12+1,\hfill \\ \hfill L_2\left(\xi \right)& =3125\xi (\xi -1)(\xi -2/5)(\xi -3/5)(\xi -4/5))/24\hfill \\ \hfill & =\left(3125{\xi }^5\right)/24-\left(4375{\xi }^4\right)/12+\left(8875{\xi }^3\right)/24-\left(1925{\xi }^2\right)/12+25\xi ,\hfill \\ \hfill L_3\left(\xi \right)& =3125\xi (\xi -1)(\xi -1/5)(\xi -3/5)(\xi -4/5))/12\hfill \\ \hfill & =-\left(3125{\xi }^5\right)/12+\left(8125{\xi }^4\right)/12-\left(7375{\xi }^3\right)/12+\left(2675{\xi }^2\right)/12-25\xi ,\hfill \\ \hfill L_4\left(\xi \right)& =\left(3125\xi \right(\xi -1\left)\right(\xi -1/5\left)\right(\xi -2/5\left)\right(\xi -4/5\left)\right)/12\hfill \\ \hfill & =\left(3125{\xi }^5\right)/12-625{\xi }^4+\left(6125{\xi }^3\right)/12-\left(325{\xi }^2\right)/2+\left(50\xi \right)/3,\hfill \\ \hfill L_5\left(\xi \right)& =-\left(3125\xi \right(\xi -1\left)\right(\xi -1/5\left)\right(\xi -2/5\left)\right(\xi -3/5\left)\right)/24\hfill \\ \hfill & =-\left(3125{\xi }^5\right)/24+\left(6875{\xi }^4\right)/24-\left(5125{\xi }^3\right)/24+\left(1525{\xi }^2\right)/24-\left(25\xi \right)/4,\hfill \\ \hfill L_6\left(\xi \right)& =\left(625\xi \right(\xi -1/5\left)\right(\xi -2/5\left)\right(\xi -3/5\left)\right(\xi -4/5\left)\right)/24\hfill \\ \hfill & =\left(625{\xi }^5\right)/24-\left(625{\xi }^4\right)/12+\left(875{\xi }^3\right)/24-\left(125{\xi }^2\right)/12+\xi .\hfill \end{array}$$

(70)

The Cartesian coordinates of control points in the Bernstein model are as follows:

$$\mathbf{Q}=\left[\begin{array}{cccccc}1& 1& {\displaystyle \frac{3+3\sqrt{2}}{4+3\sqrt{2}}}& {\displaystyle \frac{1+3\sqrt{2}}{4+3\sqrt{2}}}& {\displaystyle \frac{\sqrt{2}}{3+\sqrt{2}}}& 0\\ 0& {\displaystyle \frac{\sqrt{2}}{3+\sqrt{2}}}& {\displaystyle \frac{1+3\sqrt{2}}{4+3\sqrt{2}}}& {\displaystyle \frac{3+3\sqrt{2}}{4+3\sqrt{2}}}& 1& 1\end{array}\right]$$

(71)

The circular arc in the first quadrant is subdivided at the nodal points which correspond to the parametric vector $\Xi =\left\{0,{\displaystyle \frac{1}{5}},{\displaystyle \frac{2}{5}},{\displaystyle \frac{3}{5}},{\displaystyle \frac{4}{5}},1\right\}$. To determine them, it becomes necessary to apply the rational Bézier formula:

$$\mathbf{x}\left(\xi \right)={\displaystyle \frac{{\sum }_{j=0}^5{B}_{j,5}\left(\xi \right)w_j{\mathbf{Q}}_j}{{\sum }_{j=0}^5{B}_{j,5}\left(\xi \right)w_j}},\mathrm{with}\xi \in \Xi =\{0,\frac{1}{5},\frac{2}{5},\frac{3}{5},\frac{4}{5},1\},$$

(72)

where the coordinates ${\mathbf{Q}}_i$ are given by Eq. (71), while the corresponding weights are isolated in the following row vector:

$${\mathbf{w}}_b=\left[1,{\displaystyle \frac{3+\sqrt{2}}{5}},{\displaystyle \frac{4+3\sqrt{2}}{10}},{\displaystyle \frac{4+3\sqrt{2}}{10}},{\displaystyle \frac{3+\sqrt{2}}{5}},1\right].$$

(73)

Based on the Cartesian coordinates of nodal points on the true circular arc according to Eq. (72), and using the weights shown in Table 1, repeated below:

$${\mathbf{w}}_L=\left[1,{\displaystyle \frac{3+\sqrt{2}}{5}},{\displaystyle \frac{4+3\sqrt{2}}{10}},{\displaystyle \frac{4+3\sqrt{2}}{10}},{\displaystyle \frac{3+\sqrt{2}}{5}},1\right]???,$$

(74)

the parametric equation of the circular arc is as follows:

$$\mathbf{x}\left(\xi \right)={\displaystyle \frac{{\sum }_{i=0}^5{L}_{i,5}\left(\xi \right)w_{Li}{\mathbf{Q}}_i}{{\sum }_{i=0}^5{L}_{i,5}\left(\xi \right)w_{Li}}},\mathrm{with}\xi \in \left[0,1\right],$$

(75)

The interesting issue is the position of the nodal points in the Lagrangian model. While for quadratic interpolation (i.e., $p=2$) the nodal points are the two endpoints and the middle point, in all other situations they are arranged in a non-uniform way. For example, for current case where $p=5$, the polar angle for the nodal points (1 to 6 in counter-clockwise direction) in the first quadrant is shown in Table 4. One may observe that, although $\Delta \xi ={\displaystyle \frac{1}{5}}$ is a constant, the corresponding central angles $\Delta \alpha$ differ one another.

5. Spherical Cap

The above formulation is applicable to surfaces as well. As an example, we consider a spherical cap, which is exactly one-sixth of a unit sphere’s surface. It is known that the quadratic interpolation is not capable of accurately representing this curvilinear surface [17]. In contrast, quartic interpolation ($p=4$) in conjunction with 25 control points can accurately represent this cap [18] (see, the parametric space in Figure 5).

Regarding the southern spherical cap, the location of 25 control points and the associated weights are available through the MATLAB function southcap=rsmak(‘southcap’). With respect to Bernstein polynomials of degree $p=4$, the Cartesian coordinates and the associated weights of control points are shown in Table 5. There, one may observe that the analytical formula is not an ideal tensor product, in the sense that the weights associated with the internal control points are not products of the corresponding boundary ones (i.e., $w_{ij}\ne w_i{w}_j$).

By analogy to the one-dimensional procedure reported in the previous sections of this paper, the algorithm to derive the nodal points (on the true sperical surface) and the associated weights is as follows:

• Calculate the projected control points ${\mathbf{P}}_B^w$ and the associated weights ${\mathbf{w}}_B$. This is implemented using the abovementioned built-in MATLAB function southcap=rsmak(‘southcap’).
• Create the total transformation matrix ${\mathbf{T}}_{2D}$, of size $25\times 25$ (see below). In more detail, if the one-dimensional matrix ${\left({\mathbf{T}}_{LB}^T\right)}_{p=4}$ of Eq. (41) is shortly written as:

  $${\mathbf{T}}_{1D}=\left[\begin{array}{ccc}{t}_{11}& \dots & t_{15}\\ ⋮& \dots & ⋮\\ t_{51}& \dots & t_{55}\end{array}\right],$$

  (76)

  the 2D analog–applied to a surface patch–becomes:

  $${\mathbf{T}}_{2D}=\left[\begin{array}{ccc}{t}_{11}{\mathbf{T}}_{1D}& \dots & t_{15}{\mathbf{T}}_{1D}\\ ⋮& \dots & ⋮\\ t_{51}{\mathbf{T}}_{1D}& \dots & t_{55}{\mathbf{T}}_{1D}\end{array}\right].$$

  (77)
• Multiply the control points ${\mathbf{P}}_B^w$ by the abovementioned ${\mathbf{T}}_{2D}$, and thus determine the projected nodal points ${\mathbf{Q}}_L^w$:

  $${\mathbf{Q}}_L^w={\mathbf{T}}_{2D}{\mathbf{P}}_B^w.$$

  (78)
• Multiply the weights vector ${\mathbf{w}}_B$ by the abovementioned ${\mathbf{T}}_{2D}$, and thus derive the weights ${\mathbf{w}}_L$ associated with Lagrange polynomials:

  $${\mathbf{w}}_L={\mathbf{T}}_{2D}{\mathbf{w}}_B.$$

  (79)
• At the end, divide the projected coordinates ${\mathbf{Q}}_L^w$ (Eq. (78)) by the corresponding entries of ${\mathbf{w}}_L$ (Eq. (79)), one-to-one, and thus derive the Cartesian coordinates ${\mathbf{Q}}_L$ of nodal points on the true spherical cap. Obviously, the two vectors (${\mathbf{Q}}_L$, ${\mathbf{w}}_L$) is the final outcome of the proposed algorithm.

The implementation of the above algorithm to the southern spherical cap of a sphere of unit radius ($R=1$) gives the numerical values shown in Table 6. Based on these data, for every pair $(\xi ,\eta )\in {[0,1]}^2$ the image was found to belong to the true spherical cap (numerical error less than ${10}^{-15}$), and also coincided with the corresponding point obtained using rational Bernstein plynomials.

6. A Boundary-Value Problem

In this section we study the numerical solution in the quarter of an annulus, with internal radius $R_1=1$ and external radius $R_2=32$. The corresponding temperatures on the circular edges are $T_1={1000}^{\circ }C$ and $T_2=0^{\circ }C$, whereas the flux perpendicular to the straight edges vanish ($\partial U/\partial n=0$), as shown in Figure 6.

In terms of the radius r (with $R_1\le r\le R_2$), the analytcal solution is given by:

$$T_{\mathrm{exact}}\left(r\right)=T_1+{\displaystyle \frac{(T_2-T_1)}{log\left(\frac{R_2}{R_1}\right)}}log\left({\displaystyle \frac{r}{R_1}}\right).$$

(80)

Numerical results are presented in terms of the $L_2$ error norm in percent, which is defined as follows:

$$L_2={\displaystyle \frac{{\int }_{\Omega }{\left(T_{\mathrm{calculated}}-T_{\mathrm{exact}}\right)}^{2}d\Omega }{{\int }_{\Omega }{\left(T_{\mathrm{exact}}\right)}^{2}d\Omega }}\times 100(\%).$$

(81)

Each of the two circular arcs was divided into two equal parts using 3 nodes on the true boundary, whereas the weights associated with quadratic Lagrange polynomials were ${\mathbf{w}}_L=[1,{\displaystyle \frac{1}{2}}+{\displaystyle \frac{\sqrt{2}}{4}},1]$. The number $n_y$ of uniform sudivisions in the radial direction varied among 1 and 25 (i.e., degree $1\le p_y\le 23$), and relevant Lagrange polynomials were combined with unit weights. For any given value $n_y$, the total number of nodes in the model, on the true boundary and the interior, is $n_{\mathrm{nodes}}=3(n_y+1)$.

The calculated $L_2$ error norm (in percent) is shown in Figure 7. It may be observed that up to $p_y=23$ the convergence is monotonic, after which it (not shown) increases slightly. It is also worth noting that even for the coarsest model (i.e., $n_y=1$), the exact analytical value of the area

$$A={\displaystyle \frac{\pi \left(R_2^2-R_1^2\right)}{4}}$$

was accurately approximated by numerical performing the integral of the determinant of Jacobian matrix, i.e.,

$$A=\int dA=\int det\left(J\right)d\xi d\eta .$$

7. Discussion

It was shown that the set of weights associated with Bernstein polynomials and control points, which determine a certain curve $\mathbf{C}\left(\xi \right)$, can be successfully replaced by another set of weights associated with Lagrange polynomials of the same degree and nodal points on the same true curve.

The key point is to determine the transformation matrix, ${\mathbf{T}}_{LB}$, which relates the vector of Bernstein polynomials (B) with that of Lagrange polynomials (L). This relationship depends on the parametric location of the nodal points, thus enabling uniform and non-uniform nodal positions.

Starting with the consituents of the Bernstein-Bézier curve, i.e., the control points ${\mathbf{P}}_B$ and the weights ${\mathbf{w}}_B$, we can easily determine the projected nodal points ${\mathbf{Q}}_L^w$ and their weights ${\mathbf{w}}_L$ in the Lagrangian formulation, and eventually find the nodal points ${\mathbf{Q}}_L$ on the true curve $\mathbf{C}\left(\xi \right)$. Based on them, for any arbitrary parameter $\xi$ we can calculate the associated rational Lagrange polynomials and then apply series expansion to determine the corresponding point $\mathbf{x}\left(\xi \right)$. It was found that the relevant approximation accurately represents a circular arc (maximum error less than ${10}^{-15}$).

The same concept is applicable to a surface (curvilinear patch) as well. As an example, we considered the spherical cap which refers to one-sixth of a sphere. It is well known that this is accurately represented by a rational Bézier surface of degree $p=4$, which is not an ideal tensor product. The latter means that the weights of the interior satisfy the inequality $w_{ij}\ne w_i{w}_j$; therefore, since five control points exist in each direction, the denominator of the normalization will be a sum of 25 terms. A generic expression was presented to construct the total transformation matrix (of size $25\times 25$) in terms of the one-dimensional transformation matrix ${\mathbf{T}}_{LB}$ (of size $5\times 5$) along any edge of the quadrilateral patch. Next, the twenty-five nodal points on the true surface were determined, and for each parametric pair $(\xi ,\eta )$ the spherical cap was accurately represented (maximum deviation less than ${10}^{-15}$).

An additional advantage of using nodal points on the true boundary of the domain, is the fact that stiffness and mass matrices directly refer to the nodal values, $U_i$, and thus boundary conditions are easily imposed. Relevant numerical results–using rational Lagrange polynomials–were presented in Section 6, where excellent convergence was reported. In addition, it is worthy to mention that identical numerical results were found when rational Bernstein-Bézier polynomials were implemented.

8. Conclusions

The usual representation of smooth patches through tensor-product rational Bernstein-Bézier polynomials, can be equivalently substituted by rational Lagrange polynomials. To maintain this numerical equivalency, it is necessary to implement transformation matrices associated with the two mutually perpendicular body-fitted directions, and thus to determine nodal points on the true boundary and corresponding weights. The geometrical concept was applied to a circular arc and a spherical cap, whereas a boundary-value problem was successfully solved for the quarter of an annulus.