読者です 読者をやめる 読者になる 読者になる

シンプソン公式

perl 数値計算

シンプソン公式のプログラムを書きました。
シンプソン公式の特徴は下記の様な感じです。

  • 台形公式より精度が高い
  • 小区間[ x_n , x_{n+2} ]で、2次曲線を近似して積分する


シンプソン公式のプログラムは下記のようにしました。

sub simpson {
    my ($f, $N, $a, $b) = @_;
    my (@x, @y);

    if ($N % 2 == 1) {
        $N += 1;
    }

    my $S = 0.0;
    my $h = ($b-$a)/$N;

    for my $i (0..$N) {
        $x[$i] = $a + $i*$h;
    }

    for my $i (0..$N) {
        $x[$i] = $a + $i*$h;
    }

    for (my $i=0; $i<$N; $i=$i+2) {
        $y[$i]   = $f->($x[$i]);
        $y[$i+1] = $f->($x[$i+1]);
        $y[$i+2] = $f->($x[$i+2]);
        $S += $h*($y[$i]+4.0*$y[$i+1]+$y[$i+2])/3.0;
    }

    return $S;
}


積分する関数、分割数、積分区間を渡す関数となっています。


台形形式と比べて精度が高いという事なので、比べてみました。

use strict;
use warnings;
use Math::Trig qw( pi );

my $ans = pi() / 4.0;

my $f = sub {
    my $x = shift || 0;
    return 1.0/($x*$x+1);
};

my $N = 100;
my $a = 0.0;
my $b = 1.0;

print "PI/4 = $ans\n";
print "trapezoidal: ",trapezoidal($f, $N, $a, $b), "\n";
print "simpson    : ",simpson($f, $N, $a, $b), "\n";


sub simpson {
    my ($f, $N, $a, $b) = @_;
    my (@x, @y);

    if ($N % 2 == 1) {
        $N += 1;
    }

    my $S = 0.0;
    my $h = ($b-$a)/$N;

    for my $i (0..$N) {
        $x[$i] = $a + $i*$h;
    }

    for (my $i=0; $i<$N; $i=$i+2) {
        $y[$i]   = $f->($x[$i]);
        $y[$i+1] = $f->($x[$i+1]);
        $y[$i+2] = $f->($x[$i+2]);
        $S += $h*($y[$i]+4.0*$y[$i+1]+$y[$i+2])/3.0;
    }

    return $S;
}



sub trapezoidal {
    my ($f, $N, $a, $b) = @_;

    my (@x, @y);

    my $h = ($b-$a)/$N;
    for my $i (0..$N) {
        $x[$i] = $a + $i*$h;
    }

    my $S = 0.0;
    for my $i (0..$N) {
        $y[$i] = $f->($x[$i]);
        $y[$i+1] = $f->($x[$i+1]);
        $S = $S + $h*($y[$i] + $y[$i+1])/2.0;
    }

    return $S;
}


2つの方法で、f(x)=\frac{1}{1+x^2}を区間[0, 1]で積分させています。
この積分\frac{\pi}{4}になります。


結果は以下のようになりました。

PI/4 = 0.785398163397448
trapezoidal: 0.792893996730783
simpson : 0.785398163397438