FC2ブログ
最初
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/RombergsMethod2.txt

修正1 5/5
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/RombergsMethod3.txt

修正2 5/17

# ロンバーグ積分 Romberg's Method
# 引数 左端a 右端b 分割数 被積分関数 ($a, $b, $N, \@Integrand)
# 戻り値 ロンバーグ積分 ($RombergsMethod)
sub ROMBERGSMETHOD{
my ($a, $b, $N, $Integrand) = @_;
($a, $b) = ($b, $a) if($a > $b);
my $RombergsMethod = 0;
my @Romberg = ();
my $Sum = 0;
my $h = 0;
my $Power4 = 0;
my $Epsilon = 1.0e-10;
my $Count = @$Integrand - 1;

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

# (0, 0)
$Romberg[0] = (($b - $a) / 2) * (&INTEGRAND($a, $Integrand) + &INTEGRAND($b, $Integrand));

# 計算
for(my $i = 1; $i < $N; $i++){
$Power4 = 1;

for(my $j = $i; $j >= 0; $j--){
if($i == $j){
$Sum = 0;
$h = ($b - $a) / (2 ** $j);

for(my $k = 1; $k <= (2 ** ($i - 1)); $k++){
$Sum += &INTEGRAND($a + (((2 * $k) - 1) * $h), $Integrand);
}

# ($i, 0)
$Romberg[$j] = ($Romberg[$j - 1] / 2) + ($Sum * $h);
}else {
# 累乗
$Power4 = 4 * $Power4;

# リチャードソンの補外 Richardson's Extrapolation
$Romberg[$j] = (($Power4 * $Romberg[$j + 1]) - ($Romberg[$j])) / ($Power4 - 1);

# 収束判定
if(abs($Romberg[$j] - $Romberg[$j + 1]) < $Epsilon){
$Romberg[0] = $Romberg[$j];
$i = $N;
last;
}
}
}
}

# ロンバーグ積分 Romberg's Method
$RombergsMethod = $Romberg[0];

return $RombergsMethod;
}

# 被積分関数 Integrand
# 引数 変数 被積分関数 ($X, \@Integrand)
# 戻り値 被積分関数 ($Function)
sub INTEGRAND{
my ($X, $Integrand) = @_;
my $Function = 0;
my $Degree = @$Integrand - 1;

for(my $i = 0; $i <= $Degree; $i++){
$Function += $$Integrand[$i] * ($X ** ($Degree - $i));
}

return $Function;
}


例 x^2 定積分の近似値

use warnings;
use strict;

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

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


eval関数を使用して式から値を出せるようにした。
http://blog-imgs-36.fc2.com/a/m/a/amamiyaprog/RombergsMethod1.txt

参考URL
Romberg's method - Wikipedia, the free encyclopedia

修正1 5/5
収束判定を付けた

修正2 5/17
Power4の出し方が縦列ではなく、横列だったので
$Power4 = (4 ** $i); -> $Power4 = (4 ** ($i - $j)); に修正
$Power4 = (4 ** $i); -> $Power4 = 4 * $Power4; に修正
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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