シンプソン公式を使った数値積分です。 台形公式による数値積分と内容はほとんど一緒です。 台形公式とシンプソン公式との間の関係式を用いています。
#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)
&ref(): File not found: "wikilogo1.png" at page "シンプソン公式による数値積分"; |