シンプソン公式による数値積分

応用編

内容

シンプソン公式を使った数値積分です。 台形公式による数値積分と内容はほとんど一緒です。 台形公式とシンプソン公式との間の関係式を用いています。

ソース

ヘッダファイル

SimpsonQuad.h

#ifndef _SimpsonQuad_h
#define _SimpsonQuad_h

#include <iostream>
#include <string>
#include <octave/config.h>
//#include <octave/oct-cmplx.h>
//#include <octave/CColVector.h>
//#include <octave/CMatrix.h>

using namespace std;

//**************************************************
// Templates
//**************************************************

/**
 * A template SimpsonQuad<K,KIntegFunc>.
 * Quadrature, automatic integration by trapezoidal rule.
 * K          : K + K, double * K , K = K must be defined.
 * KIntegFunc : K KIntegFunc::eval(double t) must be defined. 
 */

template <class K, class KIntegFunc> 
class SimpsonQuad { 
  int step_min; /**< n_slice = 2^{step_min} to 2^{step_max} */
  int step_max; /**< n_slice = 2^{step_min} to 2^{step_max} */
  double eps; /**< error : |quad_new-quad_old| < 3.0 * eps */
  string calcLog; /**< log */
  /**
   * This returns |q_new-q_old|.
   * Pure virtual function. This should be implemented in subclass 
   */
  virtual double abserr(const K& q_new, const K& q_old) = 0;
  /**
   * This sets ojb -> 0.
   * Pure virtual function. This should be implemented in subclass 
   */
  virtual void initToZero(K& obj) = 0;
 public:
  /**
   * Constructor.
   * @param eps error
   * @param step_min  n_slice = 2^{min} to 2^{max}
   * @param step_max  n_slice = 2^{min} to 2^{max}
   */
  SimpsonQuad(double eps, int step_min, int step_max):
    eps(eps), step_min(step_min), step_max(step_max), 
    calcLog("object generated.\n") {}
  /**
   * Constructor
   */
  SimpsonQuad() :
    eps(1.e-5), step_min(4), step_max(12),
    calcLog("object generated.\n") {} 
  /**
   * setter of parameters (eps, step_min, step_max)
   */
  void setParameters(double eps = 1.e-5, int step_min = 4, int step_max = 12) {
    this->eps = eps;
    this->step_min = step_min;
    this->step_max = step_max;
  }
  /**
   * integration.
   * @return int_a^b f(x) dx
   */
  K do_integrate(KIntegFunc* f, double a, double b);
  string getCalcLog() { return calcLog; }
  virtual ~SimpsonQuad() {}
};

//==========
template <class K, class KIntegFunc>
K SimpsonQuad<K, KIntegFunc>::do_integrate(KIntegFunc* f,double a, double b) {
  double h; // slice
  K q_new, q_old; // Trapesoidal
  K s_new, s_old; // Simpson
  K tmp;
  int n_slice;// n_slice = 2^{step} 
  int step; 
  
  // n_slice = 2^0 = 1
  h = b - a;
  s_new = q_new = tmp = 0.5 * h * (f->eval(a) + f->eval(b));
  step = 0;
  n_slice = 1;
  
  // automatic integration : h -> h/2
  do{
    if (step > step_max) {
      calcLog += "Simpson rule does not converge.\n";
      break;
    }
    step++;
    n_slice = 2 * n_slice;    
    h = h / 2.0;
    q_old = q_new;
    s_old = s_new;
    initToZero(tmp); // tmp->0
    for (int i = 1; i < n_slice; i = i + 2) { // i = 1, 3, 5, ... n_slice-1
      tmp = tmp + f->eval(a+i*h);
    }    
    q_new = 0.5 * q_old + h * tmp;
    calcLog += "Evaluation step\n";
    s_new = (1.0/3.0) * (4.0*q_new - q_old); // Trapezoidal to Simpson
  } while ((abserr(s_new, s_old) > 1.0 * eps) || step < step_min);
  
  return s_new;
}

