三次样条插值法的C++程序实现（追赶法求解线性方程组）

三次样条插值法的C++程序实现（追赶法求解线性方程组）

求解线性方程组追赶法的程序为：

double capway(int n,double nu[],double d[],double lan[],double h[],double x[],double y[], double t){
    double beita[n-1];
    beita[0]=lan[0]/2;
    for (int i=1;i<n-1;i++){
    	beita[i]=lan[i]/(2-nu[i-1]*beita[i-1]);
	}
    double yy[n];
    yy[0]=d[0]/2;
    for (int i=1;i<n;i++){
    yy[i]=(d[i]-nu[i-1]*yy[i-1])/(2-nu[i-1]*beita[i-1]);	
	}
	double xx[n];
	xx[n-1]=yy[n-1];
	for (int i=n-2;i>=0;i--){
		xx[i]=yy[i]-beita[i]*xx[i+1];
	}
	//上面求解线性方程组成功。
	}

计算线性方程组矩阵的程序为：

double spline(int n,double x[],double y[],int s,double x1,double xn,double t){
	//下面求hn的值；
   double h[n-1];
   for (int i=0;i<n-1;i++){
   	h[i]=x[i+1]-x[i];
   } 
	double aver[n-1];
	for (int i=0;i<n-1;i++){
	aver[i]=(y[i+1]-y[i])/h[i];		
	}
	//下面判断边界条件,第一类边界条件； 
    double lan[n-1];
    double d[n];
    double nu[n-1];
    if(s==1){
    	lan[0]=1;
    	d[0]=(6/h[0])*(aver[0]-x1);
    	nu[n-2]=1;
    	d[n-1]=(6/h[n-2])*(xn-aver[n-2]);
    	cout<<d[n-1]<<"\n";
	}
	if(s==2){
		lan[0]=0;
    	d[0]=2*x1;
    	nu[n-2]=0;
    	d[n-1]=2*xn;
	}
	//下面求nu；
	for (int i=0;i<n-2;i++){
		nu[i]=h[i]/(h[i]+h[i+1]);
	}
	for (int i=1;i<n-1;i++){
		lan[i]=h[i]/(h[i-1]+h[i]);
	}
	for (int i=1;i<n-1;i++){
		d[i]=6*(aver[i]-aver[i-1])/(h[i-1]+h[i]);
	}
	 double a[n][n];
	 for (int i=0;i<n;i++){
	 	a[i][i]=2;
	 }
	 for (int i=0;i<n-1;i++){
	 	a[i+1][i]=nu[i];
	 }
	 for (int i=0;i<n-1;i++){
	 	a[i][i+1]=lan[i];
	 }	
	 cout<<"输出线性方程组矩阵"<<"\n"; 
	 for (int i=0;i<n;i++){
	 	for (int j=0;j<n;j++){
	 	cout<<a[i][j]<<"     "; 	
		 }
		 cout<<"\n";
	 }
	 cout<<"输出线性方程组系数向量\n";
	 for (int i=0;i<n;i++){
	 	cout<<d[i]<<"  ";
	 }
	 cout<<"\n";
}

总的程序为

