Intégration numérique d'une fonction compacte sur un triangle

comme le titre le suggère, j'essaie de calculer l'intégrale d'une fonction compacte (le polynôme quintique de Wendland) sur un triangle. Notez que le centre de la fonction se trouve quelque part dans l'espace 3D. J'intègre cette fonction sur un triangle arbitraire mais petit ( ). J'utilise actuellement l'intégration décrite par Dunavant, 1985 (p = 19).$area < \frac{(radius/4)^2}{2}$

Il semble cependant que ces règles de quadrature ne conviennent pas aux problèmes pris en charge de manière compacte. Ceci est soutenu par le fait que lorsque (donc une fonction qui est 1 à l'intérieur du cercle de rayon 1) sur un plan qui est discrétisé à l'aide de triangles, mes résultats (normalisés) se situent entre 1,001 et 0,897.$f(r) = [r\leq1]$

Ma question est donc la suivante: existe-t-il une règle de quadrature spécialisée pour ce type de problème? Une règle d'intégration composite d'ordre inférieur fonctionnerait-elle mieux?

Malheureusement, cette routine est vraiment critique dans mon code, donc la précision est cruciale. D'un autre côté, j'ai besoin de faire cette intégration "quelques fois" pour un seul pas de temps afin que les dépenses de calcul ne soient pas trop élevées. La parallélisation n'est pas un problème car j'exécuterai l'intégration elle-même en série.

Merci d'avance pour vos réponses.

EDIT: le polynôme quintique de Wendland est donné par avec et avec r_0 étant un vecteur arbitraire dans \ mathbb {R} ^ 3α=$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$r_0$$\mathbb{R}^3$

EDIT2: Si $\Delta$ est le triangle à deux dimensions, alors je veux calculer $\int_\Delta \omega(r) dr$ avec $\omega(r) = W(\frac{\|r-r_0\|}{h})$ . Donc $q$ dans $W$ ne sera jamais inférieur à 0. Notez que l'intégrale est une intégrale de surface sur une surface 2D dans $\mathbb{R}^3$

EDIT3: J'ai une solution analytique pour le problème 1-D (ligne). Le calcul d'un pour 2-D (triangle) pourrait également être possible.

Puisque la fonction est lisse dans , mais pas de degré fixe (dans le plan, c'est-à-dire), je suggérerais d'utiliser un schéma adaptatif simple, par exemple la règle trapézoïdale avec la méthode de Romberg , dans les deux dimensions.$q \leq 2$

Autrement dit, si votre triangle est défini par les sommets , et , et que vous avez une routine qui s'intègre le long de la ligne de à , vous pouvez faire ce qui suit (en notation Matlab):$x$z ∈ R$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 );

Dans romb, n'utilisez pas un nombre fixe de points, mais continuez à agrandir le tableau jusqu'à ce que la différence entre deux diagonales successives soit inférieure à votre tolérance requise. Puisque votre fonction est fluide, cela devrait être une bonne estimation d'erreur.

Si des parties du triangle sont en dehors du domaine de , vous pouvez essayer d'ajuster les limites d'intégration dans le code ci-dessus en conséquence.$W(q)$

Ce n'est peut-être pas le moyen le plus efficace sur le plan informatique de résoudre votre problème, mais l'adaptabilité vous donnera beaucoup plus de robustesse qu'une règle à degrés fixes.

Pour une bonne vue d'ensemble des règles de cubature, voir "R. Cools, An Encyclopédie de Cubature Formulas J. Complexity, 19: 445-453, 2003". L'utilisation d'une règle fixe peut vous donner l'avantage que certaines règles intègrent exactement les polynômes (comme le fait la quadrature gaussienne en une dimension).

Cools est également l'un des principaux auteurs de CUBPACK , un progiciel de cubature numérique.

Les règles d'intégration supposent que la fonction est localement bien approximée par un polynôme de faible degré. Votre problème n'a rien à voir avec un support compact. Les fonctions de base radiales compactes sont lisses à la limite du support et les règles de quadrature jusqu'à l'ordre de lissage peuvent être utilisées sans problème. (Les règles d'ordre supérieur n'aident pas; vous ne devriez donc probablement pas utiliser une règle qui intègre exactement les polynômes de degré 5.)

Dans votre cas, l'inexactitude vient du fait que l'hypothèse d'une bonne approximabilité polynomiale échoue dans votre cas pour les triangles proches de , même lorsqu'ils ne contiennent pas .r $r_0$$r_0$

q q r r → r 0 r $W$ est lisse en fonction de , mais est une fonction non lisse de , avec un gradient qui devient infini dans la limite . L'intégration est sur , et la fonction composite est une fonction non lisse de .$q$$q$$r$$r\rightarrow r_0$$r$$r$

Si le triangle ne contient pas , la fonction est mais cela n'aide pas car la dérivée supérieure croît très rapidement près de , et l'erreur d'une méthode d'ordre élevé est proportionnelle à une dérivée d'ordre élevé, donc très grande !C i n f r $r_0$$C^inf$$r_0$

Le remède simple consiste à diviser chaque triangle T en un nombre N_T de sous-triangles. Vous pouvez prendre loin de , et près de . Vous pouvez déterminer hors ligne la taille de pour que les triangles d'un diamètre donné et la distance de atteignent la précision souhaitée. De plus, vous ne devez utiliser que des formules d'ordre bas proches de .r 0 N T ≫ 1 r 0 N T r 0 r $N_T=1$$r_0$$N_T\gg 1$$r_0$$N_T$$r_0$$r_0$

Lorsque vous vous sur un triangle, mais que est tridimensionnel, le triangle est apparemment dans .R $r_0$$R^3$

Un remède plus rapide consisterait donc à tabuler l'intégrale pour en fonction des coordonnées du triangle (normalisé en le faisant pivoter dans un plan bidimensionnel de telle sorte qu'un sommet se trouve sur l' axe des , et en le reflétant de telle sorte qu'une seconde sommet se trouve au-dessus). Cette tabulation doit être suffisamment détaillée pour rendre une interpolation linéaire ou quadratique suffisamment précise. Mais vous pouvez utiliser la méthode lente décrite en premier pour créer ce tableau.x y $r_0=0$$xy$$x$

Une autre façon de se débarrasser du problème consiste à utiliser une fonction de base radiale compacte qui est un polynôme dans plutôt que . C'est fluide partout et facile à intégrer.$q^2$$q$