FC2ブログ

# 台形公式 Trapezoidal Rule
# 引数 左端a 右端b 分割数 被積分関数 ($a, $b, $N, \@Integrand)
# 戻り値 台形公式 ($TrapezoidalRule)
sub TRAPEZOIDALRULE{
my ($a, $b, $N, $Integrand) = @_;
($a, $b) = ($b, $a) if($a > $b);
my $TrapezoidalRule = 0;
my $h = ($b - $a) / $N;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;
my $Degree = @$Integrand - 1;
my $Count = @$Integrand - 1;

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

# 計算
for(my $i = 0; $i <= $Degree; $i++){
$SumA += $$Integrand[$i] * ($a ** ($Degree - $i));
$SumB += $$Integrand[$i] * ($b ** ($Degree - $i));
}

for(my $i = 1; $i < $N; $i++){
my $x = $a + ($h * $i);
my $tmp = 0;

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

$Sum += $tmp;
}

# 台形公式 Trapezoidal Rule
$TrapezoidalRule = ((($SumA + $SumB) / 2) + $Sum) * $h;

return $TrapezoidalRule;
}


例 x^2 定積分の近似値

use warnings;
use strict;

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

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


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

参考URL
台形公式 - Wikipedia
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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