The convergence rate and necessary-and-sufficient condition for the consistency of isogeometric collocation method

Although the isogeometric collocation (IGA-C) method has been successfully utilized in practical applications due to its simplicity and efficiency, only a little theoretical results have been established on the numerical analysis of the IGA-C method. In this paper, we deduce the convergence rate of the consistency of the IGA-C method. Moreover, based on the formula of the convergence rate, the necessary and sufficient condition for the consistency of the IGA-C method is developed. These results advance the numerical analysis of the IGA-C method.


§1 Introduction
For the integration of CAD and CAE, Hughes et. al. [22] developed the isogeometric analysis (IGA) method. Since it is based on non-linear NURBS basis functions, the IGA method can directly process the CAD models represented by NURBS, and avoid the tedious mesh transformation procedure.
Because the degree of the non-linear NURBS basis function is relatively high, it is possible to seek a numerical solution, i.e., a NURBS function, by applying the collocation method on the strong form of a differential equation. In this way, the isogeometric collocation (IGA-C) method was proposed [2]. Then unknown coefficients of the NURBS function can be determined by solving a linear system of equations, which is constructed by holding the strong form of the differential equation at some discrete points, called collocation points.
The IGA-C method is a simple and efficient method for solving the unknown coefficients of the NURBS function. A comprehensive study [31] revealed its superior behavior over the Galerkin method in terms of accuracy-to-computational-time ratio. Due to these merits, the IGA-C method has been successfully applied in some practical applications. However, the thorough numerical analysis for the IGA-C method is far from being established. Auricchio et. al. developed numerical analysis of the IGA-C method in one-dimensional case [2]. In the generic case, only some sufficient conditions were presented for the consistency and convergence of the IGA-C method [27].
In this paper, we first develop the convergence rate of the consistency of the IGA-C method, and then present the necessary and sufficient condition for the consistency of the IGA-C method, specifically, for a given boundary (or initial) problem with DT = f (refer to Eq. (1)), where D is its differential operator. Suppose T r is its numerical solution, represented by a NURBS function, and I is an interpolation operator such that If = DT r . The IGA-C method is consistency, if and only if D and I are both uniformly bounded when the knot grid size tends to 0.
The structure of this paper is as follows. In Section 1.1, some related work is briefly reviewed. After introducing some preliminaries in Section 2, an introductory example is presented in Section 3. Moreover, the convergence rate of the consistency of the IGA-C method is deduced in Section 4, and the necessary and sufficient condition is developed in Section 5. In addition, some numerical examples are presented in Section 6. Finally, Section 7 concludes the paper.

