Perlでかんたんな数値計算の入門

このエントリーをはてなブックマークに追加
ShanonAdventCalendar2018、24日目の記事をお届けします。

こんにちは。シャノンの新人エンジニアのgutsです。

今年は新人として、学ぶことの多い年でした。
何かを新しく学ぶときには、これまで自分が慣れ親しんだものと比較しながら学ぶとやりやすいかと思います。
そこで本日は、私が大学でよくやっていた数値計算のうちの一つ、積分の計算の簡単なものをPerlでやってみたいと思います。
数値計算に慣れ親しんでいるけれど、Perlは初めて、という方のニッチな(?)需要を満たしたいと思います。
具体例としては、下記の定積分を行います。
左辺の積分を実行したものが右辺になるということは、ここでは本題とずれるので、省略します。理工系の大学生・大学受験生なら教科書で見たことがあるのではないかと思います。

左辺の積分をプログラムで実行させ、その答えが計算誤差(詳細は省きますが必ず発生します)の範囲内で右辺になってくれれば、正しくプログラムできたことになります。

プログラムで積分をどうやってやるかについては、いろいろと手法がありますが、ここでは素朴で基本的な手法である「台形則」を用います。
その考え方は
1. 被積分関数を細かい部分に分割して、各部分での面積を合計する
2. その際に各部分を台形として近似する
です。
1. で関数を分割する際、分割する数が多いほど誤差が小さく、正しい結果に近づきます。
さて、積分実行プログラムをPerlで書いたものがこちらです。
最初は関数の分割数を1から開始し、計算誤差が小数点以下10桁以下になるまで、分割数を倍に増やしています。

#!/usr/bin/env perl
use Math::Trig 'pi';

my $pi_over_4 = pi / 4.0e+0;
my $result_diff = abs($pi_over_4);
my $result_integral = 0.0e+0;
my %params = (
    'left'  => 0.0e+0,# 積分区間の左端
    'right' => 1.0e+0,# 積分区間の右端
    'div'   => 1,     # 区間内分割数(初期値:1)
);
while ($result_diff > 1.0e-10) { # 誤差が10^(-1)以下になるまで繰り返す
    $params{div} *= 2;
    $result_integral = trapezoidal(\%params, \&function);
    $result_diff = abs($pi_over_4 - $result_integral);
}
print "右辺 = \t$pi_over_4\n";
print "左辺 = \t$result_integral\n";
print "誤差:   $result_diff\n";
print "分割数: $params{div}\n";

# 台形則による積分
# 引数
#   $params : パラメータを入れたhashref
#             left  : 積分区間左端
#             right : 積分区間右端
#             div   : 積分区間分割数
# 戻り値
# 積分結果
sub trapezoidal {
    my ($params, $function) = @_;
    my $sum = 0.0e+0;
    for (my $i = 0; $i < $params->{div}; $i++) {
        $sum += 0.5e+0 * ($function->(x_position($i, $params))
                            + $function->(x_position($i + 1, $params)));
    }
    $sum /= $params->{div};
    return $sum;
}

# 被積分関数
# 引数
#   $x : x
# 戻り値
#   引数xに対する関数の値
sub function {
    my ($x) = @_;
    return 1.0e+0 / (1.0e+0 + $x ** 2);
}

# x座標 (x_i)
# 引数
#   $i : 左端から数えてi番目を意味する 
#   $params : パラメータを入れたhashref
#             left  : 積分区間左端
#             right : 積分区間右端
#             div   : 積分区間分割数
# 戻り値
#   左端から数えてi番目のx座標
sub x_position {
    my ($i, $params) = @_;
    return $params->{left} + ($params->{right} - $params->{left}) * $i / $params->{div};
}

実行結果はこちら↓

右辺 =  0.785398163397448
左辺 =  0.785398163358635
誤差:   3.88130638739881e-11
分割数: 32768

右辺と左辺がほとんど一致していることがわかります。
小数点以下11桁目からは食い違っていますが、これが計算誤差です。
計算誤差は自分の目論見通りになっているので、この積分は正しく実行できていることがわかりました。 この誤差にするために必要だった分割数は、だいたい33000程度ということですね。

サブルーチンfunctionの中身や、積分区間をいろいろと変えれば、自分のやってみたい積分を行えます。
ただし、なんでもかんでも積分できるわけではなく、関数の形によってはうまくいかないことも結構ありますが、その辺の話は入門ではなさそうなので、ここでは割愛します。

最後に、数値計算は主にFortran, C言語, C++のどれかを用いて行われていることが多いです。また、最近はPythonを用いた計算も流行っているようです。私はPythonはまだ触ったことがありませんが、もっと簡単に計算を行うことができるようです。
Perlを使って数値計算というのはあまり見かけないとは思いますが、数値計算はある程度慣れているけれどPerlは初めて、という方の参考になれればと思います。

それでは引き続き、Shanon Advent Calendar 2018を引き続きお楽しみください。
次の記事
« Prev Post
前の記事
Next Post »
Related Posts Plugin for WordPress, Blogger...