#include<iostream>
#include<cstdio> 
using namespace std;
//capway为用追赶发求解线性方程组； 
double capway(int n,double nu[],double d[],double lan[],double h[],double x[],double y[], double t){
    double beita[n-1];
    beita[0]=lan[0]/2;
    for (int i=1;i<n-1;i++){
    	beita[i]=lan[i]/(2-nu[i-1]*beita[i-1]);
	}
    double yy[n];
    yy[0]=d[0]/2;
    for (int i=1;i<n;i++){
    yy[i]=(d[i]-nu[i-1]*yy[i-1])/(2-nu[i-1]*beita[i-1]);	
	}
	double xx[n];
	xx[n-1]=yy[n-1];
	for (int i=n-2;i>=0;i--){
		xx[i]=yy[i]-beita[i]*xx[i+1];
	}
	//上面求解线性方程组成功。
	int as;
	for (int i=0;i<n;i++){
		if(t>=x[i]){
			as=i;
		}
	} 
	//下面计算函数值；
	if(as!=n-1){
	double numble;
	double sum1;
	double sum2;
	double sum3;
	double sum4;
	sum1=xx[as]*(x[as+1]-t)*(x[as+1]-t)*(x[as+1]-t)/6/h[as];
	sum2=xx[as+1]*(t-x[as])*(t-x[as])*(t-x[as])/6/h[as];
	sum3=(y[as]-xx[as]*h[as]*h[as]/6)*(x[as+1]-t)/h[as];
	sum4=(y[as+1]-xx[as+1]*h[as]*h[as]/6)*(t-x[as])/h[as];
	numble=sum1+sum2+sum3+sum4;
    cout<<t<<"的三次样条插值为"; 
	printf("%f",numble);
	}
}
double spline(int n,double x[],double y[],int s,double x1,double xn,double t){
	//下面求hn的值；
   double h[n-1];
   for (int i=0;i<n-1;i++){
   	h[i]=x[i+1]-x[i];
   } 
	double aver[n-1];
	for (int i=0;i<n-1;i++){
	aver[i]=(y[i+1]-y[i])/h[i];		
	}
	//下面判断边界条件,第一类边界条件； 
    double lan[n-1];
    double d[n];
    double nu[n-1];
    if(s==1){
    	lan[0]=1;
    	d[0]=(6/h[0])*(aver[0]-x1);
    	nu[n-2]=1;
    	d[n-1]=(6/h[n-2])*(xn-aver[n-2]);
    	cout<<d[n-1]<<"\n";
	}
	if(s==2){
		lan[0]=0;
    	d[0]=2*x1;
    	nu[n-2]=0;
    	d[n-1]=2*xn;
	}
	//下面求nu；
	for (int i=0;i<n-2;i++){
		nu[i]=h[i]/(h[i]+h[i+1]);
	}
	for (int i=1;i<n-1;i++){
		lan[i]=h[i]/(h[i-1]+h[i]);
	}
	for (int i=1;i<n-1;i++){
		d[i]=6*(aver[i]-aver[i-1])/(h[i-1]+h[i]);
	}
	 double a[n][n];
	 for (int i=0;i<n;i++){
	 	a[i][i]=2;
	 }
	 for (int i=0;i<n-1;i++){
	 	a[i+1][i]=nu[i];
	 }
	 for (int i=0;i<n-1;i++){
	 	a[i][i+1]=lan[i];
	 }	
	 cout<<"输出线性方程组矩阵"<<"\n"; 
	 for (int i=0;i<n;i++){
	 	for (int j=0;j<n;j++){
	 	cout<<a[i][j]<<"     "; 	
		 }
		 cout<<"\n";
	 }
	 cout<<"输出线性方程组系数向量\n";
	 for (int i=0;i<n;i++){
	 	cout<<d[i]<<"  ";
	 }
	 cout<<"\n";
    capway(n,nu,d,lan,h,x,y,t);
}
int main(){
    int num;
    cout<<"输入节点的个数："; 
    scanf("%d",&num);
    double x[num],y[num];
    for (int i=0;i<num;i++){
    	cout<<"输入第"<<i+1<<"个x的值："; 
	    scanf("%lf",&x[i]);
		cout<<"输入第"<<i+1<<"个y的值：";
		scanf("%lf",&y[i]);
	}
	printf("输入边界条件，如果输入1，表示第一类边界条件，2表示第二类边界条件");
	int s;
	cin>>s; 
	double x1,xn; 
	printf("输入边界条件值\n");
	cout<<"输入左端边界条件值:";
	cin>>x1;
	cout<<"输入右端边界条件值:";
	cin>>xn; 
	printf("输入要计算的函数值:"); 
	double t;
	cin>>t;
	spline(num,x,y,s,x1,xn,t);
	return 0;  
}

运行结果如下：