Related work
As stated above, the IGA method [22] was proposed to advance the seamless integration of CAD and CAE, by avoiding mesh transformation. Moreover, since it has much less freedom than the traditional finite element method, the IGA method can not only save lots of computation, but also greatly improve the computational precision. Additionally, owing to the knot insertion property of the NURBS function, the shape of the computational domain represented by NURBS can be exactly kept in the mesh refinement. Due to these merits, the IGA method draws great interests in both practical applications and theoretical studies. On the one hand, the IGA method has been successfully applied in lots of simulation problems, such as elasticity [1,17], structure [13,23,33], and fluid [8,9,11], etc. On the other hand, some research on the computational aspect of the IGA method has been developed to improve the accuracy and efficiency by using reparameterization and refinement, etc. [5,7,12,24,36,37]. Recently, an optimal and totally robust multi-iterative method was developed for solving IgA Galerkin linear system [15]. For more details on the IGA method, please refer to Ref. [18] and the references therein.
Since a NURBS function has a relatively high degree, its unknown coefficients can be determined by making the strong form of the PDE hold at some collocation points, which leads to the IGA-C method [2]. Schillinger et. al. presented a comprehensive comparison between the IGA-C method and the Galerkin method, revealing that the IGA-C method is superior to the Galerkin method in terms of accuracy-to-computational-time ratio [31]. Lin et. al. developed some sufficient conditions for the consistency and convergence of the IGA-C method [27]. Moreover, Lorenzis et. al. proposed the IGA-C method for solving the boundary problem with Neumann boundary condition [16].
Besides the traditional collocation schemes, such as Greville points and Demko points [2], some new collocation schemes were developed recently. Based on the orthogonal collocation and superconvergence theory, the superconvergent point scheme was proposed, which converges in the first derivative norms at rates similar to that of the Galerkin solution [6]. Furthermore, a subset of superconvergent points is chosen to reach optimal convergence rate for odd degrees [29]. In Ref. [21], it is shown that the collocation sites that produce the Galerkin solution exactly exist, and the variational collocation scheme is designed. Moreover, Cauchy-Galerkin point scheme was presented so that collocation performed at these points can reproduce the Galerkin solution of various boundary value problems exactly [20]. Recently, Wang et al. developed the superconvergent isogeometric collocation scheme with Greville points [34].
The IGA-C method has been successfully applied in some practical applications. For instance, the IGA-C method was employed in solving Timoshenko beam problem [10] and spatial Timoshenko rod problem [4], showing that mixed collocation schemes are locking-free independently of the choice of the polynomial degrees for unknown fields. Moreover, the IGA-C method was extended to multi-patch NURBS configurations, various boundary and patch interface conditions, and explicit dynamic analysis [3]. Recently, the IGA-C method was exploited to settle the Bernoulli-Euler beam problem [30], the Reissner-Mindlin plate problem [25], Kirchhoff-Love plates and shell problem [28], Reissner-Mindlin shell problem [26], Cosserat rods and rod structure problem [35], and structural dynamics [19]. However, only very limited theoretical results for the IGA-C method were developed [2,27] currently, and the numerical analysis for the IGA-C method is still far from being established. §2 Preliminaries Suppose the IGA-C method is employed to solve the following boundary problem, where Ω ⊂ R d is a physical domain of d dimension, D : V → W is a bounded differential operator, where V and W are two Hilbert spaces, GT is a boundary condition, and f : Ω → R, g : ∂Ω → R are two given continuous functions defined on their domains. Suppose the analytical solution T ∈ C m (Ω), where m is larger than or equal to the maximum order of derivatives appearing in the operator D.
In the IGA method, the physical domain Ω is represented by a NURBS mapping, where Ω p is a parameter domain. Replacing the control points of F by unknown control coefficients, we get the representation of the numerical solution to the boundary problem (Eq. (1)), denoted as T r (η), η ∈ Ω p . Meanwhile, by the inverse mapping F −1 , the physical domain Ω can be mapped into the parameter domain Ω p , and then, the numerical solution T r is still defined on the physical domain Ω through the mapping F −1 . Additionally, by the mapping F , the function f can be defined on Ω p , and G on ∂Ω p .
In isogeometric analysis, while the physical domain of the boundary problem (Eq. (1)) is Ω, the computational domain is the parameter domain Ω p (Eq. (2)). Although the operators D and G in Eq. (1) are performed on the variables in the physical domain, the generated formulae will be transformed into the parameter domain Ω p for computation. Therefore, the functions in the function approximation problem in the IGA-C method should be considered to be defined on the parameter domain Ω p .
Definition 1 (Stable operator [32]). Let V, W be Hilbert spaces and D : V → W be a differential operator. If there exists a constant C S > 0 such that represents the domain of D, then the differential operator D is called a stable operator.
[Remark 1:] In this paper, we suppose that the L ∞ norm ∥·∥ L ∞ is equivalent to the norm ∥·∥ W in W and the norm ∥·∥ V in V. In other words, there exists nonnegative constants c v , C v , and c w , C w satisfying, is an unknown NURBS function defined on the knot grid T ρ ∈ R d , d = 1, 2, 3. Specifically, T ρ is a knot sequence in 1D case, a rectangular grid in 2D case, and a hexahedral grid in 3D case, where ρ is the knot grid size defined as the following definition.
where d(x, y) denotes the Euclidean distance between x and y. And we call ρ as the knot grid size of T ρ , which is defined as the maximum of the diameters of the knot intervals of T ρ .
The modulus of continuity [14] of the function T , denoted as ω(T, h), is defined by The modulus of continuity ω(T, h) satisfies the property [14], Definition 4. Let I ρ be an interpolation operator, and I ρ g be a spline interplant of a function g defined on the knot grid T ρ . Suppose P is a spline space composed of the splines with the same knot grid and degree as those of I ρ g. The distance of the function g to P, i.e., dist(g, P), is defined by

) §3 An introductory example
Consider the following one-dimensional boundary problem: is an analytical solution, and The physical domain [a, b] in Eq. (6) is modeled as, where B i,k (t) is a B-spline basis function of order k (degree k −1), defined on the knot sequence, Then, the numerical solution T r (t) to the boundary problem (Eq. (6)) can be generated by replacing the coefficients a (7)) by the unknowns coefficients p i i = 0, 1, · · · , N , i.e., Note that, by the inverse mapping F −1 and, In order for solving the unknown coefficients in Eq. (10) using the IGA-C method, a linear system is generated by sampling N − 1 points τ 1 , τ 2 , · · · , τ N −1 in the interval (0, 1), i.e., (11) When the knot grid size ρ = 1 N of the knot sequence G (Eq. (8)) tends to 0, it follows N → +∞.
However, because f (x(t)) is continuous on the close interval [0, 1], f (x(τ j )) is bounded. Therefore, if the linear system (Eq. (11)) has a solution, there should exist T r (t) so that the control It results in that dTr dx is also bounded when ρ → 0. All of such B-spline functions T r (t) constitute a B-spline subspace, and the first order derivative operator in Eq. (6) should be bounded on the B-spline subspace when ρ → 0. §4 The convergence rate Suppose the NURBS function T r (η), η ∈ Ω p ⊂ R d defined on the knot grid T ρ has n unknown control coefficients p i , i.e., where w i > 0 are known weights, B i (η) are the B-spline basis functions, the weight function W (η) is a known polynomial spline function, and P (η) is a polynomial spline function with n unknown control coefficients p i . Moreover, the subscript i in Eq. (12) is an index vector, According to the IGA-C method, these unknown coefficients p i can be determined by solving the following linear system of equations, where η k (k = 1, 2, · · · , n 1 ) are collocation points inside Ω p , and η l (l = n 1 + 1, · · · , n) are collocation points on ∂Ω p . Note that, throughout this paper, the operators D and G are performed on the variable in the physical domain Ω (Eq. (1)).
[Remark 2:] In this paper, we assume that the coefficient matrix of the above linear system (Eq. (13)) is of full rank and then it has a unique solution. Otherwise, the IGA-C method is invalid.
According to the result developed in Ref. [27], DT r can be represented by whereB i (η) is the result by applying the differential operator D to w i B i (η) W (η) ,W (η) is the power of W (η), andP (η) is a polynomial B-spline function with n unknowns p i . By Ref. [27],P (η) andW (η) both have the same break point sequence and the same knot intervals as T r (η). To determine these unknowns p i inP (η), let DT r (η) interpolate DT (η) = f (η) at n 1 collocation points inside the domain Ω p (refer to Eq. (13)), i.e., Note thatW (η) ̸ = 0 is a known function, Eq. (15) is equivalent tō Similarly, GT r (η) in Eq. (13) can be written as whereB i (η) are the result generated by applying the operator G to w i B i (η) W (η) ,W (η) ̸ = 0 is a known B-spline function, andP (η) is an unknown B-spline function with n unknowns p i . Then the linear equations GT r (η l ) = g(η l ) in Eq. (13) are equivalent tõ where η l ∈ ∂Ω p , l = n 1 + 1, · · · , n.
Therefore, combining Eqs. (16) and (18), the linear system (Eq. (13)) becomes Since the linear system of equations (Eq. (19)) is equivalent to Eq. (13), then the coefficient matrix of Eq. (19) is of full rank, and it also has a unique solution.
[Remark 3:] In Eq. (13), let the functions f and g vary in C m (Ω p ) and C m (∂Ω p ), respectively, and the differential operator D be fixed. In addition, let the weight function W (η) in T r (Eq. (12)) be fixed as well. Then, all the numerical solutions T r (η) (Eq. (12)) generated by the IGA-C method (Eq. (13)) constitute a linear spline space S ρ (Ω p ), where ρ is the knot grid size of T r . It should be pointed out that, all the NURBS functions T r in the linear space S ρ (Ω p ) have the same weight function W (η), the same knot grid with knot grid size ρ and the same degree. In order for ρ → 0, the knot grid of the spline functions in S ρ (Ω p ) is refined by knot insertion, thus resulting in a series of spline spaces. Moreover, because all the numerical solutions T r (η) constitute the linear space S ρ (Ω p ), all of DT r (η) compose a linear spline space S d ρ,e (Ω p ), where e is the continuity order of the splines in S d ρ,e (Ω p ). As aforementioned, DT r has the same break point sequence as that of T r (Eq. (12)), so they have the same knot grid size ρ.
The following Lemma 1 estimates the distance from a continuous function f ∈ C 0 (Ω p ) to the linear space S d ρ,e (Ω p ), i.e., dist(f, S d ρ,e ). In Ref. [14, pp.146], an inequality to estimate the distance is proposed for univariate functions, and the inequality can be extended to our case. Proof: As stated above, the NURBS functions approximating the analytical solution T constitute the linear space S ρ (Ω p ) defined on the knot grid T ρ . We select a special function from the space S ρ (Ω p ), i.e., and construct a spline function (Af )(η) to approximate the function f ∈ C 0 (Ω p ), i.e., where T is the analytical solution of Eq. (1), and Af = DT r ∈ S d ρ,e (Ω p ) (defined in Remark 3). The point sequence {τ i ∈ Ω p } is sampled in such a way that each knot interval of the knot grid T ρ contains at least one point, and τ i is in the non-zero region of B i (η).
Suppose u(η), v(η) ∈ C m (Ω p ). The function u(η) is an arbitrary function in C m (Ω p ), and, Note that |v(η)| is continuous in the close set Ω p , so |v(η)| can take its maximum value in Ω p . Namely, there exists η * ∈ Ω p such that |v(η * )| = max For an arbitrary valueη ∈ Ω p , it holds, together with Eq. (21), we have, ≤ ∥D∥ ω(T, Kρ), where K is an integer related to the degree of T r (η). Because J is the index set satisfying B i (η * ) ̸ = 0, i ∈ J, and the cardinality of J, i.e., the number of spline basis function B i (η) non-zero at a point, is related to the degree of the spline function T r (η), together with that τ i is in the non-zero region of B i (η), i ∈ J, there exists an integer K related to the degree of T r (η) such that d(η * , τ i ) < Kρ. Therefore, the last step in the above formulae holds. By Eq. (4), we get ρ,e (Ω p ), andη ∈ Ω p is an arbitrary value, it can be so chosen that ρ). Then the Lemma is proved.

