Numerische Integration der kompakt unterstützten Funktion in ein Dreieck

Wie der Titel schon sagt, versuche ich, das Integral einer kompakt unterstützten Funktion (Wendlands Quintpolynom) auf einem Dreieck zu berechnen. Beachten Sie, dass sich das Zentrum der Funktion irgendwo im 3D-Raum befindet. Ich integriere diese Funktion in ein beliebiges, aber kleines Dreieck ( ). Ich verwende derzeit die von Dunavant, 1985, beschriebene Integration (p = 19).$area < \frac{(radius/4)^2}{2}$

Es scheint jedoch, dass diese Quadraturregeln nicht für kompakt unterstützte Probleme geeignet sind. Dies wird durch die Tatsache unterstützt, dass meine (normalisierten) Ergebnisse dazwischen liegen , wenn ich (also eine Funktion, die 1 innerhalb des Kreises mit Radius 1 ist) in eine Ebene integriere, die mit Dreiecken diskretisiert wird 1,001 und 0,897.$f(r) = [r\leq1]$

Meine Frage ist also, ob es für diese Art von Problem eine spezielle Quadraturregel gibt. Würde eine zusammengesetzte Integrationsregel niedrigerer Ordnung besser funktionieren?

Leider ist diese Routine in meinem Code sehr wichtig, daher ist Präzision von entscheidender Bedeutung. Andererseits muss ich diese Integration "ein paar Mal" für einen einzelnen Zeitschritt durchführen, damit der Rechenaufwand nicht zu hoch ist. Parallelisierung ist kein Problem, da ich die Integration selbst seriell ausführen werde.

Vielen Dank im Voraus für Ihre Antworten.

EDIT: Wendlands quintisches Polynom ist gegeben durch mitα=$W(q) = [q\leq2]\frac{\alpha}{h^3}(1-\frac{q}{2})^4(2q+1)$$\alpha = \frac{21}{16\pi}$$q=\frac{\|r-r_0\|}{h}$$r_0$$\mathbb{R}^3$

EDIT2: Wenn das zweidimensionale Dreieck ist, dann möchte ich mit berechnen. . So in wird nie kleiner als 0. Hinweis sein , dass das Integral ein Oberflächenintegral über eine Oberfläche in 2-D$\Delta$$\int_\Delta \omega(r) dr$$\omega(r) = W(\frac{\|r-r_0\|}{h})$$q$$W$$\mathbb{R}^3$

EDIT3: Ich habe eine analytische Lösung für das 1-D (Line) -Problem. Möglicherweise ist auch eine Berechnung für 2-D (Dreieck) möglich.

Da die Funktion innerhalb von glatt ist , aber nicht von festem Grad (dh in der Ebene), würde ich vorschlagen , in beiden Dimensionen ein einfaches adaptives Schema zu verwenden, z. B. die Trapezregel nach Rombergs Methode .$q \leq 2$

Das heißt, wenn Ihr Dreieck durch die Eckpunkte , und ist und Sie eine Routine haben, die entlang der Linie von bis integriert ist , können Sie Folgendes tun (in Matlab-Notation):$x$$y$$z \in \mathbb R^3$romb(f,a,b)fab

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

In romb, keine feste Anzahl von Punkten, verwendet aber die Tabelle weiter wachsen , bis die Differenz zwischen zwei aufeinanderfolgenden Diagonalen unterhalb Ihrer geforderte Toleranz ist. Da Ihre Funktion reibungslos funktioniert, sollte dies eine gute Fehlerschätzung sein.

Wenn Teile des Dreiecks außerhalb der Domäne von , können Sie versuchen, die Integrationsgrenzen im obigen Code entsprechend anzupassen.$W(q)$

Dies ist möglicherweise nicht die rechnerisch effizienteste Methode zur Lösung Ihres Problems, aber die Adaptivität bietet Ihnen viel mehr Robustheit als eine Regel mit festem Grad.

Für einen guten Überblick über Kubaturregeln siehe "R. Cools, Eine Enzyklopädie der Kubaturformeln J. Complexity, 19: 445-453, 2003". Die Verwendung einer festen Regel kann Ihnen den Vorteil bieten, dass einige Regeln Polynome genau integrieren (wie dies die Gaußsche Quadratur in einer Dimension tut).

Cools ist auch einer der Hauptautoren von CUBPACK , einem Softwarepaket für numerische Kubatur.

Integrationsregeln setzen voraus, dass die Funktion lokal durch ein Polynom niedrigen Grades gut angenähert ist. Ihr Problem hat nichts mit kompakter Unterstützung zu tun. Kompakt unterstützte radiale Basisfunktionen sind an der Stützgrenze glatt, und Quadraturregeln bis zur Größenordnung der Glätte können problemlos verwendet werden. (Regeln höherer Ordnung helfen nicht; daher sollten Sie wahrscheinlich keine Regel verwenden, die Polynome des Grades 5 genau integriert.)

In Ihrem Fall ergibt sich die Ungenauigkeit aus der Tatsache, dass die Annahme einer guten Polynomnäherung in Ihrem Fall für Dreiecke in der Nähe von fehlschlägt , selbst wenn sie kein enthalten .$r_0$$r_0$

$W$ ist glatt als Funktion von , aber ist eine nicht glatte Funktion von mit einem Gradienten, der in der Grenze unendlich wird . Die Integration ist über und die zusammengesetzte Funktion ist eine nicht glatte Funktion von .$q$$q$$r$$r\rightarrow r_0$$r$$r$

Wenn das Dreieck nicht enthält , ist die Funktion aber dies hilft nicht, da die höhere Ableitung sehr schnell nahe an wächst und der Fehler einer Methode höherer Ordnung proportional zu einer Ableitung hoher Ordnung ist, daher sehr groß !$r_0$$C^inf$$r_0$

Die einfache Abhilfe besteht darin, jedes Dreieck T in eine Anzahl N_T von Unterdreiecken aufzuteilen. Sie können weit weg von und nahe . Sie können offline herausfinden, wie groß muss, damit Dreiecke mit einem bestimmten Durchmesser und Abstand von eine gewünschte Genauigkeit erreichen. Darüber hinaus sollten Sie nur Formeln niedriger Ordnung in der Nähe von .$N_T=1$$r_0$$N_T\gg 1$$r_0$$N_T$$r_0$$r_0$

Wenn Sie über ein Dreieck integrieren, aber dreidimensional ist, befindet sich das Dreieck anscheinend in .$r_0$$R^3$

Eine schnellere Abhilfe würde daher das Integral für als Funktion der Dreieckskoordinaten tabellieren (normalisiert durch Drehen in eine zweidimensionale Ebene, so dass ein Scheitelpunkt auf der Achse liegt, und so reflektieren, dass eine Sekunde Scheitelpunkt liegt darüber). Diese Tabelle muss ausreichend detailliert sein, um eine lineare oder quadratische Interpolation genau genug zu machen. Sie können jedoch die zuerst beschriebene langsame Methode verwenden, um diese Tabelle zu erstellen.$r_0=0$$xy$$x$

Eine andere Möglichkeit, das Problem zu beseitigen, besteht darin, eine kompakt unterstützte radiale Basisfunktion zu verwenden, die ein Polynom in anstelle von . Dies ist überall reibungslos und einfach zu integrieren.$q^2$$q$