W tym rozdziale zajmiemy się zadaniem obliczenia przybliżonego całek postaci:
dla funkcji ,
czy ogólniej:
dla danej wagi.
Funkcja octave'a quad()
Do obliczania przybliżonego całek jednowymiarowych tzn. typu
służy funkcja octave'a:
Jej najprostsze wywołanie to:
gdzie to początek i koniec odcinka, a
to wskaźnik (uchwyt) do funkcji
octave'a postaci:
#komendy octave'a
y=.....;
endfunction
Jeśli funkcja, której całkę chcemy obliczyć jest już wcześniej zdefiniowana, to uchwyt do niej zwraca operator octave'a: @
, np. uchwyt do funkcji zwróci komenda:
@sin
.
Na początek kilka prostych przykładów całek, które sami potrafimy obliczyć analitycznie:
- całka z
na
Otrzymaliśmy zero - jak tego oczekiwaliśmy.
-
Otrzymaliśmy wynik zgodny z oczekiwaniami.
- całka z funkcji mało regularnej
:
Sprawdźmy
c-2/3
, otrzymaliśmy, czyli zero na poziomie błędu zaokrągleń.
Czy funkcja quad()
zawsze daje dobre wyniki?
Rozpatrzmy prosty przykład
całki z funkcji z parametrem
na , por. rysunek [link], dla
Całka powinna być zawsze równa jeden.
Zdefiniujemy najpierw funkcję :
Otrzymaliśmy jeden. Powtórzmy obliczenia dla dla
:
Od pewnego momentu obliczone przybliżenia całek są równe zero zamiast jeden.
Wychwyćmy ten moment dokładniej.
Powtórzmy obliczenia dla z
.
Jest to dość gwałtowna zmiana, por. rysunek [link].
Dlaczego tak jest? Oczywiście funkcja octave'a quad(()
może obliczyć przybliżenie całki korzystając z co najwyżej skończonej ilości wartości funkcji na odcinku całkowania. Nośnik funkcji jest bardzo
mały dla względem odcinka całkowania
, stąd być może wszystkie obliczone wartości były równe zero.
Można się domyśleć, że funkcje octave'a działają na zasadzie czarnej skrzynki: podajemy parametry do funkcji, a dana funkcja, w tym przypadku quad()
powinna zwrócić przybliżenie całki.
Spróbujmy obliczyć tę całkę
dzieląc odcinek na małe pododcinki i całkując po nich tzn.
a każdą całkę możemy obliczyć przy pomocy funkcji
octave'a quad()
.
Otrzymaliśmy błąd równy
ans = 6.8505e-12
,
podczas gdy
abs(quad(g ,-1,1)-1)
zwraca błąd równy jeden.
Oczywiście znając nośnik funkcji możemy policzyć całkę po jej nośniku.
Otrzymamy:
Błąd jest równy ans = 1.1102e-16
, czyli jest na poziomie błędu zaokrągleń.
Również dla funkcji o dużej zmienności, np. silnie oscylujących, funkcja
octave'a quad(()
może zwrócić zły wynik.
Spróbujmy scałkować np. na odcinku
i porównajmy z dokładnym wynikiem
.
Tu funkcja quad()
wydrukowała na ekranie ostrzeżenie:
Warto wspomnieć, że funkcja quad()
oblicza również całki po odcinkach nieograniczonych. Za wartości w wywołaniu
quad(f,a,b)
możemy przyjąć
inf
lub odpowiednio -inf
, np. możemy scałkować funkcję po całej prostej rzeczywistej, czy jakiejś półprostej:
Czy wynik jest zgodny z oczekiwaniami?
Kwadratury złożone
Oczywiście w octave'ie możemy zaimplementować
najprostsze kwadratury złożone, np. najprostszą kwadraturę złożoną:
Implementacja w octave'ie tej kwadratury to:
#kwadratura prostokatow prawostronna
#Input: a,b, konce przedzialu calkowania
#N- ilosc punktow - domyslnie 100
#f- wskaznik do funkcji -
# zwracajacej wektor wartosci dla wektora danych
#Output: c - przyblizenie calki
h=(b-a)/N;
x=(1:N)*h;
y=f(x);
c=h*sum(y);
endfunction
Przetestujmy tę kwadraturę. Na początku na jednomianach:
prostakw(@(x) x.*x,0,1)-1/3
Wynik jest poprawny. Jeżeli zwiększymy , to być może uzyskamy lepszy wynik:
prostakw(@(x) x.*x,0,1,1000)-1/3
Błąd jest na poziomie .
Dla funkcji klasy
błąd można oszacować przez
przy czym dla funkcji w tym oszacowaniu widzimy równość.
Można pokazać, że dla dostatecznie gładkiej funkcji zachodzi:
dla pewnego punktu , a dla bardziej regularnych funkcji
dla pewnej ujemnej stałej .
Badanie rzędu kwadratury dla całki z funkcji
na
. Błąd
. (#)
Przetestujmy teraz tę kwadraturę dla ustalonej analitycznej funkcji i podwajanych
wartości , tzn. obliczymy
dla
i następnie
policzymy:
dla .
Badanie rzędu kwadratury dla całki z funkcji
na
. Błąd
. (#)
#c=2; #znana wartosc calki
e=0;
for k=1:M,
kw=prostakw(f,a,b,N);
ep=e;
e=abs(c-kw);
printf("[%d] e=%g ep/e=%6.5f \n",N,e,ep/e);
N*=2;
endfor
Wyniki dla funkcji i odcinka
są zawarte w Tabeli [link], a dla funkcji
i
całki z odcinka w Tabeli [link].
Dlaczego dla wyniki wskazują na błąd
?
Badanie rzędu kwadratury dla całki z funkcji
na
. Błąd
. (#)
Zmieńmy odcinek na i wyniki będą zgodne z oczekiwanymi, por. Tabela [link].
Normy całkowe
Zastanówmy się, jak w sposób przybliżony wyznaczyć normę typu na odcinku
dla zadanej funkcji
zdefiniowanej w octave'ie jako
Spróbujmy zdefiniować funkcję octave'a, która korzystając z quad()
wyznaczy przybliżenie normy .
Jako parametry tej funkcji chcemy podawać wskaźnik do funkcji typu y=f(x)
, , oraz końce odcinka
.
Żeby wywołać funkcję quad()
musimy podać funkcję octave'a jednoargumentową jako pierwszy argument.
Trzeba ją odpowiednio zdefiniować wewnątrz naszej funkcji. Najwygodniejszym sposobem jest użycie mechanizmu funkcji anonimowej, por Rozdział [link]:
oczywiście musi być wcześniej zdefiniowane.
Tak więc można by napisać tę funkcję następująco:
Sprawdźmy dla i prostych funkcji
np.:
Wyniki są poprawne.
Mając zaimplementowaną taką funkcję, możemy np. przeprowadzić badanie rzędów zbieżności interpolacji splajnowej w normie , zamiast w normie maksimum, tak jak w Rozdziale [link].