Furthermore, we have
where ∇T is the gradient of T , and the norm ∥·∥ E is defined as ∥η∥ E = ∥(η 1 , η 2 , · · · , η d )∥ E = √ η 2 1 + η 2 2 + · · · + η 2 d . Proof: Let x, y ∈ Ω p , d(x, y) = ∥x − y∥ E ≤ ρ, and c ∈ (0, 1). According to the mean value theorem, it follows that Then, by the definition of ω(T, ρ) (Eq. (3)), we have Moreover, we denote by I ρ an interpolation operator, which maps a continuous function to a spline function defined on the knot grid T ρ with knot grid size ρ. In other words, for the continuous function f = DT ∈ C 0 (Ω p )(refer to Eqs. (1) and (13)), we have Specifically, given an arbitrary known NURBS function T q (η) ∈ S ρ (Ω p ) expressed as where the weight function W (η) and the weight w i are the same as those in Eq. (12), two functions h(η) and h b (η) can be generated by performing the operators D and G on T q (see Eq. (1)), respectively, i.e., ) .
The following lemma holds.

Lemma 3. I ρ h = h, where h is defined in Eq. (24).
Proof: We construct an unknown NURBS function T x (η) ∈ S ρ (Ω p ) with n unknown control coefficients x i , the same knot grid and degree with T q (Eq. (23)), i.e., where the weight function W (η) and the weight w i are the same as those in Eqs. (12) and (23). The n unknown coefficients x i in T x (η) can be obtained by making DT x and GT x interpolate h(η) and h b (η) at some sampling points, respectively, similar as Eq. (13), i.e., Therefore, I ρ h = DT x . The aforementioned linear system of equations (Eq. (26)) can be rewritten as Obviously, the coefficient matrix of Eq. (27) is the same as that of the linear system (Eq. (13)), and is of full rank, too. Then the linear system of equations (Eq. (27)) has only zero solution, i.e., x i = q i , meaning that Based on Lemma 3, we can show the following lemma. Eq. (1)), and T r ∈ S ρ (Ω p ) (Refer to Eq. (12) and Remark 3) is the NURBS function approximating the analytical solution T . Then, where I ρ is the interpolation operator defined by Eq. (22), and S d ρ,e is defined as in Remark 3.
Proof: Based on Lemma 3, together with T q and h defined in Eqs. (23) and (24), respectively, we have, (explained in the following paragraph) Because T q (η) (Eq. (23)) is an arbitrary NURBS function in the spline space S ρ (Ω p ), the function h(η) = DT q (Eq. (24)) is also an arbitrary NURBS function in the linear space S d ρ,e (Ω p ). So in Eq. (30), the function h can be chosen from S d ρ,e (Ω p ) to make ∥f − h∥ W as small as possible, that is, dist(f, S d ρ,e ).

