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

10

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<(radius/4)22

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)=[r1]

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α=21W(q)=[q2]αh3(1q2)4(2q+1)α=2116πq=rr0hR 3r0R3

EDIT2: Si Δ est le triangle à deux dimensions, alors je veux calculer Δω(r)dr avec ω(r)=W(rr0h) . 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 R3

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.

Azrael3000
la source
Pourriez-vous nous donner quelques détails supplémentaires sur la fonction que vous essayez d'intégrer? Est-ce juste un polynôme? Ou un polynôme par morceaux?
Pedro
Modifié comme demandé.
Azrael3000

Réponses:

4

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.q2

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):xz R 3yzR3romb(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.

Pedro
la source
La fonction est smmoth partout sauf pour . Le voisinage de ce point pose problème. q=0
Arnold Neumaier
Ah se décomposant en deux problèmes 1-D, pas une mauvaise idée du tout. Parce qu'il y a une chose que je ne vous ai pas dit. J'ai une solution analytique en 1-D pour pouvoir remplacer le romb intérieur par une fonction analytique. Je donnerai déjà un coup +1
Azrael3000
@ArnoldNeumaier, je suis désolé, je ne vois pas comment cela est possible. Pourriez-vous expliquer?
Pedro
lisse en fonction de , mais est une fonction non lisse de , et l'intégration est sur , pour autant que j'ai compris la question. La fonction composite est donc une fonction non lisse de . q r r rqqrrr
Arnold Neumaier
1
@Pedro Je l'ai implémenté et cela fonctionne comme un charme. Nous avons également trouvé une solution analytique aujourd'hui. Mais ce n'est que pour un cas particulier qui peut être utilisé pour reconstruire le cas général. Cela signifie que nous devons effectuer une décomposition de domaine. Puisque le Romberg converge en environ 4 étapes, je pense que pour cette raison, ce sera plus rapide que d'utiliser la formule analytique. Et selon Wikipedia, nous pouvons faire encore mieux que Romberg en utilisant des polynômes rationnels. Vous trouverez votre nom dans les remerciements de mon prochain article :) A bientôt.
Azrael3000
2

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.

GertVdE
la source
Je pense que le problème ici est que la fonction est un polynôme de , mais est une fonction non linéaire dans les coordonnées spatiales. La fonction est lisse jusqu'au bord de la fonction de base, mais pas polynomiale, sauf le long des axes. qqq
Pedro
C'est bien Pedro.
Azrael3000
Ah ok. mon erreur. Désolé.
GertVdE
2

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 0r0r0

q q r r r 0 r rW 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 .qqrrr0rr

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 0r0CjenFr0

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 T1 r 0 N T r 0 r 0NT=1r0NT1r0NTr0r0

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

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 xr0=0XyX

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. qq2q

Arnold Neumaier
la source
Je pense qu'il y a un petit malentendu. J'ai mis à jour la description de ma question. En fait, dans l'intégrale ne peut jamais être inférieur à 0. Et n'est pas nécessairement contenu dans le triangle. r 0qr0
Azrael3000
Votre nouvel ajout n'a pas de sens pour moi. Si alors doit être r R 3 r 0r0R3r . Ou vous intégrez-vous sur un triangle 2D dans ? - Je n'ai pas que est dans le triangle. Je viens d'ajouter un peu plus de détails à ma réponse. R3r0
Arnold Neumaier
Oui, il est vrai que j'intègre sur un triangle 2D dans . R3
Azrael3000