C语言 - 实现三次样条插值(Spline)

 

本文代码主要参考于百度知道：金色潜鸟 - 百度知道的回答。

他的程序代码有一个很小但较为严重的bug，所以在这对代码修正。

完整代码在文末，欢迎交流讨论。

---

说明：

文中以及代码中，

x、y分别代表插值前的横纵坐标；

u、s分别代表插值后的横纵坐标。

---

我在测试中，发现他的程序中有bug，即seval()函数中，给 j 赋的初值是插值后横坐标序列（u）的长度，会导致 u 长度超过一定值后，s数组的值会产生错误。（我测试u长度为21时会出错）

正确的方法：应该是给 j 赋值插值前横坐标序列（x）长度

故应seval()函数中，添加插值前横坐标序列（x）长度，即：

double seval(int n, double u,
	double x[], double y[],
	double b[], double c[], double d[],
	int *last);

改为

double seval(int ni, double u,
	int n, double x[], double y[],
	double b[], double c[], double d[],
	int *last);

 

 

完整代码如下： 

#include "stdafx.h"
#include "spline.h"
#include <math.h>


static int spline(	int n, int end1, int end2,
					double slope1, double slope2,
					double x[], double y[],
					double b[], double c[], double d[],
					int *iflag)
{
	int nm1, ib, i, ascend;
	double t;
	nm1 = n - 1;
	*iflag = 0;
	if (n < 2)
	{ /* no possible interpolation */
		*iflag = 1;
		goto LeaveSpline;
	}
	ascend = 1;
	for (i = 1; i < n; ++i) if (x[i] <= x[i - 1]) ascend = 0;
	if (!ascend)
	{
		*iflag = 2;
		goto LeaveSpline;
	}
	if (n >= 3)
	{
		d[0] = x[1] - x[0];
		c[1] = (y[1] - y[0]) / d[0];
		for (i = 1; i < nm1; ++i)
		{
			d[i] = x[i + 1] - x[i];
			b[i] = 2.0 * (d[i - 1] + d[i]);
			c[i + 1] = (y[i + 1] - y[i]) / d[i];
			c[i] = c[i + 1] - c[i];
		}
		/* ---- Default End conditions */
		b[0] = -d[0];
		b[nm1] = -d[n - 2];
		c[0] = 0.0;
		c[nm1] = 0.0;
		if (n != 3)
		{
			c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0]);
			c[nm1] = c[n - 2] / (x[nm1] - x[n - 3]) - c[n - 3] / (x[n - 2] - x[n - 4]);
			c[0] = c[0] * d[0] * d[0] / (x[3] - x[0]);
			c[nm1] = -c[nm1] * d[n - 2] * d[n - 2] / (x[nm1] - x[n - 4]);
		}
		/* Alternative end conditions -- known slopes */
		if (end1 == 1)
		{
			b[0] = 2.0 * (x[1] - x[0]);
			c[0] = (y[1] - y[0]) / (x[1] - x[0]) - slope1;
		}
		if (end2 == 1)
		{
			b[nm1] = 2.0 * (x[nm1] - x[n - 2]);
			c[nm1] = slope2 - (y[nm1] - y[n - 2]) / (x[nm1] - x[n - 2]);
		}
		/* Forward elimination */
		for (i = 1; i < n; ++i)
		{
			t = d[i - 1] / b[i - 1];
			b[i] = b[i] - t * d[i - 1];
			c[i] = c[i] - t * c[i - 1];
		}
		/* Back substitution */
		c[nm1] = c[nm1] / b[nm1];
		for (ib = 0; ib < nm1; ++ib)
		{
			i = n - ib - 2;
			c[i] = (c[i] - d[i] * c[i + 1]) / b[i];
		}
		b[nm1] = (y[nm1] - y[n - 2]) / d[n - 2] + d[n - 2] * (c[n - 2] + 2.0 * c[nm1]);
		for (i = 0; i < nm1; ++i)
		{
			b[i] = (y[i + 1] - y[i]) / d[i] - d[i] * (c[i + 1] + 2.0 * c[i]);
			d[i] = (c[i + 1] - c[i]) / d[i];
			c[i] = 3.0 * c[i];
		}
		c[nm1] = 3.0 * c[nm1];
		d[nm1] = d[n - 2];
	}
	else
	{
		b[0] = (y[1] - y[0]) / (x[1] - x[0]);
		c[0] = 0.0;
		d[0] = 0.0;
		b[1] = b[0];
		c[1] = 0.0;
		d[1] = 0.0;
	}
LeaveSpline:
	return 0;
}


static double seval(int ni, double u,
					int n, double x[], double y[],
					double b[], double c[], double d[],
					int *last)
{
	int i, j, k;
	double w;
	i = *last;
	if (i >= n - 1) i = 0;
	if (i < 0) i = 0;
	if ((x[i] > u) || (x[i + 1] < u))//??
	{
		i = 0;
		j = n;
		do
		{
			k = (i + j) / 2;
			if (u < x[k]) j = k;
			if (u >= x[k]) i = k;
		} while (j > i + 1);
	}
	*last = i;
	w = u - x[i];
	w = y[i] + w * (b[i] + w * (c[i] + w * d[i]));
	return (w);
}


void SPL(int n, double *x, double *y, int ni, double *xi, double *yi)
{
	double *b, *c, *d;
	int iflag=0, last=0, i=0;
	b = (double *)malloc(sizeof(double) * n);
	c = (double *)malloc(sizeof(double) * n);
	d = (double *)malloc(sizeof(double) * n);
	if (!d) { printf("no enough memory for b,c,d\n"); }
	else {
		spline(n, 0, 0, 0, 0, x, y, b, c, d, &iflag);
		if (iflag == 0) 
			printf("I got coef b,c,d now\n"); 
		else 
			printf("x not in order or other error\n");
		for (i = 0; i<ni; i++) 
			yi[i] = seval(ni, xi[i], n, x, y, b, c, d, &last);

		free(b); 
		free(c); 
		free(d);
	};
}

int main() 
{
    /* Sin插值 */
	double x[7] = { 0, 1, 2, 3, 4, 5, 6 };
	double y[7] = { 0, 0.841470984807897, 0.909297426825682, 0.141120008059867, -0.756802495307928, -0.958924274663139, -0.279415498198926 };
	double u[61] = { 0 };
	double s[61] = { 0 };
	int i;

	for (i = 0; i < 61; i++)
	{
		u[i] = (double)i*0.1;
	}
	
	SPL(7, x, y, 61, u, s);

    /*
        在这里就可以打印插值后的数据
    */

	return 0;
}

---

C语言运行结果图：