IG-LSPIA: Least Squares Progressive Iterative Approximation for Isogeometric Collocation Method

Abstract

The isogeometric collocation method (IGA-C), which is a promising branch of isogeometric analysis (IGA), can be considered fitting the load function with the combination of the numerical solution and its derivatives. In this study, we develop an iterative method, isogeometric least-squares progressive-iterative approximation (IG-LSPIA), to solve the fitting problem in the collocation method. IG-LSPIA starts with an initial blending function, where the control coefficients are combined with the B-spline basis functions and their derivatives. A new blending function is generated by constructing the differences for collocation points (DCP) and control coefficients (DCC), and then adding the DCC to the corresponding control coefficients. The procedure is performed iteratively until the stop criterion is reached. We prove the convergence of IG-LSPIA and show that the computation complexity in each iteration of IG-LSPIA is related only to the number of collocation points and unrelated to the number of control coefficients. Moreover, an incremental algorithm is designed; it alternates with knot refinement until the desired precision is achieved. After each knot refinement, the result of the last round of IG-LSPIA iterations is used to generate the initial blending function of the new round of iteration, thereby saving great computation. Experiments show that the proposed method is stable and efficient. In the three-dimensional case, the total computation time is saved twice compared to the traditional method.

1. Introduction

The isogeometric analysis, proposed by Hughes et al. [1], is an important research field with a wide range of applications in Computer-aided Engineering and Computer-aided Design (CAE/CAD) due to its shared representation for CAD and CAE. The isogeometric collocation method (IGA-C) is an efficient numerical method for isogeometric analysis [2,3,4], where the numerical solution is obtained by solving a linear system formed by substituting collocation points into the strong form of a partial differential equation (PDE). In conventional IGA-C, the number of collocation points and the degrees of freedom (DoFs) of the numerical solution are the same. In the isogeometric least-squares collocation method (IGA-L), the number of collocation points is larger than the DoFs of the numerical solution [5]. However, the main drawback of those methods is that the exact solution of a PDE is usually unknown, and the approximation error of the numerical solution to the exact solution cannot be calculated. Thus, determining how many degrees of freedom are required for a numerical solution is difficult. Actually, IGA-L can be considered fitting the load function (the right-hand function of a PDE) with the combination of a numerical solution and its derivatives. In this study, we develop the least squares progressive iterative approximation method for solving the linear systems in the isogeometric collocation method (IG-LSPIA), and then show its convergence. The computational complexity in each iteration of IG-LSPIA is related only to the number of collocation points and unrelated to the DoFs of the numerical solution; thus, it is suitable for large-scale data fitting and incremental data fitting. Then, we propose an incremental approach for IG-LSPIA, namely, incremental IG-LSPIA (Incre-IG-LSPIA), where the freedom of the numerical solution can be automatically determined using the fitting precision.

As an iterative method with geometric meaning, the progressive–iterative approximation for least square fitting (LSPIA) [6,7] has been widely used in CAD [8], where a sequence of curves (surfaces) is generated iteratively to fit a given data set. The limit of the sequence is proven to be the least-squares solution of the fitting problem. Moreover, Lin et al. [9] proved that the LSPIA method is convergent whether or not the iterative matrix is singular. However, LSPIA has never been used in solving PDE. In this study, IG-LSPIA is developed for solving IGA-L. In IG-LSPIA, a sequence of blending functions is constructed to approximate the load function, and the limit of the sequence is the solution of IGA-L [5]. Specifically, during each iteration step of IG-LSPIA, we initially calculate the differences for collocation points (DCP) and then the differences for control coefficients (DCC) of the numerical solution. Subsequently, we can update the control coefficients by adding the DCC. This procedure is performed iteratively until the stop criteria are reached. Moreover, given a fitting precision, IG-LSPIA can be performed incrementally, and the approach is called Incre-IG-LSPIA. That is, if the fitting error of the current round of IG-LSPIA does not satisfy the given fitting precision, then we further raise the degree of freedom of the numerical solution using knot refinement and invoke a new round of IG-LSPIA iterations. In the new round of iterations, the result of the previous round iterations after knot refinement can be considered the initial value for the new round of iterations, saving greatly on computations. The Incre-IG-LSPIA stops when the given fitting precision of the numerical solution is reached. In this approach, the DoFs of the numerical solution can be automatically determined. The contributions of this study are summarized as follows.

• We propose the IG-LSPIA, where the convergence rate is related to the number of collocation points.
• An incremental fitting algorithm is proposed with the IG-LSPIA method, which saves computational costs.
• Using the incremental IG-LSPIA and fitting load function, the DoFs of the numerical solution can be determined automatically.

The paper is organized as follows: in Section 2, we provide a brief description of LSPIA and isogeometric collocation methods. The background of IGA-C with non-uniform rational B-splines (NURBS) is presented in Section 3. In Section 4, we present the IG-LSPIA method. The convergence analysis is demonstrated in Section 5, and we show the numerical experiments in Section 6. The conclusion and some ideas about future work are presented in Section 7.

2. Related Work

In this section, we provide a concise review of the LSPIA and isogeometric collocation method.

2.1. LSPIA

Motivated by the work of Qi et al. [10] and de Boor [11], Lin et al. [6] proposed the progressive-iterative approximation (PIA) method for fitting data points using the NURBS curve and surface. As a data-fitting algorithm without solving the linear system, its main idea is to generate a sequence of curves by iteratively adjusting the control points. Thus, the limit of the curve sequence can approach and finally interpolate the given data points. Based on the local property of the PIA, Lin [12] provided a local data fitting method for curves and surfaces. This method only updates some control points to realize the interpolation of the corresponding data points. Subsequently, the studies have mainly focused on two aspects: one is to accelerate the PIA method, and the other is to generalize the algorithm to more forms of curves and surfaces. Lu [13] and Zhang et al. [14] used the weighted difference vectors to accelerate the convergence of the PIA method. However, the former multiplies all difference vectors by the same optimal weight, whereas the latter uses the weight matrix. Chen et al. [15] set up a progressive interpolation for Catmull–Clark subdivision surfaces. Cheng et al. [16] further developed a PIA for loop subdivision surface fitting. Chen and Wang [17] designed a PIA format for triangular Bézier.

Although these techniques are intuitive and efficient, they require that the number of data points be equal to that of the control points, and they cannot solve the least-squares problem. In particular, Lin and Zhang [18] developed an extended PIA format (EPIA), and Deng and Lin [7] subsequently proposed the LSPIA to deal with such problems. In those methods, the number of data points could be much larger than that of control points, which is a common situation in practical fitting problems. Meanwhile, an incremental fitting idea is proposed, where each new round of iteration uses the previous fitting result as the initial guess. As a result, the convergence rate is enhanced. Lin et al. [9] further verified that even if a singular matrix appeared in the problem, LSPIA still converged under certain conditions. This result theoretically guarantees the applicability of the LSPIA algorithm to wider fitting problems. Zhang et al. [19,20] then applied LSPIA to generalized B-spline curves and NURBS. In addition to data fitting, PIA has been widely used in path planning [21], data compression [22], and implicit surface generation [23,24]. For more details, we refer the reader to the review [8].

