被積分関数が複素数やベクトル、行列の場合での数値積分法が見当たらなかったので、 台形公式での数値積分のプログラムを作ってみました。 共立出版株式会社 数値計算法[第2版] 小澤一文著、第6章、6.3 台形公式の自動積分法を参考にしました。 複素数Complex、ComplexColumnVector?、CompexMatrix? クラスの積分ができます。 計算の精度としてはデフォルトで小数点以下5桁程度出るはずです。 プログラムにはまだバグがたくさん残っていると思いますので、 万が一利用しようとする場合は十分気をつけてください。
TrapezoidalQuad.h -- 台形公式のアルゴリズムを記述してあります。
#ifndef _TrapezoidalQuad_h #define _TrapezoidalQuad_h #include <iostream> #include <string> #include <octave/config.h> using namespace std; //************************************************** // Templates //************************************************** /** * A template SimpsonQuad<K,KFunc>. * 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 TrapezoidalQuad { 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} */ TrapezoidalQuad(double eps, int step_min, int step_max): eps(eps), step_min(step_min), step_max(step_max), calcLog("object generated.\n") {} /** * Constructor */ TrapezoidalQuad() : 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 ~TrapezoidalQuad() {} }; //========== template <class K, class KIntegFunc> K TrapezoidalQuad<K, KIntegFunc>::do_integrate(KIntegFunc* f,double a, double b) { double h; // slice K q_new, q_old; // Trapesoidal K tmp; int n_slice;// n_slice = 2^{step} int step; // n_slice = 2^0 = 1 h = b - a; 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 += "Trapezoidal rule does not converge.\n"; break; } step++; n_slice = 2 * n_slice; h = h / 2.0; q_old = q_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"; } while ((abserr(q_new, q_old) > 1.0 * eps) || step < step_min); return q_new; } #endif
DoubleTrapezoidalQuad.h,ComplexTrapezoidalQuad.h, ComplexColumnVectorTrapezoidalQuad.h, ComplexMatrixTrapezoidalQuad.h -- TrapezoidalQuadをdouble型、Complex型 ComplexColumnVector型、ComplexMatrix型で使えるように展開してます。
sample.cpp -- サンプルプログラム
#include "DoubleTrapezoidalQuad.h" #include "ComplexTrapezoidalQuad.h" #include "ComplexColumnVectorTrapezoidalQuad.h" #include "ComplexMatrixTrapezoidalQuad.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() { DoubleTrapezoidalQuad<MyDoubleIntegFunc> dQuad; MyDoubleIntegFunc myDIntegFunc; cout << "<< DoubleTrapezoidalQuad >>" << endl; cout << dQuad.do_integrate(&myDIntegFunc, 0.0, 10.0) << endl; cout << dQuad.getCalcLog() << endl; ComplexTrapezoidalQuad<MyComplexIntegFunc> cQuad(1.e-7); MyComplexIntegFunc myCIntegFunc; cout << "<< ComplexTrapezoidalQuad >>" << endl; cout << cQuad.do_integrate(&myCIntegFunc, 0.0, 10.0) << endl; //cout << cQuad.getCalcLog() << endl; ComplexColumnVectorTrapezoidalQuad<MyComplexColumnVectorIntegFunc> cColVecQuad; cColVecQuad.setParameters(1.e-7, 4, 12); MyComplexColumnVectorIntegFunc myCColVecIntegFunc; cout << "<< ComplexColumnVectorTrapezoidalQuad >>" << endl; cout << cColVecQuad.do_integrate(&myCColVecIntegFunc, 0.0, 10.0) << endl; //cout << cColVecQuad.getCalcLog() << endl; ComplexMatrixTrapezoidalQuad<MyComplexMatrixIntegFunc> cMatrixQuad; MyComplexMatrixIntegFunc myCMatrixIntegFunc; cout << "<< ComplexMatrixTrapezoidalQuad >>" << endl; cout << cMatrixQuad.do_integrate(&myCMatrixIntegFunc, 0.0, 10.0) << endl; //cout << cMatrixQuad.getCalcLog() << endl; return 0; }
<< DoubleTrapezoidalQuad >> 0.886227 object generated. Evaluation step Evaluation step Evaluation step Evaluation step Evaluation step << ComplexTrapezoidalQuad >> (0.500007,-0.500031) << ComplexColumnVectorTrapezoidalQuad >> (0.886227,0) (0.500007,-0.500031) (10,10) << ComplexMatrixTrapezoidalQuad >> (0.886227,0) (10,10) (0.500007,-0.500031) (333.333,0)
あまり使う人はいないと思いますが、とりあえず自分用のメモとして。
class MyIntegFunc { public: ComplexMatrix eval(double x) { ComplexMatrix retval; /* hogehoge */ return retval; } };
ComplexMatrixTrapezoidalQuad<MyIntegFunc> quad;
MyIntefFunc integrand;
ComplexMatrix answer = quad.do_integrand(&integrand, 0, 10);
サンプルプログラムのように、ベクトルや行列の各成分の形がわかっている場合には、 各成分ごとにDefQuadクラスを利用した方が良いでしょう。→数値積分を行いたい
&ref(): File not found: "wikilogo1.png" at page "台形公式による数値積分"; |