C++ 用の拡張として、数学関係のライブラリをちょこっと作ってみることにしました。
今回はその第1弾で、積分の計算です。
公式に当てはめるだけなので、誰でもできることのなのですが、一応作りました。
シンプソンの公式および台形公式の2種類を用意しました。
いまのところ、ベクトルの積分と、多変数の積分は対応していません。
とりあえず計算結果を利用するだけなら、Mathematica や GNU Octave のような近似計算できる計算ソフトを利用した方が便利だと思います。
僕はそんなに数学に通なつもりはないので(いちおう工学系だし)、計算途中で数学的におかしな計算をしているかもしれません。発見したら指摘してください。
f が積分する関数。別途、定義して、その関数へのポインタを指定してください。
a, b が積分範囲です。とりあえず f > 0 で a < b ならば積分結果は正になります。
もちろん、b < a にすると、符号は逆転する・・・はずです(ためすの忘れてた・・・)
M はシンプソンの公式ならfの4階微分、台形ならfの2階微分の絶対値最大値以上の数値を入れます。
普通、そんなの計算面倒でしょうから、適当にいれても計算できます。
それでも、結果には影響するので、本当は求めた方がいいかも。
n は分割数です。nを大きくすると計算に時間がかかるようになります。
integral.h
#ifndef _INTEGRAL_H_ #define _INTEGRAL_H_ class integ_err { public: int code; integ_err(int i) { code = i; } operator int() { return code; } }; //----------------------------------------------------------------------------- // This pow is not for real number power. Use pow in the cmath. // // This is used for class T (not int, double, ...) // And this function does not check overflow. //----------------------------------------------------------------------------- template <class T> T pow(T a, int b) { if(b <= 1) return T(0); for(; b > 1; b--) { a = a * a; // a *= a; (May does not have *= operator.) } return a; } //----------------------------------------------------------------------------- // Simpson's formula using Integral Calculation // // f must be able to differentiate 4 times. // M must be over the absolute value of d^4f/dx^4(x) // // type T must be able to transform from integer and operator+, // operator-, operator*, operator/ must be defined. (if you use class) //----------------------------------------------------------------------------- template <class T> T Simpson(T (*f)(T), T a, T b, T M, int n) { T h = (b - a) / T(2 * n); T ret = T(0); for(int i = 0; i <= 2 * n; i++) { if(i == 0) ret = ret + f(a); else if(i % 2 == 1) ret = ret + (T(4) * f(a + h * T(i))); else if(i == 2 * n) ret = ret + f(b); else if(i % 2 == 0) ret = ret + (T(2) * f(a + h * T(i))); else throw (integ_err) -1; } ret = (ret * (h / T(3))) + (::pow(b - a, 5) * M / (T(2880) * ::pow(T(n), 4))); // This pow is not pow of cmath. return ret; } //----------------------------------------------------------------------------- // Trapezium formula using Integral Calculation // // f must be able to differentiate twice. // M must be over the absolute value of d^2f/dx^2(x) // // type T must be able to transform from integer and operator+, // operator-, operator*, operator/ must be defined. (if you use class) //----------------------------------------------------------------------------- template <class T> T Trapezium(T (*f)(T), T a, T b, T M, int n) { T h = (b - a) / T(n); T ret = T(0); for(int i = 0; i <= n; i++) { if(i == 0) ret = ret + f(a); else if(i == n) ret = ret + f(b); else ret = ret + (2 * f(a + h * T(i))); } ret = (ret * (h / T(2))) + (::pow(b - a, 3) * M / (T(12) * ::pow(T(n), 2))); return ret; }
(2008/06/05 追記) 変数 ret の初期化を忘れてました(滝汗
GCC (GNU C++) を使う場合に注意点
GCC の long double 型は 96 ビット(VC++ や GCC の double は 64 ビット)なので、一部の数学関数が対応していないか、なんなのかよくわかりませんが、とりあえず、計算できない場合が存在します。
計算結果が得られなかった場合、
cout << (long double) Simpson(...) << endl; // なぜかキャストが必要
は nan *1と表示されます。
*1:Not A Number の略