2.2. Isogeometric Collocation Methods

The theoretical study of the isogeometric collocation method has always been a challenge. Auricchio et al. [2] first developed a theoretical analysis of the algorithm in the 1D case. Lin et al. [3] then explored the algorithm’s consistency and convergence. Furthermore, the algorithm’s convergence order and the necessary and sufficient condition to ensure the algorithm’s consistency are provided by Lin et al. [4].

Another major bottleneck in the usage of IGA-C is that the linear systems derived from it can be typically ill-conditioned when the grid size and degree of spline increase. Remarkable progress has been made in recent years to improve the stability and convergence rate of IGA-C. It focuses mainly on selecting appropriate collocation points or methods for solving linear equations.

For example, Anitescu et al. [25] applied collocation points derived from the superconvergence theory to obtain an optimal convergence order of IGA-C. Donatelli et al. [26] introduced a fast multigrid solver based on the spectral analysis of the collocation matrix. The symbol functions capturing the spectral information could be beneficial for designing efficient preconditioners and solvers, but they are challenging to construct. In addition, Veiga et al. [27] constructed the overlapping additive Schwarz (OAS) preconditioners, resulting in a robust formulation. Cho [28] later developed optimal multilevel preconditioners for isogeometric collocation methods and showed their consistent performance under various conditions, including basis function, spline degree, and domain deformations. Lin et al. [5] proposed a least-squares collocation method that allows a more flexible number of collocation points and is shown to be more stable than the conventional IGA-C. In this work, we adopt an incremental, iterative method to reuse the information in the last round of iterations and thus shorten the total computational cost.

3. Preliminaries

In this section, we present a short description of the basics of the isogeometric collocation method and IGA-L method [5].

The following boundary value problem is considered:

