最近寫算法程序，發(fā)現(xiàn)有一些底層的程序看著貌似簡(jiǎn)單，想要把效果做的很好還是很難，之前都沒總結(jié)過。以至于之前即使用過也還得重新再寫，寫個(gè)基于C語(yǔ)言的程序集錦，供參考。

1、三次樣條插值

基本原理這里就不講了，直接上程序，復(fù)制后可以直接使用。

#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);/*在這里就可以打印插值后的數(shù)據(jù)*/return 0; }

圖1 數(shù)值實(shí)驗(yàn)結(jié)果

持續(xù)更新中。。。

總結(jié)

以上是生活随笔為你收集整理的简单算法集锦（程序）的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。