三次样条拟和.cpp

根据给定点和区间值得出三次样条拟合方程的系数

#include <stdio.h>

 

#define MAX_N 20

 

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, 

double a[], double b[], double c[], double d[]);

 

double S( double t, double Fmax, 

            int n, double x[], double a[], double b[], double c[], double d[] );

 

int main()

{
	freopen("input.txt","r",stdin);
	freopen("output.txt","w",stdout);


    int n, Type, m, i;

    double x[MAX_N], f[MAX_N], a[MAX_N], b[MAX_N], c[MAX_N], d[MAX_N];

    double s0, sn, Fmax, t0, tm, h, t;

 

    while (scanf("%d", &n) != EOF) {

       for (i=0; i<=n; i++) 

           scanf("%lf", &x[i]);

       for (i=0; i<=n; i++) 

           scanf("%lf", &f[i]);

       scanf("%d %lf %lf %lf", &Type, &s0, &sn, &Fmax);

 

       Cubic_Spline(n, x, f, Type, s0, sn, a, b, c, d);

       for (i=1; i<=n; i++)

           printf("%12.8e %12.8e %12.8e %12.8e \n", a[i], b[i], c[i], d[i]);

 

       scanf("%lf %lf %d", &t0, &tm, &m);

       h = (tm-t0)/(double)m;

       for (i=0; i<=m; i++) {

           t = t0+h*(double)i;

           printf("f(%12.8e) = %12.8e\n", t, S(t, Fmax, n, x, a, b, c, d));

       }

       printf("\n");

    }

 

    return 0;

}

 
/* Your functions will be put here */

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[])
{
	double h[MAX_N],e[MAX_N],u[MAX_N],z[MAX_N],L[MAX_N];
	for (int i=0;i<=n;i++)
		a[i]=f[i];
	for (int i=0;i<n;i++)
		h[i]=x[i+1]-x[i];
	
	if (Type==1)				  // clamped 
	{
		e[0]=(a[1]-a[0])*3.0/h[0]-3*s0;
		e[n]=3*sn-(a[n]-a[n-1])*3.0/h[n-1];
	
	}
	for	(int i=1;i<n;i++ )
		e[i]=3.0*(a[i+1]-a[i])/h[i] -3.0*(a[i]-a[i-1])/h[i-1];
						
	
	if(Type==2)                //natural
	{
		L[0]=1;
		u[0]=z[0]=0;
	}
	else                            // clamped
	{
		L[0]=2*h[0];
		u[0]=0.5;
		z[0]=e[0]/L[0];
	}
	
	for (int i=1; i<n;i++)
	{
		L[i]=2.0*(x[i+1]-x[i-1])-h[i-1]*u[i-1];
		u[i]=h[i]*1.0/L[i];
		z[i]=(e[i]-h[i-1]*z[i-1])*1.0/L[i];
	}
		
	if(Type==2)					 //natural
	{
		L[n]=1;
		z[n]=0;
		c[n]=0;
	}
	else							// clamped
	{
		L[n]=h[n-1]*(2-u[n-1]);
		z[n]=(e[n]-h[n-1]*z[n-1])*1.0/L[n];
		c[n]=z[n];
	}
	
	for(int j=n-1;j>=0;j--)
	{
		c[j]=z[j]-u[j]*c[j+1];
		b[j]=(a[j+1]-a[j])*1.0/h[j]-h[j]*(c[j+1]+2*c[j])/3.0;
		d[j]=(c[j+1]-c[j])*1.0/(3*h[j]);
	}
	
	for(int j=n;j>0;j--){
		a[j]=a[j-1];
		b[j]=b[j-1];
		c[j]=c[j-1];
		d[j]=d[j-1];
	 }

}
double S( double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[] )
{

	if ( (t>=x[0])&&(t<=x[n]) )             // range of t
	{
		
		
		int flag=1;
		double g=0; 
		int i;
		for( i=0; i<n;i++)
			if(t==x[i])
				return a[i+1];

		

		for( i=1; i<=n&&flag ;i++)
		{
			if (t<=x[i])                            // the intervals
				{	
					flag=0;
					g=t-x[i-1];
				}
		}

		i--;
		double temp=(((d[i]*g+c[i])*g+b[i])*g+a[i]);
		return temp;
	}
	
	else
		return Fmax;
		
}