$$\begin{array}{c}\left\{\begin{array}{cc}\mathcal{L}u=f,\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \mathcal{G}u=g,& \mathrm{on}\partial \mathsf{\Omega },\hfill \end{array}\right.\hfill \end{array}$$

(1)

where $u:\mathsf{\Omega }\to \mathbb{R}$ is the unknown solution, $\mathcal{L}$ is a differential operator, $\mathcal{G}$ is a boundary operator, and f and g are the continuous functions. The main idea of the isogeometric analysis algorithm is to adopt the same basis functions to represent the numerical solution $u_h$ and the geometric map G. In industrial production and applications, we usually select the NURBS basis functions, i.e.,

$$R_j\left(\xi \right)=\frac{N_{j,p}\left(\xi \right)w_j}{{\displaystyle \sum _{\widehat{j}=1}^n}{N}_{\widehat{j},p}\left(\xi \right)w_{\widehat{j}}},$$

where $N_{j,p}\left(\xi \right)$ denotes the B-spline basis function of degree p, and $w_j$ is the weight. Multidimensional NURBS basis functions can be defined analogically using tensor product

$$\begin{array}{cc}& R_{ij}(\xi ,\eta )=\frac{N_{i,p_1}\left(\xi \right)N_{j,p_2}\left(\eta \right)w_{ij}}{{\displaystyle \sum _{\widehat{i}=1}^{n_1}}{\displaystyle \sum _{\widehat{j}=1}^{n_2}}{N}_{\widehat{i},p_1}\left(\xi \right)N_{\widehat{j},p_2}\left(\eta \right)w_{\widehat{i}\widehat{j}}},\hfill \\ & R_{ijk}(\xi ,\eta ,\zeta )=\frac{N_{i,p_1}\left(\xi \right)N_{j,p_2}\left(\eta \right)N_{k,p_3}\left(\zeta \right)w_{ijk}}{{\displaystyle \sum _{\widehat{i}=1}^{n_1}}{\displaystyle \sum _{\widehat{j}=1}^{n_2}}{\displaystyle \sum _{\widehat{k}=1}^{n_3}}{N}_{\widehat{i},p_1}\left(\xi \right)N_{\widehat{j},p_2}\left(\eta \right)N_{\widehat{k},p_3}\left(\zeta \right)w_{\widehat{i}\widehat{j}\widehat{k}}},\hfill \end{array}$$

where $p_1$, $p_2$, and $p_3$ indicate the degree of the basis functions in the relevant parameters’ directions.

In the isogeometric collocation method, a finite set of collocation points are further introduced to discretize the partial differential Equation (1) directly and avoid constructing the variational form. The selection of collocation points is not unique but highly important. Greville abscissae [29] is adopted in general, but Demko abscissae [30] and superconvergent points [25] exist. Without loss of generality, we consider the 1D case in the following. Let ${\widehat{\tau }}_i,i=1,2,\dots ,m$, be the collocation points on the parametric domain ${\mathsf{\Omega }}_p$. Then, the corresponding collocation points ${\tau }_i\in \mathsf{\Omega }$ can be calculated using

$${\tau }_i=G\left({\widehat{\tau }}_i\right)=\sum _{j=1}^n{R}_j\left({\widehat{\tau }}_i\right)P_j,i=1,2,\dots ,m.$$

(2)

Furthermore, the numerical solution is defined as

$$u_h\left(\widehat{\tau }\right)=\sum _{j=1}^n{R}_j\left(\widehat{\tau }\right)u_j,$$

(3)

or

$$u_h\left(\tau \right)=\sum _{j=1}^n{R}_j\left(G^{-1}\left(\tau \right)\right)u_j,$$

where $u_j,j=1,2,\dots ,n$ is the control coefficient. The parametric domain is set as ${\mathsf{\Omega }}_p=[0,1]$ and ${\widehat{\tau }}_1=0,{\widehat{\tau }}_m=1$, and the index set of the collocation points $\mathcal{I}$ can be partitioned into two subsets, as follows:

$$\begin{array}{cc}& {\mathcal{I}}_{\mathcal{L}}=\{2,3,\dots ,m-1\},\hfill \\ & {\mathcal{I}}_{\mathcal{G}}=\{1,m\},\hfill \end{array}$$

where ${\mathcal{I}}_{\mathcal{L}}$ and ${\mathcal{I}}_{\mathcal{G}}$ denote the indices of the interior and boundary collocation points, respectively. Substituting Equation (2) into Equation (1), the numerical solution $u_h\left(\widehat{\tau }\right)$ satisfies the rules

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill \mathcal{L}{u}_h\left({\widehat{\tau }}_i\right)=f\left(G\left({\widehat{\tau }}_i\right)\right),& i\in {\mathcal{I}}_{\mathcal{L}},\hfill \\ \hfill \mathcal{G}{u}_h\left({\widehat{\tau }}_i\right)=g\left(G\left({\widehat{\tau }}_i\right)\right),& i\in {\mathcal{I}}_{\mathcal{G}},\hfill \end{array}\right.\end{array}$$

(4)

Arranging the control coefficients $u_j$ of the numerical solution $u_h\left(\widehat{\tau }\right)$ into an n-dimensional column vector, i.e., $\mathit{U}={[u_1,u_2,\dots ,u_n]}^T$, the discretized problem (4) can be reformulated in matrix form using

$$\mathit{AU}=\mathit{b},$$

(5)

with

$$\begin{array}{cc}{A}_{ij}\hfill & =A_j\left({\widehat{\tau }}_i\right)=\left\{\begin{array}{cc}\hfill \mathcal{L}{R}_j\left({\widehat{\tau }}_i\right),& i\in {\mathcal{I}}_{\mathcal{L}},\hfill \\ \mathcal{G}{R}_j\left({\widehat{\tau }}_i\right),\hfill & i\in {\mathcal{I}}_{\mathcal{G}},j=1,2,\dots ,n,\hfill \end{array}\right.\hfill \\ b_i\hfill & =b\left({\widehat{\tau }}_i\right)=\left\{\begin{array}{cc}f\left(G\left({\widehat{\tau }}_i\right)\right),\hfill & i\in {\mathcal{I}}_{\mathcal{L}},\hfill \\ g\left(G\left({\widehat{\tau }}_i\right)\right),\hfill & i\in {\mathcal{I}}_{\mathcal{G}},\hfill \end{array}\right.\hfill \end{array}$$

(6)

where $\mathit{A}$ is the collocation matrix, and $\mathit{b}$ denotes the load vector. Finally, solving the system of Equation (5), we obtain the control coefficients ${\left\{u_j\right\}}_{j=1}^n$ and further achieve the numerical solution $u_h\left(\widehat{\tau }\right)$ by substituting the control coefficients into Equation (3).

The amount of collocation points is assumed to be more than that of the control coefficients ${\left\{u_j\right\}}_{j=1}^n$ and geometric control points ${\left\{P_j\right\}}_{j=1}^n$, i.e., $m>n$. Then, the linear system of Equation (5) can be solved using the IGA-L method, that is, to find a least-squares solution satisfying

$$\underset{\mathit{U}}{min}{∥\mathit{AU}-\mathit{b}∥}^2,$$

(7)

where $\mathit{A}$ and $\mathit{b}$ are the collocation matrix and load vector, respectively.

4. IG-LSPIA Method

In this section, we present the least-squares progressive-iterative approximation for isogeometric collocation methods (IG-LSPIA). We mainly describe the 1D case, and its generalization to the multidimensional cases is straightforward.

The discretized values ${\left\{b_i\right\}}_{i=1}^m$ of load function in Equation (6) is considered an ordered point set. The collocation points on the parametric domain $\left\{0={\widehat{\tau }}_1<{\widehat{\tau }}_2<\dots <{\widehat{\tau }}_m=1\right\}$ are their corresponding parameters. Given the initial guess of the control coefficients ${\left\{u_j^{\left(0\right)}\right\}}_{j=1}^n$ with $n<m$, we obtain an initial blending function

$$U^{\left(0\right)}\left(\widehat{\tau }\right)=\sum _{j=1}^n{A}_j\left(\widehat{\tau }\right)u_j^{\left(0\right)},\widehat{\tau }\in [{\widehat{\tau }}_1,{\widehat{\tau }}_m],$$

(8)

where $A_j\left(\widehat{\tau }\right)$, $j=1,2,\dots ,n$, represents the combination of different order derivatives of the NURBS basis functions determined using the operators $\mathcal{L}$ and $\mathcal{G}$ (1), i.e.,

$$A_j\left(\widehat{\tau }\right)=\left\{\begin{array}{c}\mathcal{L}{R}_j\left(\widehat{\tau }\right),\widehat{\tau }\in {\mathsf{\Omega }}_p^{in},\hfill \\ \mathcal{G}{R}_j\left(\widehat{\tau }\right),\widehat{\tau }\in {\mathsf{\Omega }}_p^{bd}j=1,2,\dots ,n.\hfill \end{array}\right.$$

${\mathsf{\Omega }}_p^{in}$ and ${\mathsf{\Omega }}_p^{bd}$ indicate the interior and boundary of the parameter domain, respectively. Each $A_j\left(\widehat{\tau }\right)$ corresponds to the jth control coefficient. The initial control coefficients can be selected randomly. In our experiments, we choose $u_j^{\left(0\right)}=0$, $j=1,2,\dots ,n$. Because when the collocation matrix is singular, if the initial control coefficients are zero, the IG-LSPIA method converges to the M-P pseudo-inverse solution of the linear system with the minimum Euclidean norm [9].

Assume that $J_{in}$ and $J_{bd}$ are the index sets of the internal and boundary control coefficients, respectively. Without loss of generality, we further assume that the boundary control coefficients have been obtained using strong or weak imposition [31] and are fixed, i.e.,

$$u_j^{\left(k\right)}=u_j^{*},j\in J_{bd},k=0,1,2,\dots .$$

(9)

Initially, we calculate the differences for collocation points (DCP) in the first step of iteration

$$\begin{array}{cc}\hfill {\delta }_i^{\left(0\right)}& =b_i-\sum _{j=1}^n{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(0\right)}\hfill \\ & =b_i-\sum _{j\in J_{bd}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(0\right)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(0\right)},i=1,2,\dots ,m,\hfill \end{array}$$

and group all load values whose parameters fall in the local support of the jth derivatives function, i.e., $A_j\left(\widehat{\tau }\right)\ne 0$, into the jth group corresponding to the jth control coefficient. In the following, we denote the subscript of the jth group of load values using $I_j$.

Subsequently, we derive the differences for control coefficients (DCC) as follows:

$$d_j^{\left(0\right)}=\mu \sum _{h\in I_j}{A}_j\left({\widehat{\tau }}_h\right){\delta }_h^{\left(0\right)},j=1,2,\dots ,n,$$

where $\mu$ is a normalization technique to guarantee the convergence of our algorithm. Note that since the assumption (9), the boundary DCCs are zero, i.e., $d_{j\in J_{bd}}^{\left(k\right)}=0,k=0,1,2,\dots$. Then, the new control coefficients $u_j^{\left(1\right)}$ are generated by adding $d_j^{\left(0\right)}$ to $u_j^{\left(0\right)}$, i.e.,

$$u_j^{\left(1\right)}=u_j^{\left(0\right)}+d_j^{\left(0\right)},j=1,2,\dots ,n.$$

Substituting them into Equation (8), we finally obtain the new blending function, as follows:

$$U^{\left(1\right)}\left(\widehat{\tau }\right)=\sum _{j=1}^n{A}_j\left(\widehat{\tau }\right)u_j^{\left(1\right)},\widehat{\tau }\in [{\widehat{\tau }}_1,{\widehat{\tau }}_m].$$

The kth blending function, generated after the kth iteration of IG-LSPIA, is assumed to be as follows:

$$U^{\left(k\right)}\left(\widehat{\tau }\right)=\sum _{j=1}^n{A}_j\left(\widehat{\tau }\right)u_j^{\left(k\right)},\widehat{\tau }\in [{\widehat{\tau }}_1,{\widehat{\tau }}_m].$$

Then, the DCP in the $(k+1)$st iteration are obtained using

$$\begin{array}{cc}\hfill {\delta }_i^{\left(k\right)}& =b_i-\sum _{j=1}^n{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(k\right)}\hfill \\ & =b_i-\sum _{j\in J_{bd}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(k\right)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(k\right)},i=1,2,\dots ,m,\hfill \end{array}$$

(10)

Furthermore, the DCC are obtained as follows:

$$d_j^{\left(k\right)}=\mu \sum _{h\in I_j}{A}_j\left({\widehat{\tau }}_h\right){\delta }_h^{\left(k\right)},j=1,2,\dots ,n.$$

(11)

Thus, the new control coefficients are updated via the following relation:

$$u_j^{(k+1)}=u_j^{\left(k\right)}+d_j^{\left(k\right)},j=1,2,\dots ,n.$$

(12)

Consequently, the $(k+1)$st blending function is as follows:

$$U^{(k+1)}\left(\widehat{\tau }\right)=\sum _{j=1}^n{A}_j\left(\widehat{\tau }\right)u_j^{(k+1)}.$$

The above iteration process is performed until the desired fitting error is reached and a sequence of blending functions is obtained.

$$\left\{U^{\left(k\right)}\left(\widehat{\tau }\right),k=0,1,\dots \right\}.$$

(13)

5. Convergence Analysis of the Iterative Scheme

In this section, we show that the IG-LSPIA is convergent and converges to the solution of a constrained optimization problem derived from the least-squares collocation method [5].

Assume that $I_{in}$ is the index set of the internal collocation points.

Based on Equations (9), (10), and (12), the $(k+1)$st internal DCP follow the rule:

$$\begin{array}{cc}\hfill {\delta }_i^{(k+1)}& =b_i-\sum _{j=1}^n{A}_j\left({\widehat{\tau }}_i\right)u_j^{(k+1)}\hfill \\ & =b_i-\sum _{j\in J_{bd}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{(k+1)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{(k+1)}\hfill \\ & =b_i-\sum _{j\in J_{bd}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{(k+1)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)\left[u_j^{\left(k\right)}+d_j^{\left(k\right)}\right]\hfill \\ & =b_i-\sum _{j\in J_{bd}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(k\right)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(k\right)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)d_j^{\left(k\right)}\hfill \\ & =b_i-\sum _{j=1}^n{A}_j\left({\widehat{\tau }}_i\right)u_j^{\left(k\right)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)d_j^{\left(k\right)}\hfill \\ & ={\delta }_i^{\left(k\right)}-\sum _{j\in J_{in}}{A}_j\left({\widehat{\tau }}_i\right)d_j^{\left(k\right)},i\in I_{in}.\hfill \end{array}$$

(14)

According to Equations (11) and (14), the $(k+1)$st internal DCC are derived using

$$\begin{array}{cc}\hfill d_j^{(k+1)}& =\mu \sum _{h\in I_j}{A}_j\left({\widehat{\tau }}_h\right){\delta }_h^{(k+1)}\hfill \\ & =\mu \sum _{h\in I_j\cap I_{in}}{A}_j\left({\widehat{\tau }}_h\right){\delta }_h^{(k+1)}\hfill \\ & =\mu \sum _{h\in I_j\cap I_{in}}{A}_j\left({\widehat{\tau }}_h\right)\left[{\delta }_h^{\left(k\right)}-\sum _{t\in J_{in}}{A}_t\left({\widehat{\tau }}_h\right)d_t^{\left(k\right)}\right]\hfill \\ & =\mu \sum _{h\in I_j\cap I_{in}}{A}_j\left({\widehat{\tau }}_h\right){\delta }_h^{\left(k\right)}-\mu \sum _{h\in I_j\cap I_{in}}{A}_j\left({\widehat{\tau }}_h\right)\left[\sum _{t\in J_{in}}{A}_t\left({\widehat{\tau }}_h\right)d_t^{\left(k\right)}\right]\hfill \\ & =d_j^{\left(k\right)}-\mu \sum _{h\in I_j\cap I_{in}}{A}_j\left({\widehat{\tau }}_h\right)\left[\sum _{t\in J_{in}}{A}_t\left({\widehat{\tau }}_h\right)d_t^{\left(k\right)}\right],j\in J_{in}.\hfill \end{array}$$

(15)

To show the convergence of the IG-LSPIA, Equation (15) can be rewritten in matrix form as follows:

$$\begin{array}{cc}\hfill {\mathit{D}}_{in}^{(k+1)}& ={\mathit{D}}_{in}^{\left(k\right)}-\mu {\mathit{B}}^T\mathit{B}{\mathit{D}}_{in}^{\left(k\right)}\hfill \\ & =\left({\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}\right){\mathit{D}}_{in}^{\left(k\right)},\hfill \end{array}$$

(16)

with

$$\begin{array}{cc}\hfill {\mathit{D}}_{in}^{\left(k\right)}& ={\left[d_{n_1+1}^{\left(k\right)},d_{n_1+2}^{\left(k\right)},\dots ,d_n^{\left(k\right)}\right]}^T,\hfill \\ \hfill \mathit{B}& =\left[\begin{array}{cccc}{A}_{n_1+1}\left({\widehat{\tau }}_{m_1+1}\right)& A_{n_1+2}\left({\widehat{\tau }}_{m_1+1}\right)& \dots & A_n\left({\widehat{\tau }}_{m_1+1}\right)\\ A_{n_1+1}\left({\widehat{\tau }}_{m_1+2}\right)& A_{n_1+2}\left({\widehat{\tau }}_{m_1+2}\right)& \dots & A_n\left({\widehat{\tau }}_{m_1+2}\right)\\ ⋮& ⋮& ⋮& ⋮\\ A_{n_1+1}\left({\widehat{\tau }}_m\right)& A_{n_1+2}\left({\widehat{\tau }}_m\right)& \dots & A_n\left({\widehat{\tau }}_m\right)\end{array}\right],\hfill \end{array}$$

(17)

where $n_1$ and $m_1$ ($m_1>n_1$) are the numbers of the boundary control coefficients and collocation points, respectively. $n_1+1,n_1+2,\dots ,n$ are the serial numbers of internal control coefficients, ${\widehat{\tau }}_{m_1+1},$${\widehat{\tau }}_{m_1+2},\dots ,{\widehat{\tau }}_m$ are the internal collocation points. Denote $n_2$ and $m_2$ ($m_2>n_2$) are the numbers of the internal control coefficients and collocation points, respectively. We have $m=m_1+m_2$ and $n=n_1+n_2$. In addition, ${\mathit{I}}_2$ is the identity matrix with $\mathrm{rank}\left({\mathit{I}}_2\right)=n_2$, and $\mathit{B}$ is an $m_2\times n_2$ matrix. From Equation (16), the iterative matrix for internal control coefficients is ${\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}$.

Lemma 1.

If the normalization weight satisfies $0<\mu <\frac{2}{{\lambda }_{max}\left({\mathit{B}}^T\mathit{B}\right)}$, where $\mathit{B}$ is given in Equation (17), then the eigenvalues of $\mu {\mathit{B}}^T\mathit{B}$ are real numbers and satisfy $0\le \lambda \left(\mu {\mathit{B}}^T\mathit{B}\right)<2$.

Proof. 

The ${\mathit{B}}^T\mathit{B}$ matrix is generally a positive semi-definite matrix; thus, the eigenvalues are nonnegative real numbers as well as $\mu {\mathit{B}}^T\mathit{B}$, i.e.,

$$\lambda \left(\mu {\mathit{B}}^T\mathit{B}\right)\ge 0.$$

Moreover, we select $0<\mu <\frac{2}{{\lambda }_{max}\left({\mathit{B}}^T\mathit{B}\right)}$, and then we obtain the following:

$$\begin{array}{cc}\hfill \lambda \left(\mu {\mathit{B}}^T\mathit{B}\right)& \le {\lambda }_{max}\left(\mu {\mathit{B}}^T\mathit{B}\right)\hfill \\ & \le \frac{2}{{\lambda }_{max}\left({\mathit{B}}^T\mathit{B}\right)}{\lambda }_{max}\left({\mathit{B}}^T\mathit{B}\right)\hfill \\ & <2.\hfill \end{array}$$

In summary, the eigenvalues of the matrix $\mu {\mathit{B}}^T\mathit{B}$ satisfy $0\le \lambda <2$, thereby completing the proof.    □

Remark 1.

In practice, we select $\mu =\frac{2}{{∥{\mathit{B}}^T\mathit{B}∥}_{\infty }}$, which satisfies the condition in Lemma 1.

Theorem 1.

When ${\mathit{B}}^T\mathit{B}$ is nonsingular, where $\mathit{B}$ is defined in Equation (17), the iterative format (16) is convergent.

Proof. 

From Equation (16), we obtain the following equations:

$$\begin{array}{cc}\hfill {\mathit{D}}_{in}^{(k+1)}& =\left({\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}\right){\mathit{D}}_{in}^{\left(k\right)}\hfill \\ & ={\left({\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}\right)}^2{\mathit{D}}_{in}^{(k-1)}\hfill \\ & =\dots \hfill \\ & ={\left({\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}\right)}^k{\mathit{D}}_{in}^{\left(0\right)}.\hfill \end{array}$$

(18)

According to Lemma 1 and the assumptions, the eigenvalues of matrix $\mu {\mathit{B}}^T\mathit{B}$ satisfy $0<\lambda \left(\mu {\mathit{B}}^T\mathit{B}\right)<2$, resulting in the following:

$$-1<1-\lambda \left(\mu {\mathit{B}}^T\mathit{B}\right)<1,$$

Thus, the spectral radius satisfies $\rho \left({\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}\right)<1$, leading to

$$\underset{k\to \infty }{lim}{\left({\mathit{I}}_2-\mu {\mathit{B}}^T\mathit{B}\right)}^k=\mathbf{0}.$$

(19)

Substituting Equation (19) into Equation (18), we have ${lim}_{k\to \infty }{\mathit{D}}_{in}^{\left(k\right)}=\mathbf{0}$, i.e., the iterative format (16) is convergent.    □

Furthermore, by substituting Equation (10) into Equation (11) and reformulating it into the matrix form, we have

$${\mathit{D}}_{in}^{\left(k\right)}=\mu {\mathit{B}}^T\left({\mathit{b}}_{in}-{\mathit{A}}_{21}{\mathit{U}}_{bd}^{*}-\mathit{B}{\mathit{U}}_{in}^{\left(k\right)}\right),$$

(20)

where

$$\begin{array}{cc}\hfill {\mathit{U}}_{bd}^{*}& ={\left[u_1^{*},u_2^{*},\dots ,u_{n_1}^{*}\right]}^T,{\mathit{U}}_{in}^{\left(k\right)}={\left[u_{n_1+1}^{\left(k\right)},u_{n_1+2}^{\left(k\right)},\dots ,u_n^{\left(k\right)}\right]}^T,\hfill \\ \hfill {\mathit{b}}_{in}& ={\left[b_{m_1+1},b_{m_1+2},\dots ,b_m\right]}^T,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\mathit{A}}_{21}& =\left[\begin{array}{cccc}{A}_1\left({\widehat{\tau }}_{m_1+1}\right)& A_2\left({\widehat{\tau }}_{m_1+1}\right)& \dots & A_{n_1}\left({\widehat{\tau }}_{m_1+1}\right)\\ A_1\left({\widehat{\tau }}_{m_1+2}\right)& A_2\left({\widehat{\tau }}_{m_1+2}\right)& \dots & A_{n_1}\left({\widehat{\tau }}_{m_1+2}\right)\\ ⋮& ⋮& \ddots & ⋮\\ A_1\left({\widehat{\tau }}_{m_1+m_2}\right)& A_2\left({\widehat{\tau }}_{m_1+m_2}\right)& \dots & A_{n_1}\left({\widehat{\tau }}_{m_1+m_2}\right)\end{array}\right].\hfill \end{array}$$

(22)

Theorem 2.

When ${\mathit{B}}^T\mathit{B}$ is nonsingular, where $\mathit{B}$ is defined in Equation (17), the IG-LSPIA format (16) converges to the solution of the constrained least-squares collocation problem, as follows:

$$\begin{array}{c}\underset{\mathit{U}}{min}{∥\mathit{AU}-\mathit{b}∥}^2,\\ s.t\mathit{H}\mathit{U}={\mathit{U}}_{bd}^{*},\end{array}$$

(23)

where $\mathit{H}$ is an $n_1\times n$ matrix

$$\mathit{H}=\left[\begin{array}{cc}{\mathit{I}}_{\mathbf{1}}& \mathbf{0}\end{array}\right].$$

$\mathit{A}$ and $\mathit{b}$ are defined similarly to Equation (6). ${\mathit{U}}_{bd}^{*}$ consists of the known control coefficients on the boundary, and ${\mathit{I}}_1$ is the identity matrix with $\mathit{rank}\left({\mathit{I}}_1\right)=n_1$.

Proof. 

Based on Theorem 1, we have ${lim}_{k\to \infty }{\mathit{D}}_{in}^{\left(k\right)}=\mathbf{0}$. Then, when $k\to \infty$, Equation (20) leads to

$$\mu {\mathit{B}}^T\mathit{B}{\mathit{U}}_{in}^{(\infty )}=\mu {\mathit{B}}^T\left({\mathit{b}}_{in}-{\mathit{A}}_{21}{\mathit{U}}_{bd}^{*}\right).$$

Since $\mu$ is non-zero, the internal control coefficients converge to the solution of

$${\mathit{B}}^{\mathit{T}}\mathit{B}{\mathit{U}}_{in}={\mathit{B}}^T\left({\mathit{b}}_{in}-{\mathit{A}}_{21}{\mathit{U}}_{bd}^{*}\right)$$

under the IG-LSPIA algorithm.

On the other hand, we can partition the collocation matrix $\mathit{A}$ of problem (23) into

$$\mathit{A}=\left[\begin{array}{cc}{\mathit{A}}_{11}& {\mathit{A}}_{12}\\ {\mathit{A}}_{21}& \mathit{B}\end{array}\right],$$

where

$$\begin{array}{cc}\hfill {\mathit{A}}_{11}& =\left[\begin{array}{cccc}{A}_1\left({\tau }_1\right)& A_2\left({\tau }_1\right)& \dots & A_{n_1}\left({\tau }_1\right)\\ A_1\left({\tau }_2\right)& A_2\left({\tau }_2\right)& \dots & A_{n_1}\left({\tau }_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ A_1\left({\tau }_{m_1}\right)& A_2\left({\tau }_{m_1}\right)& \dots & A_{n_1}\left({\tau }_{m_1}\right)\end{array}\right],\hfill \\ \hfill {\mathit{A}}_{12}& =\left[\begin{array}{cccc}{A}_{n_1+1}\left({\tau }_1\right)& A_{n_1+2}\left({\tau }_1\right)& \dots & A_n\left({\tau }_1\right)\\ A_{n_1+1}\left({\tau }_2\right)& A_{n_1+2}\left({\tau }_2\right)& \dots & A_n\left({\tau }_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ A_{n_1+1}\left({\tau }_{m_1}\right)& A_{n_1+2}\left({\tau }_{m_1}\right)& \dots & A_n\left({\tau }_{m_1}\right)\end{array}\right],\hfill \end{array}$$

${\mathit{A}}_{\mathbf{21}}$ (22) and $\mathit{B}$ (17) are the same as the previous definitions. Subsequently, denote

$$\mathit{U}=\left[\begin{array}{c}{\mathit{U}}_{bd}\\ {\mathit{U}}_{in}\end{array}\right],\mathit{b}=\left[\begin{array}{c}{\mathit{b}}_{bd}\\ {\mathit{b}}_{in}\end{array}\right],$$

with

$$\begin{array}{cc}\hfill {\mathit{U}}_{bd}& ={\left[u_1,u_2,\dots ,u_{n_1}\right]}^T,{\mathit{U}}_{in}={\left[u_{n_1+1},u_{n_1+2},\dots ,u_n\right]}^T,\hfill \\ \hfill {\mathit{b}}_{bd}& ={\left[b_1,b_2,\dots ,b_{m_1}\right]}^T,\hfill \end{array}$$

and ${\mathit{b}}_{in}$ are load values defined in Equation (21). Then, we obtain the following equation:

$$\mathit{AU}-\mathit{b}=\left[\begin{array}{cc}{\mathit{A}}_{11}& {\mathit{A}}_{12}\\ {\mathit{A}}_{21}& \mathit{B}\end{array}\right]\left[\begin{array}{c}{\mathit{U}}_{bd}\\ {\mathit{U}}_{in}\end{array}\right]-\left[\begin{array}{c}{\mathit{b}}_{bd}\\ {\mathit{b}}_{in}\end{array}\right].$$

According to the boundary condition and because ${\mathit{A}}_{12}=0$, we have ${\mathit{A}}_{11}{\mathit{U}}_{bd}^{*}-{\mathit{b}}_{bd}=\mathbf{0}$. Then, the constrained least-squares collocation problem (23) becomes

$$\underset{{\mathit{U}}_{in}}{min}{∥{\mathit{BU}}_{in}-\left({\mathit{b}}_{in}-{\mathit{A}}_{21}{\mathit{U}}_{bd}^{*}\right)∥}^2,$$

(24)

which is equivalent to the normal equations

$${\mathit{B}}^T\mathit{B}{\mathit{U}}_{in}={\mathit{B}}^T\left({\mathit{b}}_{in}-{\mathit{A}}_{21}{\mathit{U}}_{bd}^{*}\right).$$

In conclusion, the IG-LSPIA method converges to the solution of a constrained optimization problem (23) derived from the least-squares collocation method.    □

Remark 2.

After obtaining the solution of the constrained least-squares collocation problem (23) using IG-LSPIA, i.e., ${\mathit{U}}_{in}^{(\infty )}={min}_{{\mathit{U}}_{in}}{∥{\mathit{BU}}_{in}-\left({\mathit{b}}_{in}-{\mathit{A}}_{21}{\mathit{U}}_{bd}^{*}\right)∥}^2$ and ${\mathit{U}}_{bd}={\mathit{U}}_{bd}^{*}$, we can derive the numerical solution of the partial differential Equation (1) by substituting the control coefficients

$${\mathit{U}}^{(\infty )}={\left[u_1^{(\infty )},u_2^{(\infty )},\dots ,u_{n_1}^{(\infty )},u_{n_1+1}^{(\infty )},u_{n_1+2}^{(\infty )},\dots ,u_n^{(\infty )}\right]}^T$$

with $u_j^{(\infty )}=u_j^{*}$, $j=1,2,\dots ,n_1$, into Equation (3), namely,

$$u_h\left(\widehat{\tau }\right)=\sum _{j=1}^n{R}_j\left(\widehat{\tau }\right)u_j^{(\infty )}.$$

(25)

6. Numerical Examples

In this section, we test several numerical examples to demonstrate the convergence properties and applications of the IG-LSPIA method. Particularly, we consider 1D, 2D, and 3D partial differential equations with homogenous Dirichlet boundaries. The boundary conditions are imposed directly on the control coefficients [1]. In our experiments, we use the NURBS functions to build the geometric map G (2) and the fitting error of the kth iteration is denoted with

$$E_{fit}^{\left(k\right)}=\sqrt{\frac{{\sum }_{i=1}^m{\left({\delta }_i^{\left(k\right)}\right)}^2}{m}},k=0,1,2,\dots ,$$

(26)

where ${\delta }_i^{\left(k\right)}$ (10) indicates the DCP after the kth iteration. In particular, $E_{fit}^{\left(0\right)}$ indicates the initial fitting error with a zero initial guess of control coefficients. The stopping criterion is ${∥{\mathit{U}}^{\left(k\right)}-{\mathit{U}}^{(k-1)}∥}_2<{10}^{-7}$.

When the iterations converge, but the current fitting error $E_{fit}$ does not reach the predetermined precision $ϵ$, we adopt the knot refinement strategy and then proceed with the incremental iterations. During the incremental iteration, each new round of iteration after knot refinement reuses information from the last round of iteration to avoid rebuilding the system. In summary, the Incre-IG-LSPIA is shown in Algorithm 1.

Algorithm 1 The Incremental IG-LSPIA Algorithm

Input: Initial control coefficients ${\mathit{U}}^{\left(0\right)}={\left\{u_j\right\}}_{j=1}^n$,    Discretized load functions $\mathit{b}={\left\{b_i\right\}}_{i=1}^m$,    and desired fitting accuracy $\epsilon >0$. 1: Conduct IG-LSPIA with the initial control coefficients ${\mathit{U}}^{\left(0\right)}$; 2: while $E_{fit}^{\left(k\right)}>\epsilon$do 3:   Generate the new control coefficients, control points, and corresponding weights using knot refinement; 4:   Use the new control coefficients as the initial control coefficients for the next iteration; 5:   while ${∥{\mathit{U}}^{\left(k\right)}-{\mathit{U}}^{(k-1)}∥}_2>{10}^{-7}$ do 6:     Conduct IG-LSPIA to update control coefficients; 7:     $k=k+1$; 8:   end while 9:end while Output: Numerical solution $u_h\left(x\right)$.

When the fitting error reaches the desired accuracy, we can evaluate the error between the analytical solution $u\left(\widehat{\tau }\right)$ and the current numerical solution $u_h\left(\widehat{\tau }\right)$ (25) using the uniform sampling points. Here three forms of error measures are employed, i.e., relative $L_2$ norm error, relative $L_{\infty }$ norm error, and root-mean-square error (RMSE)

$$\begin{array}{cc}\hfill E_p& =\frac{{∥\mathit{u}-{\mathit{u}}_{\mathit{h}}∥}_p}{{∥\mathit{u}∥}_p},p=2,\infty ,\hfill \\ \hfill E_R& =\sqrt{\frac{\sum _{i=1}^m{\left(u\left(\widehat{{\tau }_i}\right)-u_h\left(\widehat{{\tau }_i}\right)\right)}^2}{m}.}\hfill \end{array}$$

We implement the algorithms in C++ and run them on a PC with an Intel Core i7-10700 2.90 GHz CPU and 16 GB RAM.

6.1. 1D Tests

The first case is defined in domain $[0,1]$ with boundary conditions, as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-0.001u^{″}+10u^{\prime }+5000u=\left(5000+0.004{\pi }^2\right)sin\left(2\pi x\right)+20\pi cos\left(2\pi x\right),x\in \left(0,1\right),\hfill \\ u\left(0\right)=u\left(1\right)=0,\hfill \end{array}\right.\end{array}$$

(27)

where the analytical solution is $u\left(x\right)=sin\left(2\pi x\right)$. We initially select 17 initial control points and 600 Greville collocation points. The maximum number of iterations per round is set to 500. The geometric map G (2) is based on the cubic NURBS basis function.

In Figure 1, we provide the convergence order of the IG-LSPIA algorithm. The horizontal axis is the logarithm (base 10) of the degrees of freedom (DoFs), and the vertical axis is the logarithm (base 10) of the error. We evaluate the convergence rates under three different errors, i.e., relative $L_2$ norm error, relative $L_{\infty }$ norm error, and RMSE, denoted with $E_1$, $E_{\infty }$, and $E_R$, respectively. The decreasing trend of the three error types is consistent, indicating that the IG-LSPIA algorithm is convergent and stable.

The analytical solution to the problem (27) is illustrated in Figure 2a. We perform a knot refinement after each round of iteration to achieve the prescribed fitting accuracy, and in this example, we refine the knots five times. Figure 2b demonstrates the absolute error distribution curves of the different numbers of control coefficients obtained using knot refinement. The absolute error range gradually decreases as the count of refinements increases.

In addition, we apply the incremental algorithm to speed up the iteration after each knot refinement.

Table 1. Experimental data of Incre-IG-LSPIA algorithm with five times knot refinement. The number of collocation points is $m=600$.

# Ctrl 1	${\mathit{E}}_{\mathit{fit}}^{\left(\mathit{k}\right)}$	# Iter 2	Iteration Time (s)	Iteration Time per Step (s)	Total Time (s)
17	0.144	140	0.380	0.003	0.480
31	0.009	62	0.166	0.003	0.199
59	0.001	55	0.147	0.003	0.180
115	$6.168\times {10}^{-5}$	73	0.203	0.003	0.240
227	$1.780\times {10}^{-5}$	4	0.012	0.003	0.045
451	$5.380\times {10}^{-6}$	1	0.003	0.003	0.040

¹ Number of control coefficients. ² Number of iterations.

shows the experimental data for Incre-IG-LSPIA, i.e., the fitting error $E_{fit}^{\left(k\right)}$ (26), the number of iterations, the iteration time in seconds, the iteration time per step, and the total computation time per round iteration to reach the same stopping criterion. For Incre-IG-LSPIA, the total computation time mainly includes the computation time of the convergence coefficients $\mu$ shown in Remark 1 and the iteration time. In particular, the Incre-IG-LSPIA algorithm does not need to reassemble the collocation matrix explicitly, but it requires computing the convergence coefficients. Moreover, the average time to reassemble the collocation matrix is usually much higher than the average time to compute the convergence coefficients, approximately three times in this 1D case, i.e., 0.124 s and 0.042 s.

6.2. 2D Tests

The second case is as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& -△u+u=f,(x,y)\in \mathsf{\Omega },\hfill \\ & {u|}_{\partial \mathsf{\Omega }}=0,\hfill \end{array}\right.\end{array}$$

(28)

where

$$\begin{array}{cc}\hfill f& =(3x^4-67x^2-67y^2+3y^4+6x^2{y}^2+116)sin\left(x\right)sin\left(y\right)\\ & +(68x-8x^3-8x{y}^2)cos\left(x\right)sin\left(y\right)\\ & +(68y-8y^3-8y{x}^2)cos\left(y\right)sin\left(x\right),\end{array}$$

and $\mathsf{\Omega }$ is a quarter of an annulus domain. The analytical solution for this problem is $u(x,y)=(x^2+y^2-1)(x^2+y^2-16)sin\left(x\right)sin\left(y\right)$. We initially select $7\times 7$ initial control points and $31\times 31$ Greville collocation points. The maximum number of iterations per round is set to 5000. In 2D and 3D, we use the square of the 2-norm error between the steps of iteration as the convergence error. Moreover, we reduce the error tolerance by 0.1 after each round of iteration. As a result, each round iteration is sufficient. The diagram of $lo{g}_{10}\left(n^{1/2}\right)$ vs. $lo{g}_{10}\left(error\right)$ is illustrated in Figure 3. The figure shows that the IG-LSPIA algorithm has consistent convergence under the three errors $E_2$, $E_{\infty }$, and $E_R$.

The analytical solution to the problem (28) is shown in Figure 4a. The absolute error distribution of the initial round of iteration using $7\times 7$ control points is plotted in Figure 4b. Moreover, the results after the first and second knot refinement, i.e., $11\times 11$ and $19\times 19$ control coefficients are shown in Figur Figure 4c,d. More experimental details are demonstrated in

Table 2. Experimental data of Incre-IG-LSPIA algorithm with twice knot refinement. The number of collocation points is $m=31\times 31$.

# Ctrl 1	${\mathit{E}}_{\mathit{fit}}^{\left(\mathit{k}\right)}$	# Iter 2	Iteration Time (s)	Iteration Time per Step (s)	Total Time (s)
$7\times 7$	8.926	365	0.352	0.001	0.371
$11\times 11$	1.898	2506	2.591	0.001	2.605
$19\times 19$	0.394	5000	5.132	0.001	5.147

¹ Number of control coefficients. ² Number of iterations.

. The maximal absolute error in Figure 4 decreases gradually from 1.200 to 0.052, showing that the numerical solution is closer to the analytical solution when the fitting error of the load function is smaller. In addition, the iteration time per step of each iteration of the algorithm does not fluctuate significantly with the increase in control points, that is at 0.001 s.

6.3. 3D Tests

The third case is as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-△u+u=f,(x,y,z)\in \mathsf{\Omega },\hfill \\ {u|}_{\partial \mathsf{\Omega }}=0,\hfill \end{array}\right.\end{array}$$

(29)

with $f=(1+12{\pi }^2)sin\left(2\pi x\right)sin\left(2\pi y\right)sin\left(2\pi z\right)$. The analytical solution is $u(x,y,z)=sin\left(2\pi x\right)sin\left(2\pi y\right)sin\left(2\pi z\right)$. The 3D physical domain $\mathsf{\Omega }$ is modeled as a B-spline solid. We select $30\times 30\times 30$ collocation points and $7\times 7\times 7$ initial control coefficients.

The diagram of $lo{g}_{10}\left(n^{1/3}\right)$ v.s. $lo{g}_{10}\left(error\right)$ is illustrated in Figure 5. The RMSE error has the smallest value, but the overall trend is the same for all three errors. The error decreases rapidly after the first knot refinement and slows down the rate of increase after the second knot refinement because the number of control points gradually approaches the number of collocation points. The analytical solution to the problem (29) is plotted in Figure 6a. Furthermore, the absolute error distributions after three round iterations are shown in Figure 6b–d. The maximal absolute error is decreased from 0.032 to $4.743\times {10}^{-4}$ when the number of control coefficients arises from $7\times 7\times 7$ to $19\times 19\times 19$.

Table 3. Experimental data of Incre-IG-LSPIA algorithm with twice knot refinement. The number of collocation points is $m=30\times 30\times 30$.

# Ctrl 1	${\mathit{E}}_{\mathit{fit}}^{\left(\mathit{k}\right)}$	# Iter 2	Iteration Time (s)	Iteration Time per Step (s)	Total Time (s)
$7\times 7\times 7$	2.979	69	10.000	0.145	11.460
$11\times 11\times 11$	0.562	69	10.993	0.159	12.113
$19\times 19\times 19$	0.119	147	23.623	0.161	24.733

¹ Number of control coefficients. ² Number of iterations.

shows the data from the 3D experiments. The table shows that the average iteration time per step is approximately $0.159$ s. In particular, in the last round iteration, the conventional method requires $18.997$ s to reassemble the collocation matrix and $23.571$ s to solve the equation directly, for a total time of $42.668$ s. However, the Incre-IG-LSPIA method does not need to compute the collocation matrix explicitly but requires computing the convergence coefficient, which consumes $1.11$ s. In the end, the Incre-IG-LSPIA takes $23.623$ s to iterate and a total time of $24.733$ s, much less than solving the equations directly.

7. Discussion

In the one, two, and three-dimensional experiments, the IG-LSPIA has a consistent convergence rate under three different error metrics, i.e., relative $L_2$ norm error, relative $L_{\infty }$ norm error, and RMSE. In addition, the DoFs of the numerical solutions in IG-LSPIA are automatically determined using the fitting error rather than the approximate difference between the numerical and exact solutions. Since the exact solution is usually unknown, the proposed method is more practical than the traditional approach for engineering applications. Regarding time consumption, the cost of reassembling the collocation matrix in traditional algorithms is high, but IG-LSPIA does not require the explicit calculation of the matrix. Although the proposed algorithm involves calculating the normalization weight, the best choice for this parameter has been discussed. It only needs to be computed in the first step of each iteration round. In the one-dimensional example, this time is three times less than the time for reassembly. In the three-dimensional example, the total computation times of the last round in IG-LSPIA is $24.733$ s, while rebuilding the equations and directly solving them cost $42.668$ s. The efficiency has been increased twice. Finally, the Incre-IG-LSPIA method is adopted in our experiments to increase the convergence rate. Our method has limitations. Despite the convergence of the proposed method, the number of iterations is higher in the two-dimensional example than in the other cases.

8. Conclusions

In this study, we present the IG-LSPIA method, a novel and stable iterative format for IGA-C. Moreover, the convergence analysis of the IG-LSPIA method is provided. We discuss sufficient conditions for the iterative algorithm’s convergence and the optimal choice of the normalization weight $\mu$. To enhance iterative efficiency, we develop the Incre-IG-LSPIA to reuse the information in the last round of iterations and save iteration time to some extent. Finally, many examples are tested to demonstrate the efficiency and effectiveness of the developed IG-LSPIA method. In the three-dimensional case, twice the total computation time is saved compared to the traditional method.

For future work, an acceleration of the iterative process should be developed. The proposed method also shows the potential to solve more complicated PDEs. In addition to IGA-C, IG-LSPIA can be combined with the isogeometric Galerkin method.