and Remark 3) is the NURBS function approximating the analytical solution T . Then,
, where I ρ is the interpolation operator defined by Eq. (22), and K is an integer related to the degree of the splines in the spline space S ρ (Ω p ).
Moreover, due to Lemma 2 and 5, the convergence rate of DT r to DT when ρ → 0 is obtained as follows.

Theorem 1. Suppose the analytical solution
Here, DT, T r , I ρ , and K are delineated as in Lemma 5. In addition, if D is a stable operator (Definition 1), we can get the convergence rate of T r to T when ρ → 0. Eq. (1) is a stable differential operator, and T ∈ C 1 (Ω p ). We have,

Corollary 1. Suppose the operator D in
where C S is a positive constant, T r , I ρ , and K are delineated as in Lemma 5.

One dimensional case
In the one dimensional case, the convergence rate can be improved. In this section, suppose the operator D is a linear differential operator with constant coefficients. Lemma 6. [14, pp. 148] Let g ∈ C m (Ω p ) be a univariate function, and S d ρ,e (Ω p ) be defined as in Remark 3. It holds, where γ is a number related to the degree of the splines in S d ρ,e (Ω p ), and g ′ is the first order derivative of g.
Repeatedly using Lemma 6 leads to: (1)) is a univariate function, the linear spline space S d ρ,e (Ω p ) is defined as in Remark 3, the operator D is a linear differential operator with constant coefficients, and ν = min(m, e) ≥ 1. We have, where Γ is a number related to ν and the degree of the splines in S d ρ,e (Ω p ), and T (ν) is the ν th order derivative of T .

