三次样条插值
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了三次样条插值相关的知识，希望对你有一定的参考价值。
条件
（1）输入$x_{i},y_{i}=f(x_{i}),0\leq i\leq n$
（2）要求拟合的曲线$S(x)$满足:对于任意的$1\leq i\leq n-1$，在$x_{i}$处一阶二阶导数连续，$S(x)$ 也连续，且$S^{‘}(x_{0})=f^{‘}(x_{0})$,$S^{‘}(x_{n})=f^{‘}(x_{n})$
求解过程
设$S_{x_{j}}^{‘‘}=M_{j}$。对于区间$[x_{j},x_{j+1}]$,$S(x)$是$[x_{j},x_{j+1}]$上的线性函数，所以设$S^{‘‘}=M_{j}\frac{x_{j+1}-x}{h_{j}}+M_{j+1}\frac{x-x_{j}}{h_{j}}$,$h_{j}=x_{j+1}-x_{j}$。对$S^{‘‘}$两次积分，并利用$S(x_{j})=y_{j},S(x_{y+1})=y_{j+1}$，得到:
$S(x)=M_{j}\frac{(x_{j+1}-x)^{3}}{6h_{j}}+M_{j+1}\frac{(x-x_{j})^{3}}{6h_{j}}+(y_{j}-\frac{M_{j}h_{j}^{2}}{6})\frac{x_{j+1}-x}{h_{j}}+(y_{j+1}-\frac{M_{j+1}h_{j}^{2}}{6})\frac{x-x_{j}}{h_{j}}$
对上式求导得到：
$S^{‘}(x)=-M_{j}\frac{(x_{j+1}-x)^{2}}{2h_{j}}+M_{j+1}\frac{(x-x_{j})^{2}}{2h_{j}}+\frac{y_{j+1}-y_{j}}{h_{j}}-\frac{M_{j+1}-M_{j}}{6}h_{j}$
根据在$x_{j},1\leq j\leq n-1$处一阶导数相等可得：
$\mu_{j}M_{j-1}+2M_{j}+\lambda _{j}M_{j+1}=d_{j},1\leq j\leq n-1$
其中$\mu_{j}=\frac{h_{j-1}}{h_{j-1}+h_{j}}$,$\lambda _{j}=\frac{h_{j}}{h_{j-1}+h_{j}}$,$d_{j}=6\frac{f[x_{j},x_{j+1}]-f[x_{j-1},x_{j}]}{h_{j-1}+h_{j}}$,$f[x_{j},x_{j+1}]=\frac{y_{j+1}-y_{j}}{x_{j+1}-x{j}}$
同时对于边界条件$S^{‘}(x_{0})=f^{‘}(x_{0})$,$S^{‘}(x_{n})=f^{‘}(x_{n})$可得：$2M_{0}+M_{1}=\frac{6}{h_{0}}(f[x_{0},x_{1}]-f^{‘}(x_{0}))$,$M_{n-1}+2M_{n}=\frac{6}{h_{n-1}}(f^{‘}(x_{n})-f[x_{n-1},x_{n}])$.
#include <cstdio> #include <cstdlib> #include <algorithm> #include <cmath> #include <cassert> #include <ctime> class MclVector { public: int n; double *Mat; /** type=0: 列向量 type=1: 行向量 **/ int type; MclVector() { Mat=NULL; n=0; } MclVector(int len,double initVal=0.0) { n=len; Mat=new double[n+1]; for(int i=0;i<=n;i++) Mat[i]=initVal; type=0; } double operator[](int id) const { return Mat[id]; } double& operator[](int id) { return Mat[id]; } double length() const { double sum=0; for(int i=1;i<=n;i++) sum+=Mat[i]*Mat[i]; return sqrt(sum); } MclVector operator*(double val) const { MclVector ans=MclVector(n); for(int i=1;i<=n;i++) ans[i]=Mat[i]*val; return ans; } MclVector operator/(double val) const { MclVector ans=MclVector(n); for(int i=1;i<=n;i++) ans[i]=Mat[i]/val; return ans; } MclVector operator+(const MclVector &newVector) const { MclVector ans=MclVector(n); for(int i=1;i<=n;i++) ans[i]=Mat[i]+newVector[i]; return ans; } MclVector operator-(const MclVector &newVector) const { MclVector ans=MclVector(n); for(int i=1;i<=n;i++) ans[i]=Mat[i]-newVector[i]; return ans; } MclVector operator*=(double val) { for(int i=1;i<=n;i++) Mat[i]=Mat[i]*val; return *this; } MclVector operator/=(double val) { for(int i=1;i<=n;i++) Mat[i]=Mat[i]/val; return *this; } MclVector operator+=(const MclVector &newVector) { for(int i=1;i<=n;i++) Mat[i]+=newVector[i]; return *this; } MclVector operator-=(const MclVector &newVector) { for(int i=1;i<=n;i++) Mat[i]-=newVector[i]; return *this; } MclVector GetTranspose() const { MclVector ans=*this; ans.type=1; return ans; } void print() const { for(int i=1;i<=n;i++) printf("%8.3lf ",Mat[i]); puts(""); } }; class MclMatrix { public: int row,col; MclVector *Mat; MclMatrix() {Mat=NULL;} MclMatrix(int _row,int _col,double initVal=0.0) { row=_row; col=_col; Mat=new MclVector[row+1]; for(int i=0;i<=row;i++) Mat[i]=MclVector(col,initVal); } void setIdentityMatrix() { for(int i=1;i<=row;i++) { for(int j=1;j<=col;j++) { if(i==j) Mat[i][j]=1; else Mat[i][j]=0; } } } MclMatrix GetTranspose() const { MclMatrix ans=MclMatrix(col,row); for(int i=1;i<=ans.row;i++) { for(int j=1;j<=ans.col;j++) { ans[i][j]=Mat[j][i]; } } return ans; } void print() const { for(int i=1;i<=row;i++) Mat[i].print(); puts(""); } MclVector& operator[](int id) const { return Mat[id]; } MclVector& operator[](int id) { return Mat[id]; } MclMatrix operator*(const MclMatrix &Matrix) const { MclMatrix ans=MclMatrix(row,Matrix.col); for(int i=1;i<=row;i++) { for(int j=1;j<=Matrix.col;j++) { for(int k=1;k<=col;k++) { ans[i][j]+=Mat[i][k]*Matrix[k][j]; } } } return ans; } MclMatrix operator+(const MclMatrix &Matrix) const { MclMatrix ans=MclMatrix(row,Matrix.col); for(int i=1;i<=row;i++) { for(int j=1;j<=Matrix.col;j++) { ans[i][j]=Mat[i][j]+Matrix[i][j]; } } return ans; } MclMatrix operator-(const MclMatrix &Matrix) const { MclMatrix ans=MclMatrix(row,Matrix.col); for(int i=1;i<=row;i++) { for(int j=1;j<=Matrix.col;j++) { ans[i][j]=Mat[i][j]-Matrix[i][j]; } } return ans; } MclVector GetCol(int colId) const { MclVector ans=MclVector(row); for(int i=1;i<=row;i++) ans[i]=Mat[i][colId]; return ans; } MclVector GetRow(int rowId) const { MclVector ans=MclVector(row); for(int i=1;i<=col;i++) ans[i]=Mat[rowId][i]; return ans; } MclMatrix operator*=(const MclMatrix &Matrix) { return *this=*this*Matrix; } MclMatrix operator+=(const MclMatrix &Matrix) { return *this=*this+Matrix; } MclMatrix operator-=(const MclMatrix &Matrix) { return *this=*this-Matrix; } MclMatrix operator*(double x) const { MclMatrix ans=*this; for(int i=1;i<=row;i++) { for(int j=1;j<=col;j++) { ans[i][j]*=x; } } return ans; } }; MclMatrix vectorMulVector(const MclVector &A,const MclVector& B) { if(A.type==0) { MclMatrix ans=MclMatrix(A.n,B.n); for(int i=1;i<=A.n;i++) { for(int j=1;j<=B.n;j++) { ans[i][j]+=A[i]*B[j]; } } return ans; } else { assert(A.n==B.n); MclMatrix ans=MclMatrix(1,1); for(int i=1;i<=A.n;i++) { ans[1][1]+=A[i]*B[i]; } return ans; } } int sgn(double x) { if(x<-0.000001) return -1; if(x>0.000001) return 1; return 0; } /** 高斯消去 A为n*n B为n*1 **/ MclVector Gauss(MclMatrix A,MclVector B) { const int n=A.row; const int m=A.col; assert(m==B.n&&m==n); MclVector ans=MclVector(m); for(int i=1;i<=n;i++) { int row=i; for(int j=i;j<=n;j++) { if(sgn(A[j][i])) { row=i; break; } } if(row!=i) { for(int j=1;j<=n;j++) std::swap(A[i][j],A[row][j]); std::swap(B[i],B[row]); } for(int j=1;j<=n;j++) if(i!=j) { const double det=A[j][i]/A[i][i]; for(int k=i;k<=n;k++) A[j][k]-=det*A[i][k]; B[j]-=det*B[i]; } } for(int i=1;i<=n;i++) ans[i]=B[i]/A[i][i]; return ans; } /** n: 点个数[0,n] RangeSplitPos: 分界点 大小(n+1) SplitPosValue: 函数值 大小(n+1) LValue,RValue: 两侧的一阶导数 返回值： 每个区间的多项式 大小n **/ std::vector<Polynomial> CubicSplineInterpolation( const int n, const std::vector<double> RangeSplitPos, const std::vector<double> SplitPosValue, const double LValue, const double RValue) { #define x(i) RangeSplitPos[i] #define y(i) SplitPosValue[i] #define h(i) (RangeSplitPos[i+1]-RangeSplitPos[i]) #define f(i,j) ((y(j)-y(i))/h(i)) double *mou=new double[n+1]; double *lamd=new double[n+1]; double *d=new double[n+1]; memset(mou,0,sizeof(double)*(n+1)); memset(lamd,0,sizeof(double)*(n+1)); memset(d,0,sizeof(double)*(n+1)); lamd[0]=1; d[0]=6.0/h(0)*(f(0,1)-LValue); mou[n]=1; d[n]=6.0/h(n-1)*(RValue-f(n-1,n)); for(int j=1;j<=n-1;j++) { mou[j]=h(j-1)/(h(j-1)+h(j)); lamd[j]=h(j)/(h(j-1)+h(j)); d[j]=6.0/(h(j-1)+h(j))*(f(j,j+1)-f(j-1,j)); } MclMatrix A=MclMatrix(n+1,n+1); for(int i=1;i<=n+1;i++) { if(i>1) A[i][i-1]=mou[i-1]; if(i<n+1) A[i][i+1]=lamd[i-1]; A[i][i]=2; } MclVector B=MclVector(n+1); for(int i=1;i<=n+1;i++) B[i]=d[i-1]; delete[] mou; delete[] lamd; delete[] d; MclVector MArr=Gauss(A,B); std::vector<Polynomial> LastAns; for(int j=0;j<n;j++) { const double Mj=MArr[j+1]; const double Mj1=MArr[j+2]; Polynomial poly; double tmp=-Mj/(6.0*h(j)); poly.add(3,tmp); poly.add(2,-3*x(j+1)*tmp); poly.add(1,3*x(j+1)*x(j+1)*tmp); poly.add(0,-x(j+1)*x(j+1)*x(j+1)*tmp); tmp=Mj1/(6*h(j)); poly.add(3,tmp); poly.add(2,-3*x(j)*tmp); poly.add(1,3*x(j)*x(j)*tmp); poly.add(0,-x(j)*x(j)*x(j)*tmp); tmp=(y(j)-Mj*h(j)*h(j)/6.0)/h(j); poly.add(0,tmp*x(j+1)); poly.add(1,-tmp); tmp=(y(j+1)-Mj1*h(j)*h(j)/6.0)/h(j); poly.add(0,-tmp*x(j)); poly.add(1,tmp); LastAns.push_back(poly); } return LastAns; }
以上是关于三次样条插值的主要内容，如果未能解决你的问题，请参考以下文章