Numeryczna integracja kompaktowo obsługiwanej funkcji na trójkącie

10

jak sugeruje tytuł, próbuję obliczyć całkę kompaktowo obsługiwanej funkcji (kwintyczny wielomian Wendlanda) na trójkącie. Zauważ, że środek funkcji znajduje się gdzieś w przestrzeni 3D. Integruję tę funkcję na dowolnym, ale małym trójkącie ( area<(radius/4)22 ). Obecnie używam integracji opisanej przez Dunavant, 1985 (p = 19).

Wydaje się jednak, że te reguły kwadraturowe nie są odpowiednie dla problemów wspieranych w sposób kompaktowy. Potwierdza to fakt, że gdy całkuję f(r)=[r1] (czyli funkcja 1 w okręgu o promieniu 1) na płaszczyźnie, która jest dyskretyzowana za pomocą trójkątów, moje (znormalizowane) wyniki są pomiędzy 1,001 i 0,897.

Więc moje pytanie brzmi: czy istnieje specjalistyczna reguła kwadraturowa dla tego rodzaju problemu? Czy lepiej byłoby zastosować regułę integracji złożonej niższego rzędu?

Niestety ta procedura jest bardzo ważna w moim kodzie, więc precyzja jest kluczowa. Z drugiej strony muszę wykonać tę integrację „kilka razy” dla jednego kroku czasowego, aby koszt obliczeniowy nie był zbyt wysoki. Równoległość nie stanowi problemu, ponieważ sama integracja zostanie wykonana szeregowo.

Dzięki z góry za odpowiedzi.

EDYCJA: Kwintyczny wielomian Wendlanda daje W(q)=[q2]αh3(1q2)4(2q+1)przyα=2116π r0R3q=rr0hr0R3

EDIT2: Jeśli Δ jest dwuwymiarowym trójkątem, to chcę obliczyć Δω(r)dr z ω(r)=W(rr0h) . Więc q w W nigdy nie będzie mniejsze niż 0. Zauważ, że całka jest całką powierzchniową nad powierzchnią 2-D w R3

EDYCJA 3: Mam rozwiązanie analityczne dla problemu 1-D (linia). Możliwe jest również obliczenie jednego dla 2-D (trójkąta).

Azrael3000
źródło
Czy możesz podać nam kilka dodatkowych szczegółów na temat funkcji, którą próbujesz zintegrować? Czy to tylko wielomian? A może wielomian częściowy?
Pedro
Edytowane zgodnie z żądaniem.
Azrael3000

Odpowiedzi:

4

Ponieważ funkcja jest gładka wewnątrz , ale nie od ustalonego stopnia (w płaszczyźnie, która jest), chciałbym zaproponować za pomocą prostego systemu adaptacyjnego, np trapezów z metoda romberga w obu wymiarach.q2

To znaczy, jeśli twój trójkąt jest zdefiniowany przez wierzchołki , i , a masz procedurę, która integruje się wzdłuż linii od do , możesz wykonać następujące czynności (w notacji Matlab):y z R 3xyzR3romb(f,a,b)fab

int = romb( @(xi) romb( W , xi , y+(z-y)*(xi-x)./(z-x) ) , x , z );

W romb, nie używaj stałej liczby punktów, ale powiększaj tabelę, dopóki różnica między dwoma kolejnymi przekątnymi nie spadnie poniżej wymaganej tolerancji. Ponieważ twoja funkcja jest płynna, powinno to być dobre oszacowanie błędu.

Jeśli części trójkąta znajdują się poza domeną , możesz spróbować dostosować granice integracji w powyższym kodzie.W(q)

To może nie być najskuteczniejszy pod względem obliczeniowym sposób rozwiązania twojego problemu, ale adaptacja da ci znacznie większą niezawodność niż reguła o stałym stopniu.

Pedro
źródło
Ta funkcja jest wszędzie, z wyjątkiem . Sąsiedztwo tego punktu powoduje problemy. q=0
Arnold Neumaier
Ach, rozkładając się na dwa problemy 1-D, wcale nie jest to zły pomysł. Ponieważ jest jedna rzecz, której ci nie powiedziałem. Mam rozwiązanie analityczne w 1-D, więc mogę zastąpić wewnętrzny romb funkcją analityczną. Dam temu szansę już +1
Azrael3000
@ArnoldNeumaier, przepraszam, nie rozumiem, jak to możliwe. Czy możesz wytłumaczyć?
Pedro
gładka jako funkcja , ale jest niepłynną funkcją , a całkowanie jest ponad , o ile rozumiem pytanie. Funkcja złożona jest zatem funkcją gładką . q r r rqqrrr
Arnold Neumaier
1
@Pedro Wdrożyłem go i działa jak urok. Znaleźliśmy dziś również rozwiązanie analityczne. Ale dotyczy to tylko specjalnego przypadku, który można wykorzystać do odtworzenia ogólnego. Oznacza to, że musimy dokonać dekompozycji domen. Ponieważ Romberg zbiega się w około 4 krokach, myślę, że z tego powodu będzie on szybszy niż stosowanie formuły analitycznej. Według Wikipedii możemy robić jeszcze lepiej niż Romberg, stosując racjonalne wielomiany. Znajdziesz swoje imię w podziękowaniach za mój następny artykuł :) Pozdrawiam.
Azrael3000
2