Proof:
The proof is divided into two parts. On the one hand, when µ = 1, by Lemma 1 and 2, it follows, where Γ is a number related to the degree of the splines in S d ρ,e (Ω p ). On the other hand, when µ ≥ 2, because f = DT ∈ C m (Ω p ) is a univariate function, and D is a linear differential operator with constant coefficients, we have (DT ) (k) = DT (k) , k = 1, 2, · · · , m. By using Lemma 6 repeatedly, and denoting ν = min(m, e), it follows, where, γ 1 is a number related to the degree of the splines in S d ρ,e (denoted as deg), γ 2 is a number related to the degree of the splines in S d ρ,e−1 , i.e., deg − 1, · · · , and so on; K ν is a number related to the degree of the splines in S d ρ,e−ν+1 , i.e., deg − ν + 1. In conclusion, γ i , i = 1, 2, · · · , ν − 1, and K ν are all related to deg and ν = min(m, e), and then we denote Γ = γ 1 γ 2 · · · γ ν−1 K ν . Moreover, by Lemma 2, we have, , e), and Γ is a number related to ν and the degree of the splines in S d ρ,e (Ω p ).
Based on Lemma 4 and 7, the convergence rate for the consistency of the IGA-C method in the one-dimensional case is deduced. (1)) is a univariate function, the spline space S d ρ,e (Ω p ) is defined as in Remark 3, and the operator D is a linear differential operator with constant coefficients. We have, (m, e), and Γ is a number related to ν and the degree of the splines in S d ρ,e (Ω p ).
Moreover, if the operator D is also a stable operator (Definition 1), it holds: (1)) is a univariate function, the spline space S d ρ,e (Ω) is defined as in Remark 3, and the linear differential operator with constant coefficients D is stable (refer to Definition 1). We have, a positive constant, ν = min(m, e), and Γ is a number related to ν and the degree of the splines in S d ρ,e (Ω p ). §5 The necessary and sufficient condition In this section, we will present the necessary and sufficient condition of the consistency of the IGA-C method. Because DT = f ∈ C 0 (Ω p ) and T is continuous (Eq. (1)), we have ω(T, ρ) → 0, when ρ → 0. Based on Lemma 5, if ∥I ρ ∥ and ∥D∥ are bounded, it follows ∥DT r − DT ∥ W → 0 when ρ → 0. That is, the IGA-C method is consistent. However, since I ρ f = DT r ∈ S d ρ,e , and T r ∈ S ρ is defined on the knot grid T ρ with knot grid size ρ, the norms ∥I ρ ∥ and ∥D∥ are both related to the knot grid size ρ. Therefore, the sufficient condition for the consistency of the IGA-C method is followed.
Furthermore, the following lemma presents the necessary condition for the consistency of the IGA-C method.

