czy mogę ufać tej liczbowej potrójnej całce z Matlaba?

Ludzie informatyki:

I pierwotnie opublikowany na to pytanie w matematyce Stos Exchange i ktoś powiedział, że mogę dostać się „znacznie lepiej” odpowiedź tutaj:

Jestem nowicjuszem w metodach numerycznych i Matlabie. Próbuję oszacować następującą sumę dwóch potrójnych całek (można to oczywiście napisać prościej, ale nadal nie można jej ocenić symbolicznie (?)). Mam problem z uzyskaniem $\LaTeX$ do pracy tutaj, więc niechętnie podzieliłem go tutaj na części: chcę znaleźć sumę

\frac{2)}{((1 / 0,3) - 1)^{2)}} (\int_{1}^{1 / 0,3} \int_{1}^{r_{1}} \int_{0}^{r_{1} - r_{0}} {fa}_{1} (r_{0}, r_{1}, t) \exp (- \frac{(0,3)^{2)} t^{2)}}{4}) re t re r_{0} re r_{1}),

$\frac{2}{((1/0.3) - 1)^2}\left(\int_1^{1/0.3}\int_1^{r_1}\int_0^{r_1-r_0}F_1(r_0,r_1,t)\exp(-\frac{(0.3)^2 t^2}{4})\,dt\,dr_0\,dr_1 \right),$

\frac{2)}{((1 / 0,3) - 1)^{2)}} (\int_{1}^{1 / 0,3} \int_{1}^{r_{1}} \int_{r_{1} - r_{0}}^{r_{1} + r_{0}} {fa}_{2)} (r_{0}, r_{1}, t) \exp (- \frac{(0,3)^{2)} t^{2)}}{4}) re t re r_{0} re r_{1}),

$\frac{2}{((1/0.3) - 1)^2}\left(\int_1^{1/0.3}\int_1^{r_1}\int_{r_1-r_0}^{r_1+r_0} F_2(r_0,r_1,t)\exp(-\frac{(0.3)^2 t^2}{4})\,dt\,dr_0\,dr_1 \right),$

gdzie

F_{1} (r_{0}, r_{1}, t) = \frac{t^{2} r_{0}^{3} * (0.3)^{3}}{2 r_{1}^{3} \sqrt{π}}

$F_1(r_0,r_1,t)=\frac{t^2 r_0^3*(0.3)^3}{2r_1^3\sqrt{\pi}}$

F_{2} (r_{0}, r_{1}, t) = \frac{(0.3)^{3} π^{3 / 2} (r_{0} + r_{1} - t)^{4} (t^{2} + 2 t (r_{0} + r_{1}) - 3 (r_{1} - r_{0})^{2})^{2}}{288 (\frac{4}{3} π r_{0}^{3}) (\frac{4}{3} π r_{1}^{3})} .

$F_2(r_0,r_1,t)=\frac{(0.3)^3\pi^{3/2}(r_0+r_1-t)^4 (t^2+2t(r_0+r_1)-3(r_1-r_0)^2)^2}{288(\frac{4}{3}\pi r_0^3)(\frac{4}{3}\pi r_1^3)}.$

EDYCJA (2 marca 2013 r.): Ktoś odpowiedział, że zmusił Mathematica do wykonania całek symbolicznie. Właśnie próbowałem to zrobić (z uproszczonymi wersjami całek), a Mathematica mogła wykonać tylko dwie zewnętrzne z pierwszej, i zatrzymałem się na drugiej. Byłbym wdzięczny za pomoc. Oto co zrobiłem:

Próbowałem ocenić

przez

\int_{1}^{2)} \int_{1}^{r_{2)}} \int_{0}^{r_{2)} - r_{1}} \frac{r_{1}^{3)} t^{2)} \exp (- t^{2)})}{r_{2)}^{3)}} re t re r_{1} re r_{2)}

$\int_1^2 \int_1^{r_2} \int_0^{r_2-r_1} \frac{r_1^3 t^2 \exp(-t^2)}{r_2^3}\,dt\,dr_1\,dr_2$

Zintegruj [r1 ^ 3 / r2 ^ 3 * t ^ 2 * Exp (-t ^ 2), {t, 0, r2 - r1}, {r1, 1, r2}, {r2, 1, 2}]

i Mathematica powraca (miałem problem z tutaj, ponieważ wynik jest długi. Zerwałem na dwa równania. jeśli ktoś zna dobry sposób na wyświetlenie tego, proszę mi powiedzieć): $\LaTeX$

\int_{1}^{2)} \frac{1}{64 r {2)}^{2)}} {mi}^{- 1 - r {2)}^{2)}} (2) {mi}^{2) r 2)} (25 + r 2) (19 + 2) r 2) (1 + r 2)))) -

$\int_1^2 \frac{1}{64r2^2} e^{-1-r2^2}(2e^{2r2}(25+r2(19+2r2(1+r2)))-$

{mi}^{1 + r {2)}^{2)}} (32 r 2) (2) + r {2)}^{2)})) + \sqrt{π} (11 + 4 r {2)}^{2)} (9 + r {2)}^{2)})) Erf [1 - r 2)]) re r 2)

$e^{1+r2^2}(32r2(2+r2^2)) +\sqrt{\pi}(11+4r2^2(9+r2^2))\operatorname{Erf}[1-r2])\,dr2.$

Potem próbowałem ocenić

\int_{1}^{2)} \int_{1}^{r_{2)}} \int_{r_{2)} - r_{1}}^{r_{2)} + r_{1}} \dots

$\int_1^2\int_1^{r_2}\int_{r_2-r_1}^{r_2+r_1} \ldots \qquad \qquad \qquad$

\dots \frac{\exp (- t^{2)}) (r_{1} + r_{2)} - t)^{4} (t^{2)} + 2) t (r_{1} + r_{2)}) - 3) (r_{2)} - r_{1})^{2)})^{2)}}{r_{1}^{3)} r_{2)}^{3)}} re t re r_{} re r_{2)}

$\ldots\frac{\exp(-t^2)(r_1+r_2-t)^4(t^2+2t(r_1+r_2)-3(r_2-r_1)^2)^2}{r_1^3 r_2^3}\,dt\,dr_\,dr_2$

za pomocą

Zintegruj [(r1 + r2 - t) ^ 4 * (t ^ 2 + 2 * t * (r1 + r2) - 3 * (r2 - r1) ^ 2) ^ 2 * Exp [-t ^ 2] / r1 ^ 3 / r2 ^ 3, {r2, 1, 2}, {r1, 1, r2}, {t, r2-r1, r2 + r1}]

właśnie teraz, a Mathematica nie odpowiedziała po upływie około pół godziny (ale mam teraz problemy z siecią komputerową i mogą być winni).

[KONIEC EDYCJI 2 MARCA]

$0.007164820144202$

Wiem, że Matlab jest dobrym oprogramowaniem, ale słyszałem, że trudne są dokładne potrójne całki numeryczne, a matematycy powinni być sceptyczni, dlatego chcę w jakiś sposób zweryfikować dokładność tej odpowiedzi. Całki podają oczekiwaną wartość określonego eksperymentu (jeśli ktoś chce, mogę edytować to pytanie, aby opisać eksperyment): zaimplementowałem eksperyment w Matlabie, używając odpowiednio losowo wygenerowanych liczb, milion razy, i uśredniłem wyniki. Powtórzyłem ten proces cztery razy. Oto wyniki (przepraszam, jeśli niewłaściwie użyłem słowa „rozprawa”):

$0.007133292603256$

$0.007120455071989$

Trial 3: $0.007062595022049$

Trial 4: $0.007154940168452$

Trial 5: $0.007215000289130$

Chociaż w każdej próbie wykorzystano milion próbek, wartości symulacji zgadzają się tylko w pierwszej znaczącej cyfrze. Nie są wystarczająco blisko siebie, aby ustalić, czy potrójna całka numeryczna jest dokładna.

Czy ktoś może mi więc powiedzieć, czy mogę zaufać wynikowi „potrójnego kwadratu” i pod jakimi warunkami ogólnie można mu ufać?

Jedną z sugestii, które otrzymałem na Math Stack Exchange, było wypróbowanie innego oprogramowania, takiego jak Mathematica, Octave, Maple i SciPy. Czy to dobra rada? Czy ludzie faktycznie wykonują prace numeryczne w Mathematica i Maple? Octave jest rodzajem klonu Matlaba, więc czy mogę założyć, że używa tych samych algorytmów integracji? Jeszcze nie słyszałem o SciPy i doceniłbym wszelkie opinie na jego temat.

AKTUALIZACJA: Ktoś z Math Stack Exchange zrobił to w Maple i dostał $0.007163085468$ . Jest to zgoda na trzy znaczące liczby. To dobry znak.

Byłbym także wdzięczny za sugestie, jak wprowadzić długie, wieloliniowe wyrażenie $\LaTeX$ w stosie wymiany. Czy możesz tutaj użyć środowiska „wyrównanego”? Próbowałem i nie mogłem tego uruchomić.

matlab

— Stefan Smith
źródło

Twoje wyniki symulacji są całkowicie zgodne z wartością liczbową zwracaną przez Matlab: ich średnią z

0.00713726

$0.00713726$ jest tylko

- 1.11

$-1.11$ standardowe błędy mniejsze niż zwracane przez Matlaba. FWIW, Mathematica powraca

0.00716308537 \dots

$0.00716308537\ldots$ . Może także symbolicznie oceniać te całki pod względem wielomianów i funkcji błędów.

— whuber

@whuber Thanks. Mógłbym przysiąc, że próbowałem symbolicznie w Maple i Maple nie mógł tego zrobić. Spróbuję ponownie w Maple, a jeśli to nie zadziała, spróbuję w Mathematica. BTW, zrobiłem podobną całkę w Maple i dostałem ogromną symboliczną odpowiedź. Wydawało się, że jest to suma i różnica bardzo dużych liczb, których suma była dość mała. Podejrzewam, że błąd zaokrąglenia był prawdopodobny w ostatecznej odpowiedzi. Czy w takim problemie należy użyć symbolicznej odpowiedzi, czy po prostu wykonać całkę numerycznie?

— Stefan Smith

Zalety odpowiedzi symbolicznych stanowią kombinacje funkcji, które (często) można skutecznie obliczyć z dowolną dokładnością. Zwykle również symboliczne rozwiązanie nadaje się do szybkiego ponownego obliczenia, gdy parametry są zmienne. Z tych powodów często warto szukać symbolicznego rozwiązania.

— whuber

@ whuber: Próbowałem wykonać kilka zasadniczo równoważnych całek (zmieniając niektóre stałe i usuwając niektóre stałe multiplikatywne) w Mathematica, a Mathematica mogła wykonać tylko dwie zewnętrzne całki pierwszej całki i wydaje się, że utknęła na drugiej. Kod i wyniki opublikowałem powyżej.

— Stefan Smith

Re Edycja z 2 marca: Poprzez zredukowanie potrójnej całki symbolicznie do pojedynczej całki (w pierwszej połowie całek) osiągnąłeś wiele. Integrand zachowuje się bardzo ładnie i może być zintegrowany numerycznie z bardzo wysoką precyzją w ułamku sekundy.

— whuber

Przede wszystkim to nie oprogramowanie (a przynajmniej tak nie powinno być) decyduje o jakości rozwiązania problemu, lecz o jakości i stosowności zastosowanego algorytmu. Powinieneś sprawdzić, jakiego algorytmu używa Matlab w potrójnym kwadracie (domyślam się, że używa zagnieżdżonej adaptacyjnej kwadratury Gaussa). I powinieneś sprawdzić, jakie są wymagane tolerancje (wymagana tolerancja bezwzględna i względna). Są szanse, że domyślnie tylko o to prosi $10^{-8}$ względna precyzja.

Odpowiedź pochodząca z Maple jest prawdopodobnie udzielona przez Computer Algebra i być może mogłaby znaleźć rozwiązanie zamknięte, które następnie zostało ocenione przy użyciu zmiennoprzecinkowego podwójnej precyzji. Ma to tę zaletę, że nie przybliżasz całki za pomocą skończonego sumowania (a zatem wprowadzasz błędy aproksymacji), ale system algebry komputerowej znajdzie wyrażenie dla całki, które można następnie ocenić. Oczywiście należy zachować ostrożność przy ocenie tego wyrażenia (w celu zaokrąglenia).

Jeśli chcesz to zrobić za pomocą SciPy, musisz również skorzystać z zagnieżdżonej adaptacyjnej kwadratury gaussowskiej przy użyciu podstawowych procedur Quadpack (Piessens i in.). W Octave będziesz miał takie samo podejście. I nie byłbym zaskoczony, gdyby Matlab używał Quadpack jako silnika kwadraturowego (ponieważ jest to odniesienie).

— GertVdE
źródło

@GretVdE: Dzięki za informację. Próbowałem najpierw oszacować całkę symbolicznie, a Maple nie mógł tego zrobić (więc prawdopodobnie było to niemożliwe przy użyciu standardowych funkcji), więc poprosiłem Maple o zrobienie tego liczbowo. Nie wiem, jakiego algorytmu użył.

— Stefan Smith

@StefanSmith: można dowiedzieć się poprzez ustawienie infolevel w Maple: infolevel[`evalf/int`] := 4. Czy na pewno Mape nie może znaleźć rozwiązania zamkniętego? Całka nie wydaje się zbyt skomplikowana. Czy mógłbyś gdzieś upublicznić swój arkusz klonu?

— GertVdE

@StefanSmith: Chciałbym zamieścić kod klonu w powyższym pytaniu.

— GertVdE

Nie mogę teraz zmusić Mapla do pracy w moim systemie, ale próbowałem równoważnych całek w Mathematica, a Mathematica wykonała tylko wewnętrzne dwie pierwsze potrójne całki i utknęła na drugiej potrójnej całce. Zobacz edytowane pytanie.

— Stefan Smith