三次样条插值的C程序(很全啊)

三次样条插值C/C++程序（自己整理的）具体推导看书<<数值分析>>code：#include <iostream>using namespace std;const int MAXN = 100;int n;double x[MAXN], y[MAXN]; //下标从0..ndouble alph[MAXN], beta[MAXN], a[MAXN], b[MAXN];double h[MAXN];double m[MAXN]; //各点的一阶导数;inline double sqr(double pa) {return pa * pa;}double sunc(double p, int i) {return (1 + 2 * (p - x[i]) / (x[i + 1] - x[i])) * sqr((p - x[i + 1]) / (x[i + 1] - x[i])) * y[i] + (1 + 2 * (p - x[i + 1]) / (x[i] - x[i + 1])) * sqr((p - x[i]) / (x[i + 1] - x[i])) * y[i + 1] + (p - x[i]) * sqr((p - x[i + 1]) / (x[i] - x[i + 1])) * m[i]+ (p - x[i + 1]) * sqr((p - x[i]) / (x[i + 1] - x[i])) * m[i + 1];}int main() {int i, j;freopen("threeInsert.in", "r", stdin);scanf("%d", &n);for (i = 0; i <= n; i++) scanf("%lf%lf", &x[i], &y[i]);// scanf("%lf%lf", &m[0], &m[n]);for (i = 0; i <= n - 1; i++) h[i] = x[i + 1] - x[i];//第一种边界条件//alph[0] = 0; alph[n] = 1; beta[0] = 2 * m[0]; beta[n] = 2 * m[n];//第二种边界条件alph[0] = 1; alph[n] = 0; beta[0] = 3 * (y[1] - y[0]) / h[0]; beta[n] = 3 * (y[n] - y[n - 1] / h[n - 1]);for (i = 1; i <= n - 1; i++) {alph[i] = h[i - 1] / (h[i - 1] + h[i]);beta[i] = 3 * ((1 - alph[i]) * (y[i] - y[i - 1]) / h[i - 1] + alph[i] * (y[i + 1] - y[i]) / h[i]);}a[0] = - alph[0] / 2; b[0] = beta[0] / 2;for (i = 1; i <= n; i++) {a[i] = - alph[i] / (2 + (1 - alph[i]) * a[i - 1]);b[i] = (beta[i] - (1 - alph[i]) * b[i - 1]) / (2 + (1 - alph[i]) * a[i - 1]);}m[n + 1] = 0;for (i = n; i >= 0; i--) {m[i] = a[i] * m[i + 1] + b[i];}scanf("%lf", &xx);for (i = 0; i < n; i++) {if (xx >= x[i] && xx <= x[i + 1]) break;}printf("%lf\n", sunc(xx, i));}#include<iostream>#include<iomanip>using namespace std;const int MAX = 50;float x[MAX], y[MAX], h[MAX];//变量设置：x为各点横坐标；y为各点纵坐标；h为步长float c[MAX], a[MAX], fxym[MAX];float f(int x1, int x2, int x3)/*****************求差分函数(含三个参数)****************************/{float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);return (a - b)/(x[x3] - x[x1]);}void cal_m(int n)/***********************用追赶法求解出弯矩向量M……***************************/{float B[MAX];B[0] = c[0] / 2;for(int i = 1; i < n; i++)B[i] = c[i] / (2 - a[i]*B[i-1]);//fxym[0] = fxym[0] / 2;for(i = 1; i <= n; i++)fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);for(i = n-1; i >= 0; i--)fxym[i] = fxym[i] - B[i]*fxym[i+1];}void printout(int n);int main(){int n,i; char ch;do{cout<<"请输入已知断点个数:";cin>>n;for(i = 0; i <= n; i++){cout<<"Please put in X"<<i<<':';cin>>x[i]; //cout<<endl;cout<<"Please put in Y"<<i<<':';cin>>y[i]; //cout<<endl;}for(i = 0; i < n; i++) //求步长；其数组值较之x，y个数少一h[i] = x[i+1] - x[i];cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";int t;float f0, f1;cin>>t;switch(t){case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";//显示数据为Y0'至Yn'，即断点的一阶导数cin>>f0>>f1;c[0] = 1; a[n] = 1;fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];break;case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";//显示数据为Y0"至Yn"，即断点的二阶导数cin>>f0>>f1;c[0] = a[n] = 0;fxym[0] = 2*f0; fxym[n] = 2*f1;break;default:cout<<"不可用\n";//待定};//switchfor(i = 1; i < n; i++)fxym[i] = 6 * f(i-1, i, i+1);//调用差分函数(only!) for(i = 1; i < n; i++){a[i] = h[i-1] / (h[i] + h[i-1]);c[i] = 1 - a[i];}a[n] = h[n-1] / (h[n-1] + h[n]);cal_m(n);//调用弯矩函数（only！）cout<<"\n输出三次样条插值函数：\n";printout(n);//调用求解三次样条插值函数；函数输出cout<<"Do you to have anther try ? y/n :";cin>>ch;}while(ch == 'y' || ch == 'Y');return 0;}void printout(int n)/***************求三次样条插值函数（因已知断点个数而异）***********************/{cout<<setprecision(6);//通过操作器setprecision()设置有效位数；其为头文件<iomanip.h>所包含；括号内为参数。

