c++14で数値計算(2)---数値積分ニュートンコーツ法

# ニュートンコーツ法で数値積分したときのメモ
合成シンプソン則を使って積分しているその誤差は$\sim dx^4$程度である。
最初は合成シンプソン則だけを使おかと思ったが台形則を使う必要が出てきたので実装した。

```c++:Integrator.hh
#pragma once

#include "common.hh"
#include "SpecialFunction.hh"
#include "Solver.hh"

namespace MyMath{
  namespace Integrator{
    double trapezoidal(const std::function<double(double)>&, double x0, double x1, int N=1000);
    double simpson(const std::function<double(double)>&, double x0, double x1, int N=1000);

    using complex=std::complex<double>;
    complex trapezoidal(const std::function<complex(complex)>&, const complex &x0, const complex &x1, int N=1000);
    complex simpson(const std::function<complex(complex)>&, const complex &x0, const complex &x1, int N=1000);
  };
};
```
`std::complex<double>`と常に打つと長くなるのでusingを使った。

```c++:Integrator.cc
#include "Integrator.hh"

using namespace MyMath;

double Integrator::trapezoidal(const std::function<double(double)> &f, double x0, double x1, int N){
  const double delta=(x1-x0)/N;
  double val=0.5*(f(x0)+f(x1));
  for( int i=1; i<N; ++i ) val+=f(x0+i*delta);
  return delta*val;
}

double Integrator::simpson(const std::function<double(double)> &f, double x0, double x1, int N){
  const double delta=(x1-x0)/(2.0*N);
  double val=f(x0)+f(x1);
  for( int i=1; i<2*N; i+=2 ) val+=4.0*f(x0+i*delta);
  for( int i=2; i<2*N; i+=2 ) val+=2.0*f(x0+i*delta);
  return delta*val/3.0;
}

using complex=std::complex<double>;
complex Integrator::trapezoidal(const std::function<complex(complex)> &f, const complex &x0, const complex &x1, int N){
  const complex delta=(x1-x0)/static_cast<double>(N);
  complex val=0.5*(f(x0)+f(x1));
  for( int i=1; i<N; ++i ) val+=f(x0+static_cast<double>(i)*delta);
  return delta*val;
}

complex Integrator::simpson(const std::function<complex(complex)> &f, const complex &x0, const complex &x1, int N){
  const complex delta=(x1-x0)/(2.0*N);
  complex val=f(x0)+f(x1);
  for( int i=1; i<2*N; i+=2 ) val+=4.0*f(x0+static_cast<double>(i)*delta);
  for( int i=2; i<2*N; i+=2 ) val+=2.0*f(x0+static_cast<double>(i)*delta);
  return delta*val/3.0;
}
```
`using namespace std;`をすると自分定義のcomplexとstd::complexがぶつかった。ついつい反射的に`using namespace std;`するが当然ながら必要ないのに名前空間を汚染するのは悪である。
int型とcomplex<double>は演算が定義できていないのでdouble型にキャストしている。

## チェック用main関数
```c++:integral.cc
#include "common.hh"
#include "Integrator.hh"

using namespace std;
using namespace MyMath;

int main()
{
  cout<<"===== Integral test program ====="<<endl;
  cout<<" x  : sin(x)    : S(cos(x)dx "<<endl;
  for( double x=0; x<2*MyMath::PI; x+=0.01 ){
    cout<<setw(5)<<x<<"  "<<setw(12)<<sin(x)<<"  "<<setw(12)<<Integrator::simpson([](double x){ return cos(x); }, 0, x)<<endl;
  }
  auto func=[](std::complex<double> x){ return 1.0/x; };
  complex<double> val=Integrator::simpson(func, complex<double>(-1, -1), complex<double>(1, -1));
  val+=Integrator::simpson(func, complex<double>(1, -1), complex<double>(1, 1));
  val+=Integrator::simpson(func, complex<double>(1,  1), complex<double>(-1, 1));
  val+=Integrator::simpson(func, complex<double>(-1, 1), complex<double>(-1,-1));
  cout<<2*PI<<endl;
  cout<<"1/z residu val = "<<val<<endl;
  return 0;
}
```
$\int \cos(x)dx$と$sin(x)$の値を比べている。
複素数は$1/z$の極の周りを一周する積分をしている。これは留数定理より$2\pi i$になる。
こんなに簡単に複素積分ができるのはちょっと感動である。

## 最後に
複素数の積分は意外と解説しているサイトなどが少ないです。Mathmaticaあたりの数学ライブラリは当たり前のようにサポートしています。
誤差な
誤差$d\cos/dx\sim 1-dx^2$の誤差があるので字数がひとつ下がります。
# 参考にしたサイト
[Wikipedia:ニュートンコーツ法](https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E3%83%BB%E3%82%B3%E3%83%BC%E3%83%84%E3%81%AE%E5%85%AC%E5%BC%8F)
[複素関数の積分に対する数値積分の試み](http://www.it-hiroshima.ac.jp/institution/library/pdf/library_research44_26.pdf)