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 $X = 0;
my $Sum = 0;
my $TempA = 0;
my $TempB = 0;
my $Pi = atan2(1, 1) * 4;
my $h = (2 * $Pi) / $N;
my $HyperbolicSine = 0;
my $HyperbolicCosine = 0;
my $Pi_2 = $Pi / 2;
my $BA_2 = (($b - $a) / 2);
my $N_2 = int(($N / 2) + 0.5);
my $Count = @$Integrand - 1;

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

# 計算
for(my $i = -$N_2; $i <= $N_2; $i++){
$HyperbolicSine = &HYPERBOLICSINE($i * $h);
$HyperbolicCosine = &HYPERBOLICCOSINE($Pi_2 * $HyperbolicSine);
$X = (($BA_2 * &HYPERBOLICTANGENT($Pi_2 * $HyperbolicSine)) + (($b + $a) / 2));
$TempA = &INTEGRAND($X, $Integrand);
$TempB = &HYPERBOLICCOSINE($i * $h) / ($HyperbolicCosine * $HyperbolicCosine);

$Sum += ($TempA * $TempB);
}

# 台形公式 Trapezoidal Rule
$TrapezoidalRule = ((($b - $a) * $Pi * $h) / 4) * $Sum;

return $TrapezoidalRule;
}

# 被積分関数 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;
}

# 双曲線正弦 Hyperbolic Sine
# 引数 値 ($X)
# 戻り値 双曲線正弦 ($HyperbolicSine)
sub HYPERBOLICSINE{
my ($X) = @_;
my $HyperbolicSine = 0;

# 双曲線正弦 Hyperbolic Sine
$HyperbolicSine = (exp($X) - exp(-$X)) / 2;

return $HyperbolicSine;
}

# 双曲線余弦 Hyperbolic Cosine
# 引数 値 ($X)
# 戻り値 双曲線余弦 ($HyperbolicCosine)
sub HYPERBOLICCOSINE{
my ($X) = @_;
my $HyperbolicCosine = 0;

# 双曲線余弦 Hyperbolic Cosine
$HyperbolicCosine = (exp($X) + exp(-$X)) / 2;

return $HyperbolicCosine;
}

# 双曲線正接 Hyperbolic Tangent
# 引数 値 ($X)
# 戻り値 双曲線正接 ($HyperbolicTangent)
sub HYPERBOLICTANGENT{
my ($X) = @_;
my $HyperbolicSine = (exp($X) - exp(-$X)) / 2;
my $HyperbolicCosine = (exp($X) + exp(-$X)) / 2;
my $HyperbolicTangent = 0;

# 双曲線正接 Hyperbolic Tangent
$HyperbolicTangent = $HyperbolicSine / $HyperbolicCosine;

return $HyperbolicTangent;
}


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

有限区間(a,b)の数値積分 - 高精度計算サイト

メモ
刻み幅 my $h = (2 * $Pi) / $N; 以外
my $h = log($N) / $N;
my $h = ln($N) / $N;
などがあった
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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