for(int i = 0; i < n; i++)//所输出函数个数由所设断点个数而定{cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";float t = fxym[i]/(6*h[i]);if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";elsecout<<-t<<"*(x - "<<x[i+1]<<")^3";t = fxym[i+1]/(6*h[i]);if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";elsecout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";cout<<"\n\t";t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";elsecout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";elsecout<<" - "<<-t<<"*(x - "<<x[i]<<")";cout<<endl<<endl;}cout<<endl;}/****************提供测试数据：（来自课本55页例5 《数值分析》清华大学出版社第四版）***********//*输入327.7 4.128 4.329 4.130 3.013.0 -4.0*//*输出输出三次样条插值函数：1: [27.7 , 28]13.07*(x - 28)^3 + 0.22*(x - 27.7)^3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: [28 , 29]0.066*(29 - x)^3 + 0.1383*(x - 28)^3+ 4.234*(29 - x) + 3.962*(x - 28)3: [29 , 30]0.1383*(30 - x)^3 - 1.519*(x - 29)^3+ 3.962*(30 - x) + 4.519*(x - 29)*////////////////////////////////////////////////////// purpose：给定,值的三次样条插值多项式/////////////////////////////////////////////////////# define MAX_N 20 // 定义(x_i,y_i)的最大的维数typedef struct tagPOINT // 点的结构{double x;double y;} POINT;int main ( ){int n;int i,k;POINT points[MAX_N +1];double h[MAX_N +1],b[MAX_N +1],c[MAX_N +1],d[MAX_N +1],M[MAX_N +1]; double u[MAX_N +1],v[MAX_N +1],y[MAX_N +1];double x,p,q,S;printf ("\nInput n value："); // 输入插值点的数目scanf ("%d",&n);if (n>MAX_N){printf ("The input n is larger than MAX_N, please redefine the MAX_N. \n"); return 1;}if (n<=0){printf ("Please input a number between 1 and % d. \n",MAX_N);return 1;}// 输入插值点(x_i,y_i),M0值和Mn值printf ("Now input the (x_i,y_i),i=0,…,% d：\n",n);for (i=0;i<n;i++){scanf ("%lf %lf",&points[i].x,&points[i].y);printf("%lf %lf\n",points[i].x,points[i].y);}printf ("Now input the M[0] value：");scanf ("%lf",&M[0]);printf ("\nNow input the M[n] value：");scanf ("%lf",&M[n]);printf ("\nNow input the x value：");//输入计算三次样条插值函数的x值scanf ("%lf",&x);if (x>points[n-1]. x || x<points[0].x){printf ("Please input a number between %f and %f. \n",points[0].x,points[n].x); return 1;}// 计算M关系式中各参数的值h[0]=points[1].x-points[0].x;for (i=1;i<n;i++ ){h[i]=points[i+1].x-points[i].x;b[i]=h[i]/(h[i]+h[i-1]);c[i]=1-b[i];d[i]=6*((points[i+1].y-points[i].y)/h[i]-(points[i].y-points[i-1].y)/h[i-1])/(h[i]+h[i-1]);}// 用追赶法计算Mi,i=1,…,n-1d[1]-= c[1]*M[0];d[n-1]-=b[n-1]*M[n];b[n-1]=0;c[1]=0;v[0]=0;for ( i=1;i<n;i++){u[i]=2-c[i]*v[i-1];v[i]=b[i]/u[i];y[i]= (d[i]-c[i]*y[i-1])/u[i];}for (i=1;i<n;i++){M[n-i]=y[n-i]-v[n-i]* M[n-i+1];}// 计算三次样条插值函数在x处的值k=0;while (x>= points[k].x) k++;k=k-1;p=points[k+1].x-x;q=x-points[k].x;S=(p* p* p* M[k] +q* q* q* M[k+1]) / (6* h[k])+( p* points [k].y+q* points[k+1].y) / h[k]-h[k]*(p* M[k] +q* M[k]+q* M[k+1])/6; printf ("S(%f ) = %f\n",x,S);// 输出getchar ( );return 0;}#include"iostream.h"#include"fstream.h"#define N 100float X[N],Y[N],M[N],a,b;double f(int i,int j){double sum1=0;if(j-i==1) return (Y[j]-Y[i])/(X[j]-X[i]);else return (f(i+1,j)-f(i,j-1))/(X[j]-X[i]);}void main(){ifstream istrm("E:\\data.txt");//这里是文件的输入，请相应改成普通键盘输入int n=3;for(int i=0;i<=n;i++)istrm>>X[i]>>Y[i];istrm>>a>>b;float h[N],u[N],q[N],d[N];for(i=0;i<n;i++)h[i]=X[i+1]-X[i];q[0]=1;q[n]=0;for(i=1;i<n;i++)q[i]=h[i]/(h[i-1]+h[i]);for(i=0;i<=n;i++)u[i]=1-q[i];d[0]=(6*(f(0,1)-a))/h[0];d[n]=(6*(b-f(n-1,n)))/h[n-1];for(i=1;i<n;i++)d[i]=6*f(i-1,i+1);float m[N];for(i=0;i<=n;i++) m[i]=2;for(i=1;i<=n;i++){m[i]=m[i]-(u[i]/m[i-1])*q[i-1]; d[i]=d[i]-(u[i]/m[i-1])*d[i-1];}M[n]=d[n]/m[n];cout<<M[n]<<endl;for(i=n-1;i>=0;i--){M[i]=(d[i]-q[i]*M[i+1])/m[i]; cout<<M[i]<<endl;}}。