Ce WE, j'ai codé un peu et cherché des sources
Je passe par Gauss-Tchebychev, je code sur [-1;1], mais ça ne passe pas ma phase de test. Surestimation et vitesse de convergence douteuse, je revérifie tout, met des 2.0 à la place des 2 pour éviter les problèmes de division entière etc... galère totale. Je test des méthode modifiées, mais rien à faire, il y a un truc qui ne marche pas.
Sur wiki, y a Gauss-Kronrod (coef+ poids) pour n=15, je code ça en 5 minutes sur [-1;1], test polynômes : trop bon et fn ok, Je modifie sur [a;b] et hop du 10^-16 de précision du premier coup.
Je teste un Gauss-Legendre avec seulement 3 points + découpage en 100 de l’intervalle : trop bon comme résultats.
je regarde un peu sur internet les sources de quadrature de Gauss.
Là, je tombe sur le sources de QUADPACK pour Fortran77 le tout écrit entre 1981 et 1983.
Le fortran, c'est un peu lourd à suivre car toute variable est globale, donc, il faut systématiquement revenir à la dernière occurrence qui l'a modifiée et il faut suivre les sauts.
Finalement, il y avait à l'époque l'ago complet qui inclus tout ce dont on a besoin.
je résume grossièrement le fonctionnement :
soit f définie sur [a,b], à intégrer sur une portion [l,r]
si l n'est pas a et r n'est pas b sont loin, je suppose f régulière, je fais un Gauss-X à n et 2n points où X est un méthode où les points "n" et "2n" coïncident (Kronrod par ex), l'erreur retournée est la différence entre les deux calculs et la valeur, le calcul en 2n.
si on est en a, Gauss-Tchebychev à 12 et 24 points, mais avec la modif de f en f(x)*(x-a)^alpha, il y a une routine qui donne les nouveaux coefs à appliquer, on retourne l'erreur
si on est en b, Gauss-Tchebychev à 12 et 24 points, mais avec la modif de f en f(x)*(b-x)^beta, il y a une routine qui donne les nouveaux coefs à appliquer on retourne l'erreur
si a=oo ou b=oo, il y a une autre routine qui marche avec un développement en exp(-x)
de plus, il y a une fonction norme infinie qui permet de déterminer si f est régulière en a et b.
au final, on somme les morceaux, avec itération possible si on travaille si on veut sur l'erreur cumulée.
Un sous programme donne les poids et abscisses de Gauss-X
Comme l'intégration choisie a été [l,r] et pas [-1;1], on "passe" comme argument
a,b,l,r,alpha,beta, les vecteurs position et poids, le type d'intégration.
en retour : "a,b,l,r,alpha,beta, le type d'intégration" ne sont pas modifiés
: le vecteur position et les vecteurs poids sont renvoyé pour tous chaque types d'intégration.
Les calculs sont touts des récurrences pas très longues.
les source fortran :
http://www.netlib.org/quadpack/dqc25c.f est la première où j'ai vu cet algo, elle demande certaines autres sources comme dqmomo.f qui recalcule tous les poids.
Voilà, Aviateur, il suffisait d'être universitaire dans les années 80.