如何通过python实现三次样条插值
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了如何通过python实现三次样条插值相关的知识，希望对你有一定的参考价值。
参考技术A spline函数可以实现三次样条插值 x = 0:10; y = sin(x); xx = 0:.25:10; yy = spline(x,y,xx); plot(x,y,'o',xx,yy) 另外fnplt csapi这两个函数也是三次样条插值函数，具体你可以help一下！本回答被提问者采纳三次样条插值
百度空间关掉了因此搬家到这里。
三次样条插值函数网上搜一搜其实很多了，但代码水平参差不齐，真正好用的其实没有几个。最难受的是那些代码还被当做示例被人无数次的学习。所以还是把这段代码放出来，希望后来者可以从中收益。
最后，如果你告诉我你把这段代码用在了什么地方，我会非常高兴。
这次对三次样条插值函数做了一些修改使他的适用性更为广泛。现在三次样条插值函数可以完成2个点以上的插值。
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 %f ", 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 ", 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 ", 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 %f %f %f ", coe[i].a3, coe[i].b2, coe[i].c1, coe[i].d0); 256 free(coe); 257 }
以上是关于如何通过python实现三次样条插值的主要内容，如果未能解决你的问题，请参考以下文章