spline.c

C语言编写

#include <stdio.h>
#include <stdlib.h>
#include<math.h>


void MySpline(int N, double* x,     double* y, 
			           int R, double* xInt, double* yInt);

void chase(int N, double* b, double* c, double* d, double* a, double *p );

void main()
{
	int i;
	int N=10; //节点数
	int R = 9; //待插入点数
	double* x;    //已知节点和其函数值
	double* y;
	double* xInt;  //待插入点
	double* yInt;
	double Epsilon;           //精度

	Epsilon = 1.0e-6;
	x = (double*)calloc(N,sizeof(double));
	y = (double*)calloc(N,sizeof(double));
    xInt = (double*)calloc(R,sizeof(double));
	yInt = (double*)calloc(R,sizeof(double));
	for (i = 0; i < N; i++)
	{
		x[i]=0.1*i;
		y[i]=2.36*x[i]*x[i];
		if (i<R)
		{
			xInt[i] = x[i]+0.025;
		}
		
	}
	for (i = 0; i <N; i++)
	{
		printf("%10.6f",x[i]);
	}
	printf("\n");
	for (i = 0; i <N; i++)
	{
		printf("%10.6f",y[i]);
	}
	printf("\n\n");
	for (i = 0; i<R; i++)
	{
		printf("%10.6f",xInt[i]);
	}
	printf("\n");
	for (i = 0; i<R; i++)
	{
		printf("%10.6f",2.36*xInt[i]*xInt[i]);
	}


	MySpline(N-1,x,y,R,xInt,yInt);
	printf("\n");
	for (i = 0; i<R; i++)
	{
		printf("%10.6f",yInt[i]);
	}
	printf("\n");

}



/*--------------------------------------
三弯矩方程
Function: 针对自然边界条件,采用追赶法解三弯矩方程的三次B样条插值
Name:     MySpline
Input:      N:N+1个节点;             x:已知节点;               y:已知节点函数值
               R:R个待插值节点   xInt:待插值节点        yInt:返回插值后的函数值 
Output:    void
---------------------------------------*/
void MySpline(int N, double* x,     double* y, 
			           int R, double* xInt, double* yInt)
{
	int i,k;

	double* h;
	double* Mu;
	double* Beta;
	double* g;
	double* M;
	double dTemp;
	double* b;



	h = (double*)calloc(N+1,sizeof(double));
	Mu = (double*)calloc(N,sizeof(double));
	Beta = (double*)calloc(N,sizeof(double));
	g = (double*)calloc(N+1,sizeof(double));
	M = (double*)calloc(N+1,sizeof(double));
	b = (double*)calloc(N+1,sizeof(double));

	for (i = 0;i <=N; i++)
	{
		b[i] =2;
	}

	for (i = 1; i <= N; i++)
	{
		h[i] = x[i] - x[i-1];
	}
	//一阶导数连续, N-1个方程
	for (i = 1; i <=N-1; i++)
	{
		Beta[i] = h[i+1]/(h[i]+h[i+1]);
		Mu[i-1] = 1-Beta[i];
		dTemp = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
		g[i] = 6*dTemp/(h[i]+h[i+1]);
	}

	//自然边界条件,2个方程
	g[0] = 2.36;
	g[N] = 2.36;

    Beta[0] = 0;
	Mu[N-1] = 0;


	//追赶法解三弯矩方程
	chase(N+1,b,Beta,Mu,M,g);

	//求值
	for (i = 0; i<R; i++)
	{
		if ((xInt[i] < x[0]) ||(xInt[i] > x[N]))
		{
			printf("Wrong Interpose.\n");
			exit(1);
		}
		k = 0;
		while(xInt[i] >x[k] )
		{
			k++;
		}
		yInt[i] = M[k-1]*pow((x[k]-xInt[i]),3)/6/h[k]+M[k]*pow((xInt[i]-x[k-1]),3)/6/h[k]
			+(y[k-1]-M[k-1]*h[k]*h[k]/6)*(x[k]-xInt[i])/h[k]
			+(y[k]-M[k]*h[k]*h[k]/6)*(xInt[i]-x[k-1])/h[k];
	}
}


/*----------------------------------------------
Function:  追赶法解三弯矩方程
Name:       chase
Input:        N:方程阶数;  b:对角元素;  c:上对角元素; d:下对角元素
                  a: 待求方程的解;   p:方程组右值
----------------------------------------------*/
void chase(int N, double* b, double* c, double* d, double* a, double *p )
{
	int i;
	
	double* l;
	double* s;
	double* v;
	
	l = (double*)calloc(N,sizeof(double));
	s = (double*)calloc(N-1,sizeof(double));
	v = (double*)calloc(N,sizeof(double));
	
	l[0] = b[0];
	for (i = 0; i<=N-2;i++)
	{
		s[i]=c[i]/l[i];
		l[i+1]=b[i+1]-d[i]*s[i];
	}
	
	v[0] = p[0]/l[0];             //追
	for (i = 1; i<=N-1;i++)
	{
		v[i] = (p[i]-d[i-1]*v[i-1])/l[i];
	}
	
	a[N-1]=v[N-1];                          //赶
	for (i = N-2;i>=0; i--)
	{
		a[i] = v[i]-s[i]*a[i+1];
	}
	
	
}