内容
被積分関数が複素数やベクトル、行列の場合での数値積分法が見当たらなかったので、 台形公式での数値積分のプログラムを作ってみました。 共立出版株式会社 数値計算法[第２版] 小澤一文著、第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)
使用法
あまり使う人はいないと思いますが、とりあえず自分用のメモとして。
- 被積分関数がComplexMatrix型の場合(その他の場合は適宜読み替えてください)、
被積分関数ComplexMatrix eval(double x)をメンバ関数として持つ、
MyIntegFuncクラスを定義します。
class MyIntegFunc { public: ComplexMatrix eval(double x) { ComplexMatrix retval; /* hogehoge */ return retval; } }; - ComplexMatrixTrapezoidalQuad<MyIntegFunc>型のインスタンスquadを生成します。
ComplexMatrixTrapezoidalQuad<MyIntegFunc> quad;
- MyIntegFunc型の被積分関数用のインスタンスintegrandを生成します。
MyIntefFunc integrand;
- quadのComplexMatrix do_integrate(MyIntegFunc *integrand, double a, double b)メソッドで積分を行います。
ComplexMatrix answer = quad.do_integrand(&integrand, 0, 10);
サンプルプログラムのように、ベクトルや行列の各成分の形がわかっている場合には、 各成分ごとにDefQuadクラスを利用した方が良いでしょう。→数値積分を行いたい
コメント
| &ref(): File not found: "wikilogo1.png" at page "台形公式による数値積分"; |
- FrontPage
- Tips
- 各Classの説明
- liboctave document
最新の20件
2020-11-14 2012-07-26- 台形公式による数値積分
