b.cpp

是计算B样条函数的一个源代码,搞数值计算的时候对样条函数进行计算

#include<fstream>
#include<iostream>
#include<vector>
using namespace std;
void main(){
	
	vector<double> t,x,y,z,b;
	int n,i,j,k;
	vector<vector<double> > B;                           // 定义变量，由vector所包含的文件为所给数据分配内存


	ifstream infile("data.txt", ios::in);                //读入文件"data.txt"
	infile >> n;                                         //读入n的值
	t.resize(2*n-4);                                     //重新分配内存空间
	x.resize(n-4);
	y.resize(n-4);
	z.resize(n-4);
	b.resize(n-4);
	B.resize(n-1);
	for(i=0;i<=n-2;i++)
	B[i].resize(n-4);

	for(i=1;i<=2*n-4;i++)
	{infile >> t[i];                                      //读入结点值与函数值，由于n后有2*n-4个值，其中前n个为结点值，后面的n-4个为样条函数在greville坐标下的函数值
    cout<< t[i] <<endl;}
	for(j=1;j<=n-4;j++)
	y[j]=t[j+n];                                          //读入函数值
	infile.close();


	for(i=1;i<=n-4;i++)
	for(j=i;j<=i+4;j++)
		x[i]+=t[j]/5;                                     //建立Greville坐标


    for(j=1;j<=n-4;j++)                                   //形成插值矩阵B[i][j]
	{
	  for(i=1;i<=n-1;i++)
	  {
		 if((x[j]>=t[i])&&(x[j]<t[i+1])) B[i-1][j]=1;
		 else B[i-1][j]=0;  
	  };
	  for(k=1;k<=3;k++)
		  for(i=1;i<=n-k-1;i++)
		  {   double a,b;
			  if(t[i+k]==t[i]) a=0;
			  else a=(x[j]-t[i])/(t[i+k]-t[i]);
			  if(t[i+k+1]==t[i+1]) b=0;
		      else b=(t[i+k+1]-x[j])/(t[i+k+1]-t[i+1]);
		      B[i-1][j]=B[i-1][j]*a+B[i][j]*b;            //由B样条的递归定义求得高次B样条函数的值
		  }; 
	}


    for(i=1;i<=n-4;i++) z[i]=0;                           //设初始解为(0,0,...,0)
     for(k=1;k<=20;k++)                                  //由高斯-塞得尔迭代解方程组,k为迭代次数
	 {                                                    // Bz=y得到z[i]  
	    for(i=1;i<=n-4;i++)
		{  double p=0;
	        for(j=1;j<=n-4;j++) 
				p+=B[i-1][j]*z[j];
	          p=p-B[i-1][i]*z[i];
                z[i]=(y[i]-p)/B[i-1][i];
		};
	 };

     double xvalue;
     cout<<"Please input xvalue:"<<endl;
     cin >> xvalue;
     for(i=1;i<=n-1;i++)                                  //解出x处的b[i]=B[i](x)
	  {
		  if((xvalue>=t[i])&&(xvalue<t[i+1])) b[i]=1;
		  else b[i]=0;   
	  };
     for(k=1;k<=3;k++)
	  for(i=1;i<=n-k-1;i++)
	  {       double alpha,beta;
			  if(t[i+k]==t[i]) alpha=0;
			  else alpha=(xvalue-t[i])/(t[i+k]-t[i]);
			  if(t[i+k+1]==t[i+1]) beta=0;
		      else beta=(t[i+k+1]-xvalue)/(t[i+k+1]-t[i+1]);
		      b[i]=b[i]*alpha+b[i+1]*beta;
		 }; 


     double S=0;                                          //对z[i]*b[i]求和得到S(x)
     for(i=1;i<=n-4;i++) S+=z[i]*b[i];                     
     
	 
	 cout<<"the value of Bsphine S(x) in point x is:"<< S <<endl;
}