シンプソン公式のプログラムを書きました。
シンプソン公式の特徴は下記の様な感じです。
- 台形公式より精度が高い
- 小区間[ , ]で、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つの方法で、を区間[0, 1]で積分させています。
この積分はになります。
結果は以下のようになりました。
PI/4 = 0.785398163397448
trapezoidal: 0.792893996730783
simpson : 0.785398163397438