#endif

ComplexSimpsonQuad.h, DoubleSimpsonQuad.h,ComplexColumnVectorSimpsonQuad.h, ComplexMatrixSimpsonQuad.h

サンプルプログラム

sample.cpp -- サンプルプログラム

#include "DoubleSimpsonQuad.h"
#include "ComplexSimpsonQuad.h"
#include "ComplexColumnVectorSimpsonQuad.h"
#include "ComplexMatrixSimpsonQuad.h"

using namespace std;

//=====
class MyDoubleIntegFunc { 
public:
  double eval(double t) { return exp(-t*t); };
};

//=====
class MyComplexIntegFunc {
public:
 Complex eval(double t) {
    Complex const IU(0.0, 1.0);
    return exp(-(1.0+IU)*t);
  }
};

//=====
class MyComplexColumnVectorIntegFunc {
public:
  ComplexColumnVector eval(double t) {
    Complex const IU(0.0, 1.0);
    ComplexColumnVector retval(3);
    retval(0) = exp(-t*t);
    retval(1) = exp(-(1.0+IU)*t);
    retval(2) = 1.0+IU;
    return retval;
  }
};


//=====
class MyComplexMatrixIntegFunc {
public:
  ComplexMatrix eval(double t) {
    Complex const IU(0.0, 1.0);
    ComplexMatrix retval(2,2);
    retval(0,0) = exp(-t*t)       ; retval(0,1) = 1.0+IU;
    retval(1,0) = exp(-(1.0+IU)*t); retval(1,1) = t*t;
    return retval;
  }
};

int main()
{ 
  DoubleSimpsonQuad<MyDoubleIntegFunc> dQuad;  
  MyDoubleIntegFunc myDIntegFunc;
  cout << "<< DoubleSimpsonQuad >>" << endl;
  cout << dQuad.do_integrate(&myDIntegFunc, 0.0, 10.0) << endl;  
  cout << dQuad.getCalcLog() << endl;
  
  ComplexSimpsonQuad<MyComplexIntegFunc> cQuad(1.e-7);
  MyComplexIntegFunc myCIntegFunc;
  cout << "<< ComplexSimpsonQuad >>" << endl;
  cout << cQuad.do_integrate(&myCIntegFunc, 0.0, 10.0) << endl;
  //cout << cQuad.getCalcLog() << endl;

  ComplexColumnVectorSimpsonQuad<MyComplexColumnVectorIntegFunc> cColVecQuad;
  cColVecQuad.setParameters(1.e-7, 4, 12);
  MyComplexColumnVectorIntegFunc myCColVecIntegFunc;  
  cout << "<< ComplexColumnVectorSimpsonQuad >>" << endl;
  cout << cColVecQuad.do_integrate(&myCColVecIntegFunc, 0.0, 10.0) << endl;
  //cout << cColVecQuad.getCalcLog() << endl; 

  ComplexMatrixSimpsonQuad<MyComplexMatrixIntegFunc> cMatrixQuad;
  MyComplexMatrixIntegFunc myCMatrixIntegFunc;  
  cout << "<< ComplexMatrixSimpsonQuad >>" << endl; 
  cout << cMatrixQuad.do_integrate(&myCMatrixIntegFunc, 0.0, 10.0) << endl;
  //cout << cMatrixQuad.getCalcLog() << endl;

  return 0;
}

出力

<< DoubleSimpsonQuad >>
0.886227
object generated.
Evaluation step
Evaluation step
Evaluation step
Evaluation step
Evaluation step
Evaluation step

<< ComplexSimpsonQuad >>
(0.500007,-0.500031)
<< ComplexColumnVectorSimpsonQuad >>
(0.886227,0)
(0.500007,-0.500031)
(10,10)

<< ComplexMatrixSimpsonQuad >>
 (0.886227,0) (10,10)
 (0.500006,-0.500031) (333.333,0) 

コメント