多項式近似

多項式近似

多項式近似を行いましょう．

1D有限区間(GSL)

GSLでは, Cubic-splineが簡単に利用できる. データ点のx座標は単調増大である必要がある.

#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
...
class CSpline {
  std::vector& X;
  gsl_interp_accel *acc=nullptr;
  gsl_spline *spline=nullptr;
public:
  CSpline(std::vector& _X):X(_X){
    acc= gsl_interp_accel_alloc();
    spline = gsl_spline_alloc(gsl_interp_cspline, X.size());
  }
  ~CSpline(){
    gsl_spline_free (spline);spline=nullptr;
    gsl_interp_accel_free (acc);acc=nullptr;
  }
  void Set(std::vector& Y){
    if (Y.size()!=X.size())
      throw(std::runtime_error("Your fucking data does not match X-size"));
    gsl_spline_init(spline, X.data(), Y.data(), X.size());
  }
  double Get(double X){return gsl_spline_eval (spline,X, acc);
}

これで, 

  std::vector<double> myX(5),myY(5);
myX[0]=....; myY[0]=...  // データを設定
CSpline CS(myX);         //Xデータを設定
CS.Set(myY);             //Yデータを交換
auto y=CS.Get(1.2);      //補間を実行

こんな感じです．結果は

　

　．．．まあ5点でsinを1周期は最後がダメだよな．知ってた

なお, GSLは, データ範囲外の X の値を与えると「おまえ，補間，つったろ？つったよな？まじかお前，範囲外だよ死ねよ今すぐコラ」とか泣き叫んで働こうとしない．

1D有限区間(BOOST)

BOOSTも頑張ってますよ.

Barycentric_Rational

こいつは, 不均等な格子で与えられたデータに有効です．均等格子では，ただの遅くて精度が悪い近似です．

#include <boost/math/interpolators/barycentric_rational.hpp>
...
std::vector<double> x,y;
...
boost::math::barycentric_rational<double> CS(x.data(),y.data(),x.size());
auto y=CS(1.2)

こんな感じですね. 

BOOSTの奴は少々とぼけており，範囲外の X を与えても，「まあ・・・こんな感じですかね」的な値を返す（下図）

しかしまあなんだな. barycentric_rationalでは, 毎回コンストラクトするのはまあ許しても, std::vectorを毎回作り直していたら，時間がかかって仕方がないので，高速にはならない．と思っていたが, Boost1.70以降で仕様変更されていた.

今まで通り

#include <boost/math/interpolators/barycentric_rational.hpp>
...
boost::math::barycentric_rational<double> CS(x.data(),y.data(),x.size());
auto y=CS(1.2)

で呼び出すと，配列xとyのコピーを作成して利用する．ので，配列の作成に時間がかかるわけだ. そこで

boost::math::barycentric_rational<double> CS(std::move(x),std::move(y));
std::cout << x.size() << std::endl;

と呼び出すのだ！ するとx, yの中身がCSに食われて失われる！つまり, x, y配列サイズは0になってしまうのだ．それでも一向に困らない，という人もいるだろうから，それで配列作成の時間を省略できる.

だが, 「おいらは，これを繰り返し使うんだよ．せっかく確保した配列エリアが，CSと同時に滅亡するのを座視するわけにはいかないんだ！」という人も多いでしょう（普通そうだと思う）．それは

x=CS.return_x();
y=CS.return_y();

で取り返すことが可能です(もちろん内臓を取り出されてしまったCSは骸となり動作不能に陥ります)．返却されるのは配列領域だけで, 値は破壊し尽くされています. 配列を再度ご利用の前には再設定が必要です.

Cardinal cubic B-spline

均等格子に対して極めて高速で精度が高い補間を実行します.

これは

Cardinal quadratic B-spline

Cubicでは嫌な場合があるらしいので, これが活躍する局面も存在するらしいです.

Cardinal quintic B-spline

CubicでもQuadraticでもダメな時にはこれらしいです．

Whittaker-Shannon interpolator

無限回微分可能な近似関数を, 極めて低速に構成する方法です.

Catmull-Rom splines

CGでは美しいらしいです．

Bézier polynomials

CGでは美しいらしいです．

cubic Hermite interpolation

微分方程式の基本解を用いてみます.

modified Akima interpolation

Cubic Hermiteを何か改良したらしい．

PCHIP interpolation

単調増加・単調減少の関数に限定で有効です．

Bilinear Uniform Interpolation

双線形の関数では，何かいいことがあることがあるそうです．

2D有限領域

あ

球面

あ