三次样条插值

  1 /**********************************************************
  2 ** 三次样条插值函数.C
  3 ** 进行自然边界，夹持边界，非扭结边界条件下的插值
  4 **
  5 ** 2012-05-25 将函数从原来的需要至少4个插值点才能计算
  6 **            扩展成只需要2个插值点就可以完成计算
  7 **            其他一些改变
  8 ** 2008-12-27 创建函数
  9 **********************************************************/
 10 #include <stdio.h>
 11 #include <stdlib.h>
 12 
 13 #define ABS(p) ( ((p) > 0) ? (p) : -(p) )
 14 #define SWAP(x, y, temp) (temp) = (y); (y) = (x); (x) = (temp);
 15 
 16 typedef enum _condition
 17 {
 18     NATURAL, // 自然边界
 19     CLAMPED, // 夹持边界
 20     NOTAKNOT // 非扭结边界
 21 }condition;
 22 
 23 typedef struct _coefficient
 24 {
 25     double a3;
 26     double b2;
 27     double c1;
 28     double d0;
 29 }coefficient;
 30 
 31 typedef struct _point
 32 {
 33     coefficient *coe;    // 插值结果系数矩阵
 34     double *xCoordinate; // 插值点横坐标
 35     double *yCoordinate; // 插值点纵坐标
 36     double f0;           // 在夹持条件下的最左边点的导数值
 37     double fn;           // 在夹持条件下的最右边点的导数值
 38     int num;             // 插值点数
 39     condition con;       // 边界条件
 40 }point;
 41 
 42 
 43 /*
 44 ** 资源不足时函数返回 -1
 45 ** 插值点数小于2时，函数返回 -2
 46 ** 计算正确完成函数返回插值点的数量 n
 47 */
 48 int spline(point *point)
 49 {
 50     double *x = (*point).xCoordinate, *y = (*point).yCoordinate;
 51     int n = (*point).num - 1; // 循环上限
 52     int i;                    // 循环辅助变量
 53     coefficient *coe = (*point).coe;
 54     condition con = (*point).con;
 55     double *h, *d;
 56     double *a, *b, *c, *f, *m;
 57     double temp;
 58 
 59     if (n < 1)
 60     {return -2;}
 61 
 62     h = (double *)malloc( (6 * n + 4) * sizeof(double) ); /* 0,       1,...,(n-1)          */
 63     if (h == NULL)
 64     {return -1;}
 65     d = h + n;                                            /* 0,       1,...,(n-1)          */
 66     a = d + n;                                            /* 特别使用,1,...,(n-1),       n */
 67     b = a + (n + 1);                                      /*        0,1,...,(n-1),       n */
 68     c = b + (n + 1);                                      /*        0,1,...,(n-1),特别使用 */
 69     f = c + (n + 1);                                      /*        0,1,...,(n-1),       n */
 70     m = b;
 71 
 72 
 73     /* 计算 h[] d[] */
 74     for (i = 0; i < n; i++)
 75     {
 76         h[i] = x[i + 1] - x[i];
 77         d[i] = (y[i + 1] - y[i]) / h[i];
 78         /* printf("%f\t%f\n", h[i], d[i]); */
 79     }
 80     /**********************
 81     ** 初始化系数增广矩阵
 82     **********************/
 83     a[0] = (*point).f0;
 84     c[n] = (*point).fn;
 85     /* 计算 a[] b[] c[] f[] */
 86     switch(con)
 87     {
 88     case NATURAL:
 89         f[0] = 0;
 90         f[n] = 0;
 91         a[0] = 0;
 92         c[n] = 0;
 93         c[0] = 0;
 94         a[n] = 0;
 95         b[0] = 1;
 96         b[n] = 1;
 97         break;
 98 
 99     case CLAMPED:
100         f[0] = 6 * (d[0] - a[0]);
101         f[n] = 6 * (c[n] - d[n - 1]);
102         a[0] = 0;
103         c[n] = 0;
104         c[0] = h[0];
105         a[n] = h[n - 1];
106         b[0] = 2 * h[0];
107         b[n] = 2 * h[n - 1];
108         break;
109 
110     case NOTAKNOT:
111         f[0] = 0;
112         f[n] = 0;
113         a[0] = h[0];
114         c[n] = h[n - 1];
115         c[0] = -(h[0] + h[1]);
116         a[n] = -(h[n - 2] + h[n - 1]);
117         b[0] = h[1];
118         b[n] = h[n - 2];
119         break;
120     }
121 
122     for (i = 1; i < n; i++)
123     {
124         a[i] = h[i - 1];
125         b[i] = 2 * (h[i - 1] + h[i]);
126         c[i] = h[i];
127         f[i] = 6 * (d[i] - d[i - 1]);
128     }
129     /* for (i = 0; i < n+1; i++)   printf("%f\n", c[i]); */
130 
131     /*************
132     ** 高斯消元
133     *************/
134     if (n > 2)
135     {
136         /* 第0行到第(n-3)行的消元 */
137         for(i = 0; i <= n - 3; i++)
138         {
139             /* 选主元 */
140             if ( ABS(a[i + 1]) > ABS(b[i]) )
141             {
142                 SWAP(a[i + 1], b[i], temp);
143                 SWAP(b[i + 1], c[i], temp);
144                 SWAP(c[i + 1], a[i], temp);
145                 SWAP(f[i + 1], f[i], temp);
146             }
147             temp = a[i + 1] / b[i];
148             a[i + 1] = 0;
149             b[i + 1] = b[i + 1] - temp * c[i];
150             c[i + 1] = c[i + 1] - temp * a[i];
151             f[i + 1] = f[i + 1] - temp * f[i];
152         }
153     }
154     if (n >= 2)
155     {
156         /* 最后3行的消元 */
157         /* 第(n-2)行的消元 */
158         /* 选主元 */
159         if ( ABS(a[n - 1]) > ABS(b[n - 2]) )
160         {
161             SWAP(a[n - 1], b[n - 2], temp);
162             SWAP(b[n - 1], c[n - 2], temp);
163             SWAP(c[n - 1], a[n - 2], temp);
164             SWAP(f[n - 1], f[n - 2], temp);
165         }
166         /* 选主元 */
167         if ( ABS(c[n]) > ABS(b[n - 2]) )
168         {
169             SWAP(c[n], b[n - 2], temp);
170             SWAP(a[n], c[n - 2], temp);
171             SWAP(b[n], a[n - 2], temp);
172             SWAP(f[n], f[n - 2], temp);
173         }
174         /* 消元 */
175         temp = a[n - 1] / b[n - 2];
176         a[n - 1] = 0;
177         b[n - 1] = b[n - 1] - temp * c[n - 2];
178         c[n - 1] = c[n - 1] - temp * a[n - 2];
179         f[n - 1] = f[n - 1] - temp * f[n - 2];
180         /* 消元 */
181         temp = c[n] / b[n - 2];
182         c[n] = 0;
183         a[n] = a[n] - temp * c[n - 2];
184         b[n] = b[n] - temp * a[n - 2];
185         f[n] = f[n] - temp * f[n - 2];
186     }
187     /* 最后2行的消元 */
188     /* 第(n-1)行的消元 */
189     /* 选主元 */
190     if ( ABS(a[n]) > ABS(b[n - 1]) )
191     {
192         SWAP(a[n], b[n - 1], temp);
193         SWAP(b[n], c[n - 1], temp);
194         SWAP(f[n], f[n - 1], temp);
195     }
196     /* 消元 */
197     temp = a[n] / b[n-1];
198     a[n] = 0;
199     b[n] = b[n] - temp * c[n - 1];
200     f[n] = f[n] - temp * f[n - 1];
201 
202     /******************
203     ** 回代求解 m[]
204     ******************/
205     m[n] = f[n] / b[n];
206     m[n - 1] = (f[n - 1] - c[n - 1] * m[n]) / b[n-1];
207     for (i = n - 2; i >= 0; i--)
208     {
209         m[i] = ( f[i] - (m[i + 2] * a[i] + m[i + 1] * c[i]) ) / b[i];
210     }
211     /* for (i = 0; i < n+1; i++)   printf("%f\n", m[i]); */
212 
213     /***********************
214     ** 计算插值函数所有参数
215     ***********************/
216     for (i = 0; i < n; i++)
217     {
218         coe[i].a3 = (m[i + 1] - m[i]) / (6 * h[i]);
219         coe[i].b2 = m[i] / 2;
220         coe[i].c1 = d[i] - (h[i] / 6) * (2 * m[i] + m[i + 1]);
221         coe[i].d0 = y[i];
222     }
223     free(h);
224     // 计算正确完成返回插值点的数量
225     return n + 1;
226 }
227 
228 
229 void main()
230 {
231     /* x[]内不能出现重复数据 */
232     double x[] = { 0, 0.5, 0.75, 1.25, 2.5,  5, 7.5, 10, 15,
233                   20,  25,   30,   35,  40, 45,  50, 55, 60,
234                   65,  70,   75,   80,  85, 90,  95, 100};
235     double y[] = {     0, 0.6107, 0.7424, 0.9470, 1.3074, 1.7773,    2.1,
236                   2.3414, 2.6726, 2.8688, 2.9706, 3.0009, 2.9743, 2.9015,
237                   2.7904, 2.6470, 2.4762, 2.2817, 2.0662, 1.8320, 1.5802,
238                   1.3116, 1.0263, 0.7239, 0.4033, 0.0630};
239     int n = sizeof(x) / sizeof(double);
240     coefficient *coe;
241     int i;
242     point p;
243     coe = (coefficient *)malloc((n - 1) * sizeof(coefficient));
244     p.xCoordinate = x;
245     p.yCoordinate = y;
246     /* 下面两个参数f0和fn只在夹持条件下有效，其他条件下这两个值会被忽略 */
247     p.f0 = 0;// 在夹持条件下的左边点的导数值
248     p.fn = 0;// 在夹持条件下的右边点的导数值
249     p.num = n;
250     p.con = NOTAKNOT;
251     p.coe = coe;
252     spline(&p);
253     
254     for (i = 0; i < n - 1; i++)
255         printf("%f\t%f\t%f\t%f\n", coe[i].a3, coe[i].b2, coe[i].c1, coe[i].d0);
256     free(coe);
257 }