Proof:
We employ the method of proof by contradiction to show that DT r is bounded when ρ → 0.
The consistency of the IGA-C method means that DT r → DT = f, when ρ → 0. (32) By contradiction, suppose DT r is not uniformly bounded when ρ → 0, i.e., ∥DT r ∥ W → ∞, when ρ → 0. Because f is continuous, it is bounded on its domain Ω p ∪ ∂Ω p . However, DT r is unbounded when ρ → 0. This violates the consistency condition (Eq. (32)). So the hypothesis is not true, DT r is uniformly bounded when ρ → 0. That is, there exists a positive constant C r such that ∥DT r ∥ W ≤ C r , when ρ → 0. Therefore, we have ∥D∥ = sup and (refer to Eq. (22)) It means that the interpolation operator I ρ (Eq. (22)) and the differential operator D (Eq. (1)) are both uniformly bounded when ρ → 0.
Based on Lemmas 8 and 9, the necessary and sufficient condition for the consistency of the IGA-C method is followed.
∥f ∥ L ∞ is uniformly bounded when the knot grid size sequence ρ k → 0.
In Fig. 1, the diagrams of where the degrees of the numerical solution T r range from 2 to 7. It can be seen from the diagrams in Fig. 1 that, when k → ∞ and ρ k → 0, as an indicator of ∥D∥ L ∞ , the ratio Therefore, it is uniformly bounded as ρ k → 0, (k → ∞), which validates Theorem 3.
Similarly, the diagrams of Fig. 2. It can be seen that the ratio ∥f ∥ L ∞ , as an indicator of ∥I ρ ∥, is also uniformly bounded as ρ k → 0, (k → ∞).
The physical domain Ω in Eq. (34) is a quarter of an annulus, which is represented by a cubic NURBS patch with 4 × 4 control points. The control points and weights of the cubic NURBS patch are listed in Tables 1 and 2, respectively. The knot vectors of the cubic NURBS patch along u− and v−direction are, respectively, 0 0 0 0 1 1 1 1, 0 0 0 0 1 1 1 1. To make the knot grid size tend to 0, we uniformly insert knots in the interval (0, 1) along u− and v−directions, respectively. So, the knot grid sizes are ρ k = 1 k+1 , k = 0, 1, 2, · · · . Fig. 3 shows the diagrams ∥DTr∥ L ∞ ∥Tr∥ L ∞ v.s. lg ρ k for the case of two-dimensional source problem (Eq. (34)), where the degrees of the numerical solution T r range from 2 to 7. Similar as the case of one-dimensional problem, ∥DTr∥ L ∞ ∥Tr∥ L ∞ tend to a limit when ρ k → 0(k → ∞). So they are all uniformly bounded when ρ k → 0 (k → ∞).  In this paper, we developed the convergence order for the consistency and convergence of the IGA-C method, and then, deduced the necessary-and-sufficient condition for the consistency of the IGA-C method. Specifically, suppose D is the differential operator of a boundary value problem with DT = f (Eq. (1)), a NURBS function T r is the numerical solution, and I ρ is an interpolation operator such that I ρ f = DT r . First, the formula of the convergence order for the consistency of the IGA-C method is developed, which includes the norms of the operator D and I ρ . Then, the necessary-and-sufficient condition for the consistency of the IGA-C method is deduced. That is, the IGA-C method is consistent if and only if D and I ρ are both uniformly bounded when ρ → 0. These results will advance the numerical analysis of the IGA-C method.