FC2ブログ

# 指数積分 En(x) Exponential Integral En
# 引数 変数 変数 ($N, $X)
# 戻り値 指数積分 En(x) ($ExponentialIntegralEn)
sub EXPONENTIALINTEGRALEN{
my ($N, $X) = @_;
my $ExponentialIntegralEn = 0;
my $SimpsonsRule = 0;
my $h = 0.1;
my $i = 1;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;
my $PrevSum = 0;
my $Epsilon = 1.0e-20;

if($X == 0){
my $Temp = 0;

if($N <= 1){
$Temp = "Inf";
}else {
$Temp = (($N - int($N)) == 0 ? (1 / ($N - 1)) : 0);
}

# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEn = $Temp;

return $ExponentialIntegralEn;
}

do {
# 刻み幅 $h は適当
my $T = 1 + ($i * $h);
my $Number = ($i % 2 == 0 ? 2 : 4);
my $Temp = &INTEGRAND($N, $X, $T);

$PrevSum = $Sum;
$Sum += $Number * $Temp;
$i++;
} while(($Sum - $PrevSum) >= $Epsilon);

$SumA = &INTEGRAND($N, $X, 1);
$SumB = &INTEGRAND($N, $X, (($h * $i) - $X));

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 指数積分 En(x) Exponential Integral
$ExponentialIntegralEn = $SimpsonsRule;

return $ExponentialIntegralEn;
}

# 被積分関数 Integrand
# 引数 変数 変数 変数 ($N, $X, $T)
# 戻り値 被積分関数 ($Integrand)
sub INTEGRAND{
my ($N, $X, $T) = @_;
my $Integrand = 0;

# 被積分関数 Integrand
$Integrand = exp(-($X * $T)) / ($T ** $N);

return $Integrand;
}


参考URL
指数積分 - Wikipedia
指数積分En(x) - 高精度計算サイト

一言
Xが0の場合は、このやり方でいいのかな?
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

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

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

検索フォーム


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