三次样条插值算法c语言

三次样条插值2008年12月20日 星期六 19:44MATLAB的样条插值函数默认的设置，是一种称为 not-a-knot 的边界条件

not-a-knot（非扭结）边界条件，也就是强制使第一个点的三次导数和第二点的三次导数一样；最后一个点的三次导数和倒数第一个点一样。

费了半天的劲，自己写出了spline函数，具备与MATLAB一样的数值精度！

/**************************
** 三次样条插值函数.C
** 进行自然边界，
**     夹持边界，
**     非扭结边界
** 条件下的插值
** 2008-12-27 slx 保留所有权利
**************************/
#include 
#include 

typedef enum _condition
{
    NATURAL
    CLAMPED
    NOTaKNOT
}condition;

typedef struct _coefficient
{
    double a3;
    double b2;
    double c1;
    double d0;
}coefficient;

typedef struct _point
{
    coefficient *coe;
    double *xCoordinate;
    double *yCoordinate;
    double f0;
    double fn;
    int num;
    condition con;
}point;


int spline（point *point）
{
    double *x = （*point）.xCoordinate *y = （*point）.yCoordinate;
    int n = （*point）.num - 1;
    coefficient *coe = （*point）.coe;
    condition con = （*point）.con;
    double *h *d;
   double *a *b *c *f *m;
    double temp;
    int i;
    h = （double *）malloc（ n * sizeof（double） ）; /* 01--（n-1）n */
    d = （double *）malloc（ n * sizeof（double） ）; /* 01--（n-1）n */
    
    a = （double *）malloc（ （n + 1） * sizeof（double） ）; /* 特别使用1--（n-1）       n */
    b = （double *）malloc（ （n + 1） * sizeof（double） ）; /*        01--（n-1）       n */
    c = （double *）malloc（ （n + 1） * sizeof（double） ）; /*        01--（n-1）特别使用 */
    f = （double *）malloc（ （n + 1） * sizeof（double） ）; /*        01--（n-1）       n */
    m = b;
    if（f == NULL）
    {
        free（h）;
        free（d）;
        free（a）;
        free（b）;
        free（c）;
        return 1;
    }
    /* 计算 h[] d[] */
    for （i = 0; i < n; i++）
    {
        h[i] = x[i + 1] - x[i];
        d[i] = （y[i + 1] - y[i]） / h[i];
        /* printf（“%f\t%f\n“ h[i] d[i]）; */
    }
    /**********************
    ** 初始化系数增广矩阵
    **********************/
    a[0] = （*point）.f0;
    c[n] = （*point）.fn;
    /* 计算 a[] b[] d[] f[] */
    switch（con）
    {
        case NATURAL:
            f[0] = a[0];
            f[n] = c[n];
            a[0] = 0;
            c[n] = 0;
            c[0] = 0;
            a[n] = 0;
            b[0] = 1;
            b[n] = 1;
        break;
        
        case CLAMPED:
            f[0] = 6 * （d[0] - a[0]）;
            f[n] = 6 * （c[n] - d[n - 1]）;
            a[0] = 0;
            c[n] = 0;
            c[0] = h[0];
            a[n] = h[n - 1];
            b[0] = 2 * h[0];
            b[n] = 2 * h[n - 1];
        break;
        
        case NOTaKNOT:
            f[0] = 0;
            f[n] = 0;
            a[0] = h[0];
            c[n] = h[n - 1];
            c[0] = -（h[0] + h[1]）;
            a[n] = -（h[n - 2] + h[n - 1]）;
            b[0] = h[1];
            b[n] = h[n - 2];
        break;
    }

    for （i = 1; i < n; i++）
    {
        a[i] = h[i - 1];
        b[i] = 2 * （h[i - 1] + h[i]）;
        c[i] = h[i];
        f[i]