Gauss legendre法とは

高精度、高速の数値積分法。関数の補間にルジャンドル関数を用いて積分する(らしい)。
ニュートンコーツ法が等間隔の離散化を用いて誤差がその微小区間によるのに対しGauss〜法は分割数に依存してnが大きくなるほど急速に良くなる。
そのためよく使われる使い方はnとn+1の時を比べてそれが認められる誤差以下なら計算を打ち切ると言った使い方をするようだ。

legendre関数(多項式)

legendre関数$P_n(x)$は最大次$x^n$の$-1\sim1$で直交性を持つ多項式。この直交性があるから良い関数の補間が得られる。その他の数学的には球対象な系、$r,\theta,\phi$に分離できる系の微分方程式を解くとほぼ必ず出てくる。

SpecialFunction.hh

#pragma once

#include "common.hh"

namespace MyMath{
  namespace SpecialFunction{
    std::function<double(double)> legendre(int n);
    using complex=std::complex<double>;
    std::function<complex(complex)> legendreComplex(int n);
  };
};

SpecialFunction.cc

#include "SpecialFunction.hh"
#include "Util.hh"

using namespace MyMath;
//factrialは別途実装 n!
std::function<double(double)> SpecialFunction::legendre(int n){
  int max= n%2==0 ? n/2 : (n-1)/2;
  std::vector<double> cofficients(max+1);
  for( int k=0; k<=max; k++ ){
    cofficients[k]=factrial(2*n-2*k)/(factrial(k)*factrial(n-k)*factrial(n-2*k)*pow(2, n));
    if( k%2!=0 ) cofficients[k]*=-1;
  }

  return [cofficients, n](double x){
    double val=0;
    for( size_t k=0; k<cofficients.size(); k++ ) val+=cofficients[k]*pow(x, n-2*k);
    return val;
  };
}

using complex=std::complex<double>;
std::function<complex(complex)> SpecialFunction::legendreComplex(int n){
  int max= n%2==0 ? n/2 : (n-1)/2;
  std::vector<double> cofficients(max);
  for( int k=0; k<max; k++ ){
    cofficients[k]=factrial(2*n-2*k)/(factrial(k)*factrial(n-k)*factrial(n-2*k)*pow(2, n));
    if( k%2!=0 ) cofficients[k]*=-1;
  }

  return [cofficients, n](complex x){
    complex val=0;
    for( size_t k=0; k<cofficients.size(); k++ ) val+=cofficients[k]*pow(x, n-2*k);
    return val;
  };
}

legendre関数は多項式でその係数は漸化式で計算できるのでパラメーターを計算してstd::functionで返す。階乗はベタ書きしても良かったが可読性が下がるので関数化した(ベタ書きしてバグらせた)

Gauss legendre法のパラメーター計算

実装は以前紹介したIntegrator(積分器)に書く計算自体はメモリに保存し複数回しないようにしたがユーザーに公開する必要はないのでソースファイルに書く。
パラメーターはlegendre関数の0点とその時の重さである。
legendre関数の0点は$-1\sim1$の間ラフに調査し、符号が変化した時点でnewton法で探している。
$n\sim25$ぐらいでnewton法が収束しなくなった...、エラー設定が厳しいせいだと思われる。(そこまで使うことは多分ないだろうが)

Integrator.cc

#include "Integrator.hh"

using namespace MyMath;

namespace{
  std::vector<std::vector<double>> legendre_0points(1, std::vector<double>(1, 1.0));
  std::vector<std::vector<double>> gaussLegendre_weight(1, std::vector<double>(1, 1.0));

  void searchGaussLegendrePoint(int N){
    std::cout<<"===== gauss Legendre method point search ====="<<std::endl;
    std::cout<<"Searched points : "<<legendre_0points.size()<<std::endl;
    for( int n=static_cast<int>(legendre_0points.size()); n<=N; n++ ){
      std::function<double(double)> legendre=SpecialFunction::legendre(n);
      std::vector<double> zero_point;
      double x0=-1.0;
      while( x0<1.0 ){
        double x1=x0+0.2/n;
        if( legendre(x0)*legendre(x1)<0 ) zero_point.push_back(Solver::newton(legendre, x1));
        x0=x1;
      }

      std::vector<double> weight(zero_point.size());
      for( size_t i=0; i<zero_point.size(); i++ ){
        weight[i]=2.0*(1.0-zero_point[i]*zero_point[i])/pow(n*SpecialFunction::legendre(n-1)(zero_point[i]), 2);
      }
      legendre_0points.push_back(zero_point);
      gaussLegendre_weight.push_back(weight);
    }
  }
};

一応パラメーターは5次まで目視で確認した。なんかイマイチな実装なきがする探せばネットにパラメーターもあるだろうしベタ書きでもいいかなと思いだした。

Integrator.cc

double Integrator::gaussLegendre(const std::function<double(double)> &f, double x0, double x1, int n)
{
  if( x0==x1 ) return 0;
  double c1=0.5*(x1-x0);
  double c0=0.5*(x1+x0);

  if( static_cast<int>(gaussLegendre_weight.size())<=n ) throw std::runtime_error("Integrator::gaussLegendre  exceed max size");

  double val=0;
  for( int i=0; i<n; i++ ) val+=f(c1*legendre_0points[n][i]+c0)*gaussLegendre_weight[n][i];

  return c1*val;
}

double Integrator::gaussLegendre(const std::function<double(double)> &f, double x0, double x1, double error){
  searchGaussLegendrePoint(25);
  double init=gaussLegendre(f, x0, x1, 0);
  double next=gaussLegendre(f, x0, x1, 1);
  int n=2;
  while( fabs(init-next)>error ){
    init=next;
    next=gaussLegendre(f, x0, x1, n);
    n++;
  }
  return next;
}

計算部分(本体)の実装、変数変換$t=\frac{a-b}{2}x+\frac{a+b}{2}$で積分範囲を$[a:b]\rightarrow[-1:1]$に変更している。一つ実装で注意する点はx0==x1で実数の完全な一致判定をしています。その場合は積分は必ず0です。

確認は前のニュートンコーツ法用のものに追加で出力して確認しました。

参考にしたサイトなど

Wikipedia:ガウス求積
 “ホームレスこ～よ～”:ガウスルジャンドル積分---パラメーターのexcelファイルあり