ガウス ルジャンドル積分公式の係数の計算
long int、関数オブジェクト、テンプレート、例外など
通常は、double 型を用いていても、公式の係数を求めるときはdouble 型より１つ精度の高いlong int
型で計算を行って求め、それをdouble 型の計算で用いると精度がより安心である。定数としてdouble 型の精度を超える桁数の数値を書いても、コンパイラーによって四捨五入されてdouble
型に変換される。係数を求める例として、ガウス・ルジャンドルの積分公式における係数を求めるプログラムを作成してみた。分点が n 個である積分公式の分点
は、次式の根として与えられる（岡本安晴「Delphiで学ぶデータ分析法」ＣＱ出版社、1998、pp. 236-238）。
分点の重みは次式によって与えられる。
上式により
の場合のプログラムを作成した。関数オブジェクト、テンプレート、例外の生成（下記のプログラムでは、例外は生成されないが）も使用してみた。プログラムは圧縮ファイルSource.zipとしてまとめたので、圧縮ファイル名Source.zipをクリックダウンロード解凍すればエディタで開くことができる。プログラムをFedora23で実行すると以下のようになる。
Source.cppの内容は以下のとおりである。
#include <iostream>
#include <iomanip>
#include <sstream>
#include <cmath>
using namespace std;
inline long double sqr(
long double v ){
return
v * v;
}
class CP {
public:
int n;
CP(){}
CP(int vn): n(vn){}
~CP(){}
void
set_n(int vn) {
n
= vn;
}
long
double pfunc(long double x, int
k) {
if
(k > 1) {
return
((2.0 * k - 1.0) * x * pfunc(x, k - 1)
-
(k - 1.0) * pfunc(x, k - 2)) / k;
}
else
if (k == 1) {
return
x;
}
else
if (k == 0) {
return
1.0;
}else{
ostringstream sout;
sout << "Error at pfunc...k
= " << k;
throw
sout.str();
}
}
long
double operator()(long double x) {
return
pfunc(x, n);
}
};
class CDP {
public:
int n;
CDP()
{}
CDP(int vn) : n(vn)
{}
~CDP()
{}
void
set_n(int vn) {
n
= vn;
}
long
double pfunc(long double x, int
k) {
if
(k > 1) {
return
((2.0 * k - 1.0) * x * pfunc(x, k - 1)
-
(k - 1.0) * pfunc(x, k - 2)) / k;
}
else
if (k == 1) {
return
x;
}
else
if (k == 0) {
return
1.0;
}
else
{
ostringstream sout;
sout << "Error at pfunc...k
= " << k;
throw
sout.str();
}
}
long
double dpfunc(long double x, int
k) {
if
(k > 1) {
return
((2.0 * k - 1.0) * (pfunc(x, k - 1) + x * dpfunc(x, k
- 1))
-
(k - 1) * dpfunc(x, k - 2)) / k;
}
else
if (k == 1) {
return
1.0;
}
else
if (k == 0) {
return
0.0;
}
else
{
ostringstream sout;
sout << "Error at dpfunc...k
= " << k;
throw
sout.str();
}
}
long
double operator()(long double x) {
return
dpfunc(x, n);
}
};
// Bisection法による根の計算
template<class TF>
void TBisection(TF
& f,
// 関数
long
double s,
// fが増加関数のとき s > 0.0
long
double c,
// 根に対する関数値
long
double L_b, long double U_b, // 根の存在区間の下限と上限
long
double & Root,
// 根 f(Root) = c
long
double acc)
// 精度
{
if
((f(L_b) - c) * (f(U_b) -
c) > 0.0) {
ostringstream sout;
sout << "Invalid combination of L_b and U_b at TBisection.";
throw
sout.str();
}
do
{
long
double m = 0.5 * (L_b + U_b);
long
double v = f(m);
if
(s * (v - c) > 0.0) {
U_b = m;
}
else
{
L_b = m;
}
}
while (abs(U_b - L_b) >=
acc);
Root
= 0.5 * (L_b + U_b);
}
// Newton法による根の計算
template<class TF, class TDF>
void TNewton( long double & x, // 根の初期値、戻り値は根 f(x)=c
long double c, // 関数値
TF & f, // 関数
TDF & df, // 導関数
long double acc) // 精度
{
int NStep = 0;
int ck = 0;
long
double x0;
do
{
x0
= x;
x
= x0 + (c - f(x0)) / df(x0);
}
while (abs(x - x0) >= acc);
return;
}
int main() {
try
{
CP
p(9);
CDP
dp(9);
long
double xpos0 = -1.0 + 0.01;
cout << setprecision(19)
<< scientific;
for
(int istep = 0; istep < 9; istep++) {
cout << "istep =
" << istep << endl;
long
double ck_v = p(xpos0);
long
double xpos1 = xpos0;
long
double v;
do
{
xpos1
+= 0.01;
v
= p(xpos1);
if
(v * ck_v < 0.0) break;
}
while (true);
long
double xi;
TBisection<CP>(p, v - ck_v, 0.0, xpos0, xpos1, xi,
1e-7);
TNewton<CP, CDP>( xi, 0.0, p, dp,
1.0e-19);
cout << "xi = " << setw(30) << left << xi;
long
double wi = 2.0 / ((1.0 - sqr(xi))
* sqr(dp(xi)));
cout << "wi = "
<< wi << endl;
xpos0
= xi + 0.01;
}
}
catch
(string s){
cout << s << endl;
}
return
0;
}