Aby uzyskać dobry przegląd reguł kubatury, patrz „R. Cools, An Encyclopaedia of Cubature Formulas J. Complexity, 19: 445-453, 2003”. Korzystanie ze stałej reguły może dać przewagę polegającą na tym, że niektóre reguły dokładnie integrują wielomiany (tak jak kwadratura Gaussa w jednym wymiarze).

Cools jest także jednym z głównych autorów CUBPACK , pakietu oprogramowania do kubatury numerycznej.

GertVdE
źródło
Myślę, że problem polega na tym, że funkcja jest wielomianem , ale jest funkcją nieliniową we współrzędnych przestrzennych. Funkcja jest gładka do krawędzi funkcji bazowej, ale nie jest wielomianowa, z wyjątkiem wzdłuż osi. qqq
Pedro
To prawda Pedro.
Azrael3000
ah ok. mój błąd. Przepraszam.
GertVdE
2

Reguły integracji zakładają, że funkcja jest lokalnie dobrze przybliżona przez wielomian niskiego stopnia. Twój problem nie ma nic wspólnego z kompaktową obsługą. Kompaktowo obsługiwane funkcje podstawy radialnej są gładkie na granicy podparcia, a zasady kwadratury do rzędu gładkości mogą być używane bez problemów. (Reguły wyższego rzędu nie pomagają, dlatego prawdopodobnie nie powinieneś używać reguły, która integruje dokładnie wielomiany stopnia 5).

W twoim przypadku niedokładność wynika z faktu, że założenie dobrej przybliżalności wielomianowej zawodzi w twoim przypadku dla trójkątów w pobliżu , nawet jeśli nie zawierają one .r 0r0r0

q q r r r 0 r rW jest gładkie w funkcji , ale jest niepłynną funkcją , z gradientem, który staje się nieskończony w granicy prawo w prawo . Całkowanie jest ponad , a funkcja złożona jest funkcją nie gładką .qqrrr0rr

Jeśli trójkąt nie zawiera , funkcją jest ale to nie pomaga, ponieważ wyższa pochodna rośnie bardzo szybko blisko , a błąd metody wysokiego rzędu jest proporcjonalny do pochodnej wysokiego rzędu, a zatem bardzo duży !C i n f r 0r0Cinfr0

Prostym lekarstwem jest podzielenie każdego trójkąta T na liczbę N_T pod trójkątów. Możesz wziąć daleko od , a blisko . Możesz dowiedzieć się offline, jak duże musi być dla trójkątów o danej średnicy i odległości od aby osiągnąć pożądaną dokładność. Ponadto powinieneś używać formuł niskiego rzędu w pobliżu .r 0 N T1 r 0 N T r 0 r 0NT=1r0NT1r0NTr0r0

Gdy całkujesz przez trójkąt, ale jest trójwymiarowy, trójkąt najwyraźniej znajduje się w .R 3r0R3

W związku z tym szybszy środek zaradczy przedstawiłby całkę dla jako funkcję współrzędnych trójkąta (znormalizowanych przez obrócenie jej w 2-wymiarową płaszczyznę taki sposób, że jeden wierzchołek leży na osi i odzwierciedlenie jej w taki sposób, że druga wierzchołek leży nad nim). Ta tabela musi być wystarczająco szczegółowa, aby interpolacja liniowa lub kwadratowa była wystarczająco dokładna. Ale możesz użyć powolnej metody opisanej w pierwszej kolejności, aby utworzyć tę tabelę.x y xr0=0xyx

Innym sposobem na pozbycie się problemu jest użycie kompaktowo obsługiwanej funkcji podstawy radialnej, która jest wielomianem w zamiast . Jest to płynne wszędzie i łatwe do zintegrowania. qq2q

Arnold Neumaier
źródło
Myślę, że jest małe nieporozumienie. Zaktualizowałem opis mojego pytania. W rzeczywistości w całce nigdy nie może być mniejsza niż 0. I niekoniecznie jest zawarte w trójkącie. r 0qr0
Azrael3000
Twój nowy dodatek nie ma dla mnie sensu. Jeśli to musi być . Czy integrujesz ponad trójkątem 2D w ? - Nie zakładałem, że jest w trójkącie. Właśnie dodałem za chwilę więcej szczegółów do mojej odpowiedzi. r R 3 r 0r0R3rR3r0
Arnold Neumaier
Tak, to prawda, że ​​całkuję na trójkącie 2D w . R3
Azrael3000