定义:令U={u0,u1,…,um}是一个单调不减的实数序列,即ui≤ui+1,i=0,1,…,m-1。其中,ui称为节点,U称为节点矢量,用Ni,p(u)表示第i个p次(p+1阶)B样条基函数,其定义为
由此可知:
(1)Ni,0(u)是一个阶梯函数,它在半开区间u∈[ui,ui+1)外都为零;
(2)当p>0时,Ni,p(u)是两个p-1次基函数的线性组合;
(3)计算一组基函数时需要事先制定节点矢量U和次数p;
(4)定义式中可能出现0/0,我们规定0/0=0;
(5)Ni,p(u)是定义在整个实数轴上的分段多项式函数,但我们一般只对它在区间[u0,um]上的部分感兴趣;
(6)半开区间[ui,ui+1)称为第i个节点区间(knot span),它的长度可以为零,因为相邻节点可以是相同的;
(7)计算p次基函数的生成过程生成一个如下形式的三角形阵列:
为了书写方便,我们通常将Ni,p(u)写为Ni,p。
性质:
(1)(局部支撑性)如果u∉[ui,ui+p+1),则Ni,p(u)=0。
(2)在任意给定的节点区间[uj,uj+1)内,最多p+1个Ni,p是非零的,它们是Nj-p,p,…,Nj,p。
(3)(非负性)对于所有的i,p和u,有Ni,p(u)≥0。
(4)(规范性)对于任意的节点区间[ui,ui+1),当u∈[ui,ui+1)时
(5)(可微性)在节点区间内部,Ni,p(u)是无限次可微的。
(6)除p=0的情况外,Ni,p(u)严格地达到最大值一次。
#include<iostream>
int FindSpan(int n,int p,double u,double* U)
/*计算参数u所在区间的下标返回参数u所在节点区间的下标,即u所在的区间[ui,ui+1]的iU=[u0,u1,u2,...,um]n:节点数组最大下标-1,n=m-1p:次数u:参数U:节点数组*/
{if(u==U[n+1]){return n;}int low=p;int high=n+1;int mid=(low+high)/2;while(u<U[mid] or u>=U[mid+1]){if(u<U[mid]){high=mid;}else{low=mid;}mid=(low+high)/2;}return mid;
}
void BasisFuns(int i,double u,int p,double* U,double* N)
/*计算非零B样条基函数的值i:参数u所在节点区间的下标,即u所在的区间[ui,ui+1]的ip:次数u:参数U:节点数组N:B样条基函数值数组N(i,i-p),...,N(i,p)*/
{double temp=0.0;double saved=0.0;double left[p+1];double right[p+1];N[0]=1.0;for(int j=1;j<=p;j++){left[j]=u-U[i+1-j];right[j]=U[i+j]-u;saved=0.0;for(int r=0;r<j;r++){temp=N[r]/(right[r+1]+left[j-r]);N[r]=saved+right[r+1]*temp;saved=left[j-r]*temp;}N[j]=saved;}
}
int main()
{int n=10;const int p=2;double u=5.0/2;double U[]={0,0,0,1,2,3,4,4,5,5,5};double N[p+1];BasisFuns(4,u,p,U,N);std::cout<<FindSpan(n,p,u,U)<<std::endl;for(int i=0;i<=p;i++){std::cout<<N[i]<<std::endl;}system("pause");
}