C三次样条插值函数(自然边界)

/****************函数说明*********************///pxpy为已知的数据点，xs为要插值的x坐标，最终会得到xs坐标下的y值usingSystem;usingSystem.Collections.Generic;usingSystem.Linq;usingSystem.Text;namespacespline{classProgram{staticvoidMain(string[]args){point[]points=newpoint[13];double[]px={64,304,544,1035,1502,2061,2540,2939,3498,3897,4456,4696,4936};double[]py={4.663,4.476,4.969,4.853,4.766,5.2415,4.795,4.7055,5.4565,5.5515,5.411,5.8385,6.7085};for(inti=0;i13;i++){points[i]=newpoint();points[i].x=px[i];points[i].y=py[i];}point.DeSortX(points);double[]xs={64,304,544,1023,1502,1981,2460,2939,3418,3897,4376,4616,4856};splineInsertPoint(points,xs,1);Console.ReadLine();}staticdouble[]splineInsertPoint(point[]points,double[]xs,intchf){intplength=points.Length;double[]h=newdouble[plength];double[]f=newdouble[plength];double[]l=newdouble[plength];double[]v=newdouble[plength];double[]g=newdouble[plength];//三转角法的待定一阶系数法。基本方程为l(i)m(i-1)+2m(i)+u(i)m(i+1)=g(i),i=1,2,...,n-1//目前有n-1个方程，n+1个未知量（从m0到mn），需要添加两个边界条件for(inti=0;iplength-1;i++)//共有plength(n+1)个点，h总有plength-1(n)个值{h[i]=points[i+1].x-points[i].x;}for(inti=1;iplength-1;i++){l[i]=h[i]/(h[i-1]+h[i]);v[i]=h[i-1]/(h[i-1]+h[i]);g[i]=3*(l[i]*f[i-1]+v[i]*f[i]);//i=1,2,3,...,n-1}//求解矩阵，用追赶法double[]U=newdouble[plength];//上三角矩阵U,plength=n+1阶矩阵double[]L=newdouble[plength];//下三角矩阵Ldouble[]m=newdouble[plength];//用于保存m0-mn的一阶系数数值double[]Y=newdouble[plength];//中间矩阵Y/***********************矩阵分界，求解待定的一阶系数，m1到m(n-1)*************************/U[0]=2;v[0]=1;for(inti=1;iplength;i++){L[i]=l[i]/U[i-1];U[i]=2-v[i-1]*L[i];}Y[0]=g[0];for(inti=1;iplength;i++){Y[i]=g[i]-L[i]*Y[i-1];}m[plength-1]=Y[plength-1]/U[plength-1];for(inti=plength-2;i=0;i--){m[i]=(Y[i]-v[i]*m[i+1])/U[i];}/*****************************************************************************************/intxlength=xs.Length;//记录需要插值的点的个数double[]insertRes=newdouble[xlength];//用于返回插值点的Y坐标for(inti=0;ixlength;i++){intj=0;for(j=0;jplength;j++){if(xs[i]points[j].x)break;}j=j-1;//Console.WriteLine(j);//表示要插值的点数组的坐标在第j+1组区间if(j==-1||j==plength-1){if(j==-1)thrownewException(插值下边界超出);if(j==plength-1&&xs[i]==points[j].x)//这种情况是，需要插值的点正好等于一系列点x坐标最大的那个点insertRes[i]=points[j].y;elsethrownewException(插值下边界超出);}else{//根据每个分段的三次多项式求值insertRes[i]=(h[j]+2*(xs[i]-points[j].x))*Math.Pow(xs[i]-points[j+1].x,2)*points[j].y/Math.Pow(h[j],3)+(h[j]+2*(points[j+1].x-xs[i]))*Math.Pow(xs[i]-points[j].x,2)*points[j+1].y/Math.Pow(h[j],3)+(xs[i]-points[j].x)*Math.Pow(xs[i]-points[j+1].x,2)*m[j]/Math.Pow(h[j],2)+(xs[i]-points[j+1].x)*Math.Pow(xs[i]-points[j].x,2)*m[j+1]/Math.Pow(h[j],2);Console.WriteLine(f(+xs[i]+)=+insertRes[i]);}}returninsertRes;}}classpoint{publicdoublex=0;publicdoubley=0;publicpoint(){x=0;y=0;}//-------写一个排序函数，使得输入的点按顺序排列，是因为插值算法的要求是，x轴递增有序的---------publicstaticpoint[]DeSortX(point[]points){intlength=points.Length;doubletemx,temy;for(inti=0;ilength-1;i++){for(intj=0;jlength-i-1;j++)if(points[j].xpoints[j+1].x){temx=points[j+1].x;points[j+1].x=points[j].x;points[j].x=temx;temy=points[j+1].y;points[j+1].y=points[j].y;points[j].y=temy;}}returnpoints;}}}