Matlab中的pchiptx.m 修改为C语言 三次样条插值

 

应为实验室项目原因 所以改写了pchiptx.m这个函数/*修改Matlab中的pchiptx.m三次样条插值*/#include <iostream>#include <math.h>using namespace std;float XT[100],YT[100],X[10000],Y[10000];//XT YT 已知的点集合 X Y 需要插入点集的点集float d[100]={0};float c[100],b[100];int sign( float x){ if ( x > 0) { return 1; } else if ( x==0) { return 0; } else { return -1; }}float pchipend( float h1, float h2, float del1, float del2){ float d = ( ( 2 * h1 + h2) * del1 - h1 * del2) / ( h1 + h2); if( sign(d) != sign(del1)) { d = 0; } else if ( ( sign(del1) != sign(del2)) && abs( d) > abs( 3*del1)) { d = 3 * del1; } return d;}float pchiptx( float xt[], float yt[], float x, int pointNumber)//XT YT 已知的点集合用于计算两点之间的三项式 x是用于取得点的y值{ float diffX[100], diffY[100], v, s, w1, w2; int i, k ; for ( i = 0; i < pointNumber - 1; i++) { diffX[i] = xt[i+1] - xt[i]; diffY[i] = ( yt[i+1] - yt[i]) / diffX[i]; } for ( i = 0; i <= pointNumber-2; i++)//不确定循环次数??? { if ( sign( diffY[i]) * sign( diffY[i+1]) > 0 ) { k = i + 1; w1 = 2 * diffX[k] + diffX[k-1]; w2 = diffX[k] + 2 * diffX[k-1]; d[k] = ( w1 + w2) / ( w1 / diffY[k-1] + w2 / diffY[k]); } } d[0] = pchipend( diffX[0], diffX[1], diffY[0], diffY[1]); for( i = 2; i < pointNumber ; i++) { d[i] = pchipend( diffX[i-1], diffX[i-2], diffY[i-1], diffY[i-2]); } for( i = 0; i < pointNumber; i++ ) { c[i] = ( 3 * diffY[i] - 2 * d[i] -d[i+1] ) / diffX[i]; b[i] = ( d[i] - 2 * diffY[i] + d[i+1] ) / ( diffX[i]*diffX[i]); } for ( i = 0; i < pointNumber - 1; i++) { if( xt[i] <= x && xt[i+1] >= x) { k = i; break; } } s = x - xt[k]; v = yt[k] + s*(d[k] + s*(c[k] + s*b[k])); return v;} int main(){ /*float pchiptx( float xt[], float yt[], float x, int pointNumber); int sign( float x); float pchipend( float h1, float h2, float del1, float del2);*/ FILE *out; out = fopen("c:/11111.txt","w+"); XT[0] = 0.388379; YT[0] = 0.248085; XT[1] = 0.322661; YT[1] = 0.263178; XT[2] = 0.276978; YT[2] = 0.290349; XT[3] = 0.231787; YT[3] = 0.328792; XT[4] = 0.205338; YT[4] = 0.375406; XT[5] = 0.187918; YT[5] = 0.440494; XT[6] = 0.175814; YT[6] = 0.511144; XT[7] = 0.175994; YT[7] = 0.582219; float temp = YT[0]; for( int i = 0; ; i++) { if( temp > YT[7]) { break; } Y[i] = temp; X[i] = pchiptx( YT, XT, temp, 8); fprintf(out,"%f %f/n",X[i],Y[i]); temp += 0.001; } fclose(out);} 

其中 输入的第一组数据必须从小到大