obliczanie całek numerycznie

Mamy metodę, która wylicza całki z zadanej funkcji, np. takie:
exp(sin(x)), sin(x)/x, itd.

przykładowe wyniki dla sinx/x, od 0 do pi - kolejne kroki:

2: 1.847490777552169, er = -0.153
3: 1.852007377958264, er = 0.00452
4: 1.851936351344189, er = -7.1e-05
5: 1.851937056818966, er = 7.05e-07
6: 1.851937051957923, er = -4.86e-09
7: 1.851937051982544, er = 2.46e-11
8: 1.851937051982497, er = -4.71e-14
9: 1.851937051982458, er = -3.89e-14
10: 1.851937051982440, er = -1.78e-14
11: 1.851937051982448, er = 7.55e-15
12: 1.851937051983066, er = 6.18e-13
13: 1.851937051980151, er = -2.91e-12

gdzie: er = różnica pomiędzy kolejnymi krokami
zatem prawdopodobnie 11 daje najlepszy wynik:
1.851937051982448,

a teraz druga próba - z krokiem 1/2, co chyba powinno dać lepszy wynik:
2: 1.851695958423204, er = -0.0339
3: 1.851937974460681, er = 0.000242
4: 1.851937049740691, er = -9.25e-07
5: 1.851937051986264, er = 2.25e-09
6: 1.851937051982457, er = -3.81e-12
7: 1.851937051982470, er = 1.24e-14
8: 1.851937051982428, er = -4.13e-14
9: 1.851937051982596, er = 1.68e-13
10: 1.851937051982235, er = -3.61e-13

teraz mamy najlepszy w kroku 7: 1.851937051982470,
a poprzednio było w 11: 1.851937051982448;
różnica: 1.851937051982470 - 1.851937051982448 = 0.000000000000022 = 2.2e-14

ale który jest lepszy?

poprawny wynik: ?

......

int exp(sin(x)^2/x^2) dx; 0..pi =
2: 5.263772126141796, er = 1.4e-05
3: 5.263769320925950, er = -2.81e-06
4: 5.263769321647102, er = 7.21e-10
5: 5.263769321648034, er = 9.33e-13
6: 5.263769321648027, er = -7.11e-15
7: 5.263769321647993, er = -3.38e-14
8: 5.263769321647951, er = -4.26e-14
9: 5.263769321648339, er = 3.88e-13

?

Jaką metodą liczysz tę całkę?

Brakuje - 1 w 13-tej linii.

no to może takie coś:

exp(-x-exp(-x^2)); 0 do 1:

2: 0.288753667797922, er = 0.000195
3: 0.288753841408115, er = 1.74e-07
4: 0.288753841430690, er = 2.26e-11
5: 0.288753841430665, er = -2.48e-14
6: 0.288753841430665, er = -7.22e-16
7: 0.288753841430665, er = 7.77e-16
8: 0.288753841430661, er = -4.33e-15
9: 0.288753841430676, er = 1.47e-14
10: 0.288753841430662, er = -1.37e-14
11: 0.288753841430613, er = -4.9e-14
12: 0.288753841430808, er = 1.95e-13

Simpson dla n= 999 przedziałów daje: 0.288753841430666

który wynik jest poprawny? :)

Obliczasz przybliżenie całki i pytasz się który wynik jest poprawny...

Jak chcesz absolutny wynik to musisz policzyć ją analitycznie.

Analitycznie to, np.: f = x^5 możesz sobie obliczyć.

a teraz oblicz taki skecz:
y' = -2x y;
y = ?

Potrzeba mi lepszej precyzji: 30 cyfr chociaż!

Ten skecz, to gdzie ma puentę?
y'=-2xy
dy/dx = -2xy
dy/y=-2xdx
(lny)' = -2xdx
y= exp(-x^2) + C ?

@kwalifika: zdajesz sobie sprawę, że WolframAlpha jest nadal za darmo?
I że większość problemów analitycznych na poziomie liceum można tam rozwiązać?
https://www.wolframalpha.com/input/?i=y%27+%3D+-2x+y

Co tam wolfram... zabawka dla uczniów.

Błąd zaokrągleń rośnie chyba proporcjonalnie do: sqrt(n), n - liczba operacji.

dla Simpsona n = 1000 daje 32-50 eps,
a ten drugi błąd jest typu: h^4/80, zatem dla n = 1000, mamy tu ~10^-14,
co jest tym samym dla eps = 1e-16, czyli taki jest limit dokładności: n = 1000, powyżej będzie gorzej - większy błąd zamiast mniejszy!

Porównajmy to z tą drugą metodą: tam jest n^2 operacji na krok, ale rząd rośnie w każdym kroku = 2n;
Zatem z obcięcia dla n = 8 otrzymamy aż h^16, co dla 10-ciu przedziałów osiąga już limit = 0 w zakresie double;
a wtedy błąd zaokrągleń wynosi jakieś: 4n x eps = 4 x 8 eps = 32 eps.

OK.
Gdyby to obliczać na 30 cyfrach, wówczas Simpson potrzebuje aż: n = 1000000 - cały milion!
Natomiast 'Richardsona' n = 20 byłoby limitem, co daje z 500-1000 operacji, zaledwie.

Simpson obliczałby to aż 1000 razy dłużej!!

.........

Reasumując.
mamy tu zależności nakładów obliczeń zależnie od żądanej precyzji - liczby cyfr c:
P. 2^c -> tak działa Simpson, i chyba dowolna metoda stałego rzędu!
K. c^2 -> tak działa ekstrapolacja wg optymalnego schematu, czyli: 1, 2, 3, 4 ..., tu rząd rośnie liniowo: 2n.

przy podwojeniu precyzji otrzymamy wzrost liczby obliczeń:
P. 2^(2c) = 4^c = (2^c)^2 = n^2, strasznie dużo obliczeń!
K. (2c)^2 = 4c^2, tu mamy jedynie 4 krotnie więcej obliczeń

c = 16 (cyfr precyzji) ->
P (od prymityw) daje: 2^16 = 65000+ !
K : 16^2 = 256

widać różnicę?
trochę słabo...

c = 30 (cyfr precyzji) ->
P: 2^30 = 1000000000+ ponad miliard
K : 30^2 = 900, słownie: dziewięćset! :)

teraz lepiej?

Pozostaje zagadka: czy istniej szybsza metoda, od tej K = made by Kwalifika. :)

@kwalifika: W sumie przykro się patrzy, że masz 125 postów, i zaledwie 10 punktów reputacji. Czyżby nikt nie doceniał geniuszu?

Trzeba było się pojawić w piątek, akurat dyskutowaliśmy liczne ważkie sprawy - myślę że jeszcze ze dwa głębsze i już byśmy szukali dokładnych rozwiązań przybliżeniami.

enedil napisał(a):

@kwalifika: W sumie przykro się patrzy, że masz 125 postów, i zaledwie 10 punktów reputacji. Czyżby nikt nie doceniał geniuszu?

OP dawniej udzielał się jako @wil (koleś od oprogramowania do generowania krzyżówek, charakteryzujący się wręcz niebywałą ostrością umysłu), ale jakieś dwa lata temu za swoje zachowanie dostał perm bana. :D

Poznałem to po tym wątku, bo dawniej jegomość masowo tworzył tego typu tematy, w których popisywał się swoją wiedzą (jakiej jakości wiedzą to już pozostawiam wam do oceny), a każdego innego użytkownika traktował jak idiotę, który nie zna się na niczym. Często kończyło się to wycieczkami osobistymi i raportowaniem jego wypocin, aż miarka się przebrała.

No, można się rozejść. Sugeruję jak najszybciej zapomnieć o tym wątku i użytkowniku, bo szkoda na to czasu.