Friday, June 10, 2016

całkowanie numeryczne

W ramach zaliczenia pewnego przedmiotu trzeba było sobie wybrać zadanie i wykonać projekt. Wybrałem całkowanie numeryczne, bo uznałem to za ciekawy temat. Napisałem względnie prosty skrypt wyznaczający przybliżoną wartość całki oznaczonej funkcji jednej zmiennej czterema metodami. Projekt miał za zadanie wypisać obliczenia, błąd bezwzględny oraz czas wykonywania każdej z metod i wyniki miały zostać przedyskutowane. Do jednej metody (Romberga) wykorzystałem zewnętrzną bibliotekę oraz przyjąłem wartości zwracane przez tę funkcję biblioteczną za wzorcowe, ponieważ, według moich obserwacji, zdaje się ona podawać wynik najbliższy analitycznemu (kiedy to tylko możliwe).

Skrypt wykorzystuje 4 metody:
  • prostokątów
  • trapezów
  • Simpsona (parabol)
  • Romberga

Pomijając metodę Romberga, w przypadku 3 pozostałych metod przedział, w którym całka ma być wyznaczona, jest dzielony na N podprzedziałów. Skrypt zakłada przy tym, że funkcja, którą mu podamy istnieje i że jest ciągła w każdym punkcie przedziału < a ; b >, gdzie a oznacza dolną granicę całkowania a b - górną.

Metoda prostokątów i coraz dokładniejsze dopasowanie wraz ze wzrostem podziałów:


Credit: Khurram Wadee (https://en.wikipedia.org/wiki/File:Rectangle_rule.gif)

Metoda trapezów:



Skrypt przyjmuje 4 parametry: dolną (a) oraz górną (b) granicę całkowania, ilość podprzedziałów (N) oraz wzór funkcji jednej zmiennej.
Napisany w Perlu. Nie wrzucam tu wszystkich wzorów - jeżeli interesuje Cię ten temat, z łatwością odszukasz je w Internecie. Chciałbym bardziej skupić się na przykładowym rezultacie działania.


(Przypominam, że przyjmuję w swoich rozważaniach wartość metody Romberga jako wzorzec. Mowa o całkowaniu numerycznym, więc operujemy na przybliżeniach.)

Nasunęło mi się kilka wniosków w związku z powyższymi danymi:
  • pomijając metodę biblioteczną Romberga, metoda Simpsona daje najbardziej dokładny wynik
  • czas liczenia (za wyłączeniem metody Romberga) jest mniej więcej liniowy
  • z każdym 10-krotnym zwiększeniem ilości podprzedziałów wyniki stają się o 2 rzędy wielkości dokładniejsze
  • przy znacznie zwiększonej liczbie podziałów metoda prostokątów jest nieznacznie szybsza od pozostałych

To tylko mały "skrawek" zagadnienia. Funkcja, którą wykorzystałem jako przykład jest względnie gładka (zob. wykres). Problemem byłoby policzenie całki oznaczonej "stromej" funkcji, która "bardzo skacze" i do czego wymagana byłaby duża ilość podziałów, co wiązałoby się z wydłużeniem czasu trwania obliczeń. Wykorzystywane podejścia dzielą przedział na równe podprzedziały, bez badania przebiegu zmienności funkcji i "inteligentnego" ich ustalania lub dobierania metody w zależności od tego, jaki funkcja "ma mniej więcej kształt".

Do kodu nie roszczę sobie żadnych praw, możesz z niego korzystać w dowolny sposób, przy czym nie ponoszę odpowiedzialności za to, co nim zrobisz. Na pewno może być napisany lepiej, ale pierwotne założenie (przeważnie prostoty) spełnia.
#!/usr/bin/perl
use Math::Integral::Romberg qw(integral);
use POSIX;
use Time::HiRes qw(time);
use Scalar::Util qw(looks_like_number);

use constant PI => 4 * atan2(1,1);
use constant e => 2.71828182845904523536;

sub wrapper {
 my ($x) = @_;
 eval $funkcja;
}

sub rectangles {
 $lower = $_[0];
 $upper = $_[1];
 $intervals = $_[2];
 $step = ( $upper - $lower ) / $intervals;

 for($x = $lower; $x < $upper; $x += $step){
  $a = $step;
  $b = eval $funkcja;
  $sum += $a*$b;
 }
 return $sum;
}

sub trapezoids {
 $lower = $_[0];
 $upper = $_[1];
 $intervals = $_[2];
 $step = ( $upper - $lower ) / $intervals;

 $temp_sum = 0;
 $overall_sum = 0;
 $iterator = 0; #zliczanie "przebytych" podzialow; ich dlugosc nie musi (i najczesciej nie bedzie) liczba calkowita
 for($x = $lower; $x < $upper; $x += $step, $iterator++){
  $temp_sum = eval $funkcja;
  $overall_sum += ( $iterator==0 ? ($temp_sum/2) : ( $iterator==($intervals-1) ? ($temp_sum/2) : ($temp_sum) ) );
 }
 $overall_sum *= $step;
 return $overall_sum;
}

sub simpsons {
 $lower = $_[0];
 $upper = $_[1];
 $intervals = $_[2];
 $step = ( $upper - $lower ) / $intervals;

 $iterator = 0;
 $k = 0;
 $overall_sum = 0;
 for($x = $lower; $x <= $upper; $x += $step, $iterator++){
  $temp_sum = eval $funkcja;
  if($iterator==0 || $iterator==$intervals){
   $overall_sum += $temp_sum;
  } else {
   $k ^= 1;
   $overall_sum += ($k+1) * 2 * $temp_sum;
  }
 }
 $overall_sum *= $step/3;
 return $overall_sum;
}

sub parse_constants {
 $tempstr = $_[0];
 $tempstr =~ s/pi/PI/g;
 $tempstr =~ s/Pi/PI/g;
 $tempstr =~ s/pI/PI/g;
 #$tempstr =~ s/e/e/g;
 #print "\n $tempstr \n";
 $tempstr = eval $tempstr;
 return $tempstr;
}

sub call_redirector {
 $start = time;
 $sum = $_[0]($_[2], $_[3], $_[4]);
 $end = time;
 ${$_[1]} = ( $end - $start ) * 1000;
 return $sum;
}

sub wypisz_wynik_metody {
 print " metoda ", ($_[0]), "\t: ", (${$_[1]});
 if ($cmpval != 0) {
  $bl_bezwzgl = abs( ${$_[1]} - $cmpval );  
  print "\tblad bezwzgledny: $bl_bezwzgl";
  if ($bl_bezwzgl == 0) { print "\t\t"; }
 }
 print "\t\tczas: ", (${$_[2]}), " ms\n\n ";
 return;
}

$n_args = @ARGV;
if ($n_args < 4) { die "Za malo parametrow.\n"; }

$lower = parse_constants(shift); #$ARGV[0]
$upper = parse_constants(shift); #$ARGV[1]
$intervals = shift; #$ARGV[2]

$not_numeric1 = "Podana wartosc ";
$not_numeric2 = " zdaje sie nie byc wartoscia liczbowa.\n";

if(!looks_like_number($lower)) { die "$not_numeric1" . "dolnej granicy" . "$not_numeric2"; }
if(!looks_like_number($upper)) { die "$not_numeric1" . "gornej granicy" . "$not_numeric2"; }
if(!looks_like_number($intervals)) { die "$not_numeric1" . "ilosci interwalow" . "$not_numeric2"; }

unless ($intervals > 0) { die "Ilosc interwalow musi byc dodatnia.\n"; }

for($i = 0; $i < $n_args - 3; $i++ ) { #n_args minus 3, bo 3 parametry zostały już zdjęte z ARGV
 if ($i==0) {
  $funkcja = shift;
 } else {
  $funkcja = $funkcja . ' ' . shift;
 }
}
$_funkcja = $funkcja;
our $funkcja =~ s/x/\$x/g;

print "\n dolna granica calkowania: $lower\n";
print " gorna granica calkowania: $upper\n";
print " ilosc interwalow : $intervals\n";
print " funkcja: $_funkcja\n";

#tak sie sklada, ze wszystkie wlasciwe funkcje liczace przyjmuja po 3 parametry
$czas_prostokaty = 0;
$suma_prostokaty = call_redirector(\&rectangles, \$czas_prostokaty, $lower, $upper, $intervals);

$czas_trapezy = 0;
$suma_trapezy = call_redirector(\&trapezoids, \$czas_trapezy, $lower, $upper, $intervals);

$czas_simpson = 0;
$suma_simpson = call_redirector(\&simpsons, \$czas_simpson, $lower, $upper, $intervals);

$czas_romberg = 0;
$suma_romberg = call_redirector(\&integral, \$czas_romberg, \&wrapper, $lower, $upper);

our $cmpval = $suma_romberg; #funkcja biblioteczna zwraca bardzo dokladne wartosci; jej wyjscie posluzy jako wartosc wzorcowa

print "\n przyblizona wartosc calki oznaczonej:\n\n ";
wypisz_wynik_metody("Romberga", \$suma_romberg, \$czas_romberg);
wypisz_wynik_metody("prostokatow", \$suma_prostokaty, \$czas_prostokaty);
wypisz_wynik_metody("trapezow", \$suma_trapezy, \$czas_trapezy);
wypisz_wynik_metody("Simpsona", \$suma_simpson, \$czas_simpson);

No comments:

Post a Comment