阿克玛插值(2D)

最新推荐文章于 2024-08-11 14:33:36 发布

yliang

最新推荐文章于 2024-08-11 14:33:36 发布

阅读量2k

点赞数

文章标签： iterator less exception distance insert header

本文链接：https://blog.csdn.net/yliang/article/details/236840

版权

该博客介绍了如何实现二维阿克玛插值，包括插入点、获取插值值和导数值的方法。使用了`Point2D`类来存储坐标，并通过`InterpolationSmooth`类进行插值计算，利用`lower_bound`查找指定位置的插值点，并应用多项式插值公式。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

#ifndef INTERPOLATION_HEADER

#define INTERPOLATION_HEADER

#include<cmath>

#include<set>

#include<vector>

#include<algorithm>

using namespace std;

#define DELTA 1e-6

#define VERY_CLOSE(x,y) fabs((x)-(y))<DELTA

template<typename T>

class Point2D

{

public:

Point2D(T _X=0,T _Y=0)

:X(_X),Y(_Y)

{

}

T X,Y;

};

bool MyFunc(const Point2D<double>& P1, const Point2D<double>& P2)

{

return P1.X+DELTA<P2.X;

}

class Point2D_X_Less_Delta

{

public:

template<typename T>

bool operator()(const Point2D<T>& P1, const Point2D<T>& P2)

{

return P1.X+DELTA<P2.X;

}

};

template<typename T>

ostream & operator <<(ostream &os,const Point2D<T>& P)

{

return os<<P.X<<"/t"<<P.Y<<endl;

}

template<typename T>

istream & operator >>(ostream &is,Point2D<T>& P)

{

return is>>P.X>>P.Y;

}

class InterpolationSmooth

{

public:

InterpolationSmooth(){}

template<typename InputIterator>

InterpolationSmooth(InputIterator iter1,InputIterator iter2):

Pts(iter1,iter2)

{

}

void Insert(const Point2D<double> &P)

{

Pts.insert(P);

}

template<typename InputIterator>

void Insert(InputIterator iter1,InputIterator iter2)

{

Pts.insert(iter1,iter2);

}

double GetValue(double X)

{

if (X<(*Pts.begin()).X)

{

throw exception("Iuput is too small!");

}

set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter=lower_bound(Pts.begin(),Pts.end(),Point2D<double>(X,0),MyFunc);

if (iter==Pts.begin())

{

iter++;

}

int k=distance(Pts.begin(),iter);

k--;

double Xkp1=(*iter).X;

iter--;

double s0=(*iter).Y;

double s1=G(k);

double s2=(U(k)*3-G(k)*2-G(k+1))/(Xkp1-(*iter).X);

double s3=(G(k)+G(k+1)-U(k)*2.0)/((Xkp1-(*iter).X)*(Xkp1-(*iter).X));

s1*=X-(*iter).X;

s2*=(X-(*iter).X)*(X-(*iter).X);

s3*=(X-(*iter).X)*(X-(*iter).X)*(X-(*iter).X);

return s0+s1+s2+s3;

}

double GetDerivativeValue(double X)

{

if (X<(*Pts.begin()).X)

{

throw exception("Iuput is too small!");

}

set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter=lower_bound(Pts.begin(),Pts.end(),Point2D<double>(X,0),MyFunc);

if (iter==Pts.begin())

{

iter++;

}

int k=distance(Pts.begin(),iter);

k--;

double Xkp1=(*iter).X;

iter--;

double s1=G(k);

double s2=(U(k)*3-G(k)*2-G(k+1))/(Xkp1-(*iter).X);

double s3=(G(k)+G(k+1)-U(k)*2.0)/((Xkp1-(*iter).X)*(Xkp1-(*iter).X));

s2*=(X-(*iter).X)*2.0;

s3*=(X-(*iter).X)*(X-(*iter).X)*3.0;

return s1+s2+s3;

}

set<Point2D<double>,Point2D_X_Less_Delta >::iterator begin()

{

return Pts.begin();

}

set<Point2D<double>,Point2D_X_Less_Delta >::iterator end()

{

return Pts.end();

}

size_t size()

{

return Pts.size();

}

private:

double U(int _Pos)

{

int Num=(int)(Pts.size());

if(Num<=2 ||_Pos>Num) throw exception("size is too small!");

if (_Pos>=Num-1)

{

return U(_Pos-1)*2.0-U(_Pos-2);

}

if (_Pos<0)

{

return U(_Pos+1)*2.0-U(_Pos+2);

}

set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter1=Pts.begin();

set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter2=Pts.begin();

advance(iter1,_Pos);

advance(iter2,_Pos+1);

return ((*iter2).Y-(*iter1).Y)/((*iter2).X-(*iter1).X);

}

double G(int _Pos)

{

double Uk =U(_Pos),Ukp1=U(_Pos+1),Ukm1=U(_Pos-1),Ukm2=U(_Pos-2);

if ( VERY_CLOSE( Uk ,Ukp1) && VERY_CLOSE(Ukm1,Ukm2) )

{

return ( Uk +Ukm1)/2.0;

}

return (fabs(Ukp1-Uk)*Ukm1+fabs(Ukm1-Ukm2)* Uk )/(fabs(Ukp1-Uk)+fabs(Ukm1-Ukm2));

}

set<Point2D<double>,Point2D_X_Less_Delta > Pts;

};

#endif