Comment calculer numériquement les résidus?

15

J'ai besoin de calculer l'intégrale suivante: Où est une matrice (cinétique d'une seule particule et énergie potentielle exprimée en base), est une matrice qui dépend de (une seule particule plusieurs- la fonction du corps vert) et l'intégrale du contour est un demi-cercle gauche. L'intégrande a des pôles sur l'axe réel négatif et son évaluation est coûteuse. Quelle est la façon la plus efficace de calculer une telle intégrale?

12πiCf(E)dE
f(E)=Tr((h+E)G(E))
hGEf(E)

Voici mes recherches jusqu'à présent:

1) J'utilise l'intégration gaussienne, mon chemin d'intégration est un rectangle. J'ai fixé le côté gauche et droit (c'est-à-dire la largeur) et joué avec la hauteur (au-dessus et en dessous de l'axe réel) de telle sorte que pour l'ordre d'intégration donné, j'obtiens la plus grande précision. Par exemple pour l'ordre 20, si la hauteur est trop grande, la précision diminue (évidemment), mais si elle est trop petite, elle diminue également (ma théorie est qu'il faut de plus en plus de points autour des pôles à mesure que la hauteur augmente 0). Je me suis installé avec une hauteur optimale de 0,5 pour ma fonction.

2) Ensuite, j'ai défini le côté droit du rectangle à E0, généralement E0 = 0, mais cela pourrait être E0 = -0,2 ou quelque chose de similaire.

3) Je commence à déplacer le côté gauche du rectangle vers la gauche et pour chaque étape, je fais la convergence de l'ordre d'intégration pour m'assurer que mon intégrale est complètement convergée pour chaque rectangle. En augmentant la largeur, j'obtiens finalement une valeur convergente dans la limite du demi-cercle gauche infini.

Le calcul est vraiment lent et pas très précis pour les grandes largeurs. Une amélioration consiste à simplement partitionner la grande largeur en "éléments" et à utiliser l'intégration gaussienne sur chaque élément (comme dans FE).

Une autre option serait d'intégrer un petit cercle autour de chaque pôle et de le résumer. Problèmes:

a) Comment trouver numériquement les pôles de la fonction ? Il doit être robuste. La seule chose que je sais, c'est qu'ils sont sur l'axe réel négatif. Pour certains d'entre eux (mais pas tous), je connais également une assez bonne estimation initiale. Existe-t-il une méthode qui fonctionne pour toute fonction analytique ? Ou cela dépend-il de la forme réelle de ?f ( E ) f ( E )f(E)f(E)f(E)

b) Une fois que nous connaissons les pôles, quel schéma numérique est le meilleur pour intégrer le petit cercle autour de lui? Dois-je utiliser l'intégration gaussienne sur un cercle? Ou dois-je utiliser une distribution uniforme des points?

Une autre option pourrait être qu'une fois que je connais les pôles grâce à a), il pourrait y avoir un moyen semi-analytique d'obtenir les résidus sans avoir besoin de l'intégration complexe. Mais pour l'instant, je serais heureux d'optimiser simplement l'intégration des contours.

Ondřej Čertík
la source
1
Avez-vous vérifié le livre "Méthodes numériques pour l'inversion de transformation de Laplace" de Cohen (2007)? L'IIRC, Robert Piessens (de la renommée QUADPACK) a également travaillé sur ce sujet.
GertVdE

Réponses:

7

Je peux offrir une suggestion pour votre première question: si vous savez que vos pôles sont quelque part le long de l'axe réel, vous pouvez les localiser assez efficacement en utilisant l' interpolation / approximation rationnelle . Cela revient à trouver des polynômes et q ( x ) tels quep(x)q(X)

F(X)p(X)q(X)

XF(X) q(X)

L'interpolation / approximation rationnelle peut être une chose délicate, mais j'ai récemment co-écrit un article sur un algorithme stable pour les calculer à l'aide du SVD. L'article contient du code Matlab implémentant l'algorithme, et une version plus étendue de celui-ci est disponible en tant que fonction ratinterpdans le projet Chebfun , dont je suis l'un des développeurs.

Pour votre deuxième question, ce document peut être utile.

Pedro
la source
Merci pour tous les conseils! Voici le code netlib.org/toms/579 du papier Bengt Fornberg. Malheureusement, il y a un bug numérique, car c'est la sortie que j'obtiens : gist.github.com/2942970#file_output . Je vais donc devoir le réimplémenter ou le déboguer. Le lien Chebfun me donne 404 (je l'ai essayé il y a quelques mois avec les mêmes résultats, alors peut-être que cela ne fonctionne tout simplement pas aux États-Unis).
Ondřej Čertík
@ OndřejČertík: Je n'ai jamais utilisé le code TOMS 579 moi-même, donc je ne sais pas quoi vous dire sur les erreurs. Quant à la page d'accueil de Chebfun, pourriez-vous essayer de la "googler" et de voir si cela fonctionne?
Pedro
Google trouve la page d'accueil de Chebfun et affiche les versions mises en cache. Mais quand je clique sur la page, voici
Ondřej Čertík
Essayez un autre navigateur? Ou d'un autre FAI. Le site Web fonctionne bien à partir d'ici (aux États-Unis.)
Costis
J'ai essayé Firefox et Chrome. Il doit donc par mon FAI. Bizarre.
Ondřej Čertík