Całkowanie funkcji metodą Simpsona(3/8)

Dyskusje dotyczące tworzenia makropoleceń, pisania skryptów oraz programowania przy użyciu UNO
Zbych
Posty: 2
Rejestracja: wt gru 19, 2017 12:21 am

Całkowanie funkcji metodą Simpsona(3/8)

Post autor: Zbych »

Dzień dobry forumowiczom

Znam 3 metody całkowania: prostokątów(kwadratów), trapezów, Simpsona(1/3). Kolejna następna jest coraz dokładniejsza i to sprawdziłem. Jest jeszcze 4-ta metoda dokładniejsza od pozostałych mianowicie metoda Simpsona(3/8). Popełniam jednak jakiś błąd którego nie mogę znaleźć.

Jeśli ktoś pomoże znaleźć mój błąd byłbym bardzo wdzięczny i załączę plik w którym porównuję powyższe cztery metody.

Pozdrawiam wszystkich Zbyszek
OpenOffice 4 Win 7
Awatar użytkownika
antekg
Posty: 18
Rejestracja: śr sie 25, 2010 6:18 pm
Lokalizacja: Warszawa

Re: Całkowanie funkcji metodą Simpsona(3/8)

Post autor: antekg »

Zbyszku, nie ma pliku! Z ciekawością bym spojrzał.
AOO 4.1.4 (Polish) na Windows 10 (64b) / AOO 4.1.0 na Windows Vista / LibreOffice na Mageia Linux
Zbych
Posty: 2
Rejestracja: wt gru 19, 2017 12:21 am

Re: Całkowanie funkcji metodą Simpsona(3/8)

Post autor: Zbych »

Witam serdecznie. Wszystkiego najlepszego, zdrowia, obfitych Świąt Bożego Narodzenia i pomyślności w roku 2018.

Jeśli odkryjesz mój błąd daj znać.

Dołączam plik porównujący metody prostokątów, trapezów, Simpsona(1/3) i Simpsona(3/8). Jak widać Simpson(3/8) jest prawie najgorszy a powinien być najlepszy.

Pozdrawiam
Zbyszek
Załączniki
-6x^2+6x.ods
(51.87 KiB) Pobrany 209 razy
OpenOffice 4 Win 7
Jan_J
Posty: 4558
Rejestracja: pt maja 22, 2009 1:20 pm
Lokalizacja: Wrocław

Re: Całkowanie funkcji metodą Simpsona(3/8)

Post autor: Jan_J »

Przy podziale przedziału [0,1] na 10 części metoda 1/3 wymaga 5 kroków długości 0,2, co zostało osiągnięte zamieszczeniem poprawnej formuły w co drugiej komórce. Trochę ręczne sterowanie, ale nie ma błędu.
Zauważ, że formułę Simpsona da się przedstawić jako średnią ważoną metody trapezów i (centralnej) metody prostokątów, obu obliczonych z dwukrotnie dłuższym krokiem między sąsiednimi węzłami. Z tym, że centralnego wzoru prostokątów nie używałeś w swoim przykładzie.

Na przykład obliczenia w 5 krokach na przedziale [0,1] dają:
trapezy: St = 0,2 * (f(0,0)/2 + f(0,2) + f(0,4) + f(0,6) + f(0,8) + f(1,0)/2)
centralne prostokąty: Sp = 0,2 * (f(0,1) + f(0,3) + f(0,5) + f(0,7) + f(0,9))
simpson: Ss = 0,2 * (f(0,0) + 4*f(0,1) + 2*f(0,2) + 4*f(0,3) + 2*f(0,4) + 4*f(0,5) + 2*f(0,6) + 4*f(0,7) + 2*f(0,8) + 4*f(0,9) + f(1,0)) / 6 =
= 0,2 * (f(0,0)/2 + f(0,2) + f(0,4) + f(0,6) + f(0,8) + f(1,0)/2) * 2/6 +
+ 0,2 * (f(0,1) + f(0,3) + f(0,5) + f(0,7) + f(0,9)) * 4/6 =
= (St + 2*Sp) / 3
Ten związek ma charakter uniwersalny, zachodzi dla każdej funkcji i dla każdej liczby węzłów.

Metoda 3/8 w przypadku 10 kroków jednostkowych jest niewykonalna, bo trzeba je blokować nie po 2, tylko po 3. Popatrz na ostatnią formułę w F11: sięga daleko za wyniki tablicowania Twojej funkcji. Jeżeli chcesz robić test porównawczy na wspólnej siatce węzłów, niech ich będzie 12, a nie 10.
Drugi błąd, to wielkość kroku zastosowana we wzorze 3/8. Implementujesz wersję, w której dx jest długością jednostkowej działki (np. a3-a2), ale w formule wstawiasz długość całego podprzedziału na którym wyliczany jest składnik sumy (np. a5-a2), czyli 3 razy za dużo.

Po korekcie obu tych błędów wyniki stają się poprawne.

Metody Simpsona 1/3 i 3/8 powinny dawać wyniki dokładne dla wielomianów do 3 stopnia włącznie. Co w przypadku 3/8 nie dziwi (interpolacja wielomianem 3. stopnia), a w przypadku 1/3 jest bardzo dobrym wynikiem -- błąd ma rząd o 1 mniejszy niż wynika to z samego stopnia interpolacji.
Ta sama właściwość, ale o dwa stopnie niżej, dotyczy także metody centralnych prostokątów. Systematyczna analiza utrzymania dokładności przybliżenia całki przy korzystaniu z możliwie małej liczby wywołań funkcji prowadzi do tzw. kwadratur Gaussa. Zaś skorzystanie z zależności między przybliżeniami różnych rzędów dla wyliczenia przybliżonej wartości granicznej znane jest jako metoda Romberga.
JJ
LO (7.6|24.2) ∙ Python (3.12|3.10) ∙ Unicode 15 ∙ LᴬTEX 2ε ∙ XML ∙ Unix tools ∙ Linux (Rocky|CentOS)
ODPOWIEDZ