FC2ブログ

# ガウス-ルジャンドル数値積分 Legendre-Gauss Quadrature
# 引数 左端a 右端b 分割数 被積分関数 ($a, $b, $N, \@Integrand)
# 戻り値 ガウス-ルジャンドル数値積分 ($LegendreGaussQuadrature)
sub LEGENDREGAUSSQUADRATURE{
my ($a, $b, $N, $Integrand) = @_;
($a, $b) = ($b, $a) if($a > $b);
my $LegendreGaussQuadrature = 0;
my @NewtonsMethod = ();
my $LegendrePolynomial = 0;
my $X = 0;
my $Sum = 0;
my $AB2 = ($a + $b) / 2;
my $BA2 = ($b - $a) / 2;
my $Degree = @$Integrand - 1;
my $Count = @$Integrand - 1;

# 配列数と分割数の確認
if(($Count < 0) || ($N < 2)){
return 0;
}

# 計算
for(my $i = 1; $i <= $N; $i++){
my $tmp = 0;
# ニュートン法 Newton's Method
@NewtonsMethod = &NEWTONSMETHOD($N, $i);
$X = $NewtonsMethod[0];
$LegendrePolynomial = $NewtonsMethod[1];
# 重み
my $Weight = (2 * (1 - ($X * $X))) / (($N * $LegendrePolynomial) * ($N * $LegendrePolynomial));

$X = $AB2 + ($BA2 * $X);
for(my $j = 0; $j <= $Degree; $j++){
$tmp += $$Integrand[$j] * ($X ** ($Degree - $j));
}

$Sum += $Weight * $tmp;
}

# ガウス-ルジャンドル数値積分 Legendre-Gauss Quadrature
$LegendreGaussQuadrature = $BA2 * $Sum;

return $LegendreGaussQuadrature;
}

# ニュートン法 Newton's Method
# 引数 次数 位置 ($N, $K)
# 戻り値 ニュートン法 ($NewtonsMethod)
sub NEWTONSMETHOD{
my ($N, $K) = @_;
my @NewtonsMethod = ();
my @LegendrePolynomial = ();
my $Pi = atan2(1, 1) * 4;
my $X = cos(($Pi * ($K - 0.25)) / ($N + 0.5));
my $PrevX = 0;
my $Limit = 100;
my $Epsilon = 1.0e-20;

# 計算
for(my $i = 0; $i < $Limit; $i++){
# ルジャンドル多項式 Legendre Polynomial
@LegendrePolynomial = &LEGENDREPOLYNOMIAL($N, $X);
my $Num = $LegendrePolynomial[0];
my $Den = (($N * ($LegendrePolynomial[1] - ($X * $LegendrePolynomial[0]))) / (1 - ($X * $X)));

# 近似解
$PrevX = $X;
$X = $X - ($Num / $Den);
# ニュートン法 Newton's Method
# [0] 近似解 [1] 一つ前のルジャンドル多項式の値
$NewtonsMethod[0] = $X;
$NewtonsMethod[1] = $LegendrePolynomial[1];

# 収束判定
last if(abs($X - $PrevX) < $Epsilon);
}

return @NewtonsMethod;
}

# ルジャンドル多項式 Legendre Polynomial
# 引数 次数 値 ($N, $x)
# 戻り値 ルジャンドル多項式 (@LegendrePolynomial)
sub LEGENDREPOLYNOMIAL{
my ($N, $x) = @_;
my @LegendrePolynomial = ($x , 1);

# 計算
for(my $i = 1; $i < $N; $i++){
my $tmp = ((((2 * $i) + 1) * $x * $LegendrePolynomial[0]) - ($i * $LegendrePolynomial[1])) / ($i + 1);

# [0] 最新 [1] 一つ前
# ルジャンドル多項式 Legendre Polynomial
$LegendrePolynomial[1] = $LegendrePolynomial[0];
$LegendrePolynomial[0] = $tmp;
}

return @LegendrePolynomial;
}


例 x^2 定積分の近似値

use warnings;
use strict;

my $a = 0;
my $b = 1;
my $N = 2;
my @Integrand = (1,0,0);

my $LegendreGaussQuadrature = &LEGENDREGAUSSQUADRATURE($a, $b, $N, \@Integrand);
print "$LegendreGaussQuadrature\n";


eval関数を使用して式から値を出せるようにした。
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/LegendreGaussQuadrature1.txt
ルジャンドル多項式
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/LegendrePolynomial1.txt

参考URL
ガウスールジャンドルの数値積分
ガウス-ルジャンドル数値積分 - 高精度計算サイト
Legendre多項式に関する漸化式
ルジャンドル多項式(グラフ) - 高精度計算サイト
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

Author:雨宮
Firefoxを使用しているので気づかなかったけど、IE6でソースコードを上手くコピーできない

5/3
携帯用ならIE6でもソースコードをコピーできる
携帯用

検索フォーム


あわせて読みたいブログパーツ
一寸先は闇 RSS