Discussion :

[Ankhotep]

Bonjour, et merci de l'attention que vous porterez à ma question

Je cherche un algorithme pour calculer les (coefficients de) polynômes orthogonaux pour un certain poids f. Ma première idée était d'utiliser l'algorithme de Gram-Schmidt, mais celui-ci s'avère trop lent.

En fait, le but est de calculer des déterminants de Hankel d'une fonction f. La première manière de le faire est simplement de créer la matrice de Hankel (qui nécessite donc 2n-1 intégrations (n étant l'ordre du déterminant de Hankel)) et ensuite d'en prendre le déterminant.

Mais grâce à un théorème, on sait que ce même déterminant est en fait un produit des n-1 coefficients dominants des n-1 premiers polynômes orthogonaux de f. Je me demande donc s'il est possible de calculer le déterminant de Hankel de cette manière. Seulement pour l'instant (avec Gram-Schmidt), cela requiert de l'ordre de n^2 intégrations, ce qui est beaucoup plus qu'avec la manière "simple" et directe de le faire.

La raison pour laquelle je cherche à utiliser les polynômes orthogonaux est la suivante : quand n devient grand (en fait à partir de n ~= 20), le déterminant de Hankel devient énorme, ce qui engendre des erreurs de calculs bien trop importantes. Hors mon but est d'ensuite diviser par une constante qui est aussi énorme. Je voulais donc "diluer" cette division à travers le produit des coefficients orthogonaux...

Malheureusement Gram-Schmidt est vraiment trop long (~= 10 secondes pour n = 6 contre ~= 2 pour la méthode "simple").

La question est donc : y a-t-il un algorithme efficace (de l'ordre de n en fait) pour calculer les n premiers polynômes orthogonaux de f (f pouvant être absolument immonde :p (lisse, tout de même)).
Merci d'avance

[Prof]

Bonjour.

La réponse à votre question est : OUI.

Les polynômes orthogonaux sont connus à un facteur multiplicatif près.
Si on les prend orthonormaux, il y a unicité ; mais la normalisation introduit des racines carrées qui compliquent les calculs, alors qu'elles ne servent à rien.
Il est plus simple de prendre P(n) = X^n + a(n)X^(n-1) + b(n)X^(n-2) + ... sans chercher à normaliser P(n).
Dans ce cas, la matrice de passage entre la base 1,X,X^2,...,X^n et la base P(0),P(1),P(2),...,P(n) est triangulaire supérieure avec des 1 sur la diagonale.
Il en résulte que le déterminant de la matrice du produit scalaire dans la base 1,X,...,X^n est égal au déterminant de la matrice du produit scalaire dans la base P(0),P(1),...,P(n).
Or la matrice du produit scalaire dans la base P(0),P(1),P(2),...,P(n) est diagonale,
avec a(i,i) = <P(i),P(i)>.
Le déterminant de Hankel cherché est donc le produit des <P(i),P(i)>.
Tout revient à calculer simplement les polynômes P(i) pour i= 0,1,...,n.

C'est ici qu'intervient une remarque importante.
L'article de Wikipédia consacré aux polynômes orthogonaux signale une propriété intéressante :
il y a une relation de récurrence simple entre P(n), P(n-1) et P(n-2).
( voir http://fr.wikipedia.org/wiki/Polyn%C3%B4mes_orthogonaux )

Avec les polynômes P(n) sous la forme P(n) = X^n + ..., la relation de récurrence se met sous la forme
P(n) = (X + alpha(n))P(n-1) + beta(n)P(n-2).
Il ne reste plus qu'à calculer les deux coefficients alpha(n) et beta(n).
Une identification des coefficients de X^(n-1) et X^(n-2) donne leur valeur en fonction de a(n), a(n-1), b(n), b(n-1).
Puis la valeur de P(n) donnée par l'algorithme de Gram-Schmidt permet d'exprimer alpha(n) et beta(n) à l'aide de produits scalaires.

Voici l'algorithme correspondant en MAPLE :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 
> ps:=(f,g)->int(f*g*w,x=a..b);
 
> Poly:= proc(n) 
            local P,a,b,pp,produit,k,alpha,beta;
            a[0]:=  0 ;
            b[0]:=  0 ;
            P[0]:=  1 ;
            pp[0]:= ps(P[0],P[0]) ;
            a[1]:=  - ps(x,P[0])/pp[0] ;
            b[1]:=  0 ;
            P[1]:=  x+a[1] ;
            pp[1]:= ps(P[1],P[1]) ;
            produit:= pp[0]*pp[1] ;
            for k from 2 to n do
              a[k]:=  -ps(x^k,P[k-1])/pp[k-1] ;
              b[k]:=  a[k]*a[k-1]- ps(x^k,P[k-2])/pp[k-2] ;
              alpha:= a[k]-a[k-1] ;
              beta:=  b[k]-b[k-1] - alpha*a[k-1] ;
              P[k]:=  (x + alpha)*P[k-1]+beta*P[k-2] ;
#              print(sort(expand(P[k]))) ; # pour vérification des polynômes
              pp[k]:= ps(P[k],P[k]) ;
              produit:= produit*pp[k] ;
            od ;
            RETURN(produit) ;
         end :

Attention : il y a les variables globales a et b, qui interviennent dans le produit scalaire, et les variables locales a et b dans la fonction Poly.
Afin de les distinguer, la déclaration des variables locales est donc indispensable.

Les variables P,a,b,pp sont des tableaux car elles sont utilisées à plusieurs endroit.
pp[i] vaut <P(i),P(i)> et est l'un des facteurs du produit cherché.
Les variables alpha et beta étant utilisées de façon immédiate n'ont pas besoin d'être stockées dans des tableaux.
D'ailleurs on pourrait les supprimer, mais le code en serait moins compréhensible.
Sous MAPLE, une fonction retourne par défaut la valeur de la dernière instruction ; RETURN n'est donc pas indispensable mais je l'ai placé pour faciliter la traduction de cet algorithme dans un autre langage.
( je ne connais pas MATLAB )

J'ai vérifié que cet algorithme donnait les bons polynômes, en contrôlant avec les polynômes orthogonaux classiques ( Legendre, Tchebicheff, Laguerre, Hermite ).
Puis j'ai neutralisé la vérification en utilisant le caractère # qui signale un commentaire.

Cet algorithme donne les polynômes orthogonaux avec un nombre d'intégrations en O(n), alors que la méthode de Gram-Schmidt les donne avec O(n^2) intégrations..
Plus précisément, il donne P(0),P(1),...,P(n) en utilisant 1 + 2 + 3(n-2) = 3n-3 intégrations.
C'est ce que vous souhaitiez.

Exemples :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
> w:=1: a:=-1: b:=1: Poly(6);
 
   2147483648 / 7888446990683634375
 
> w:=1: a:=-1: b:=1: Poly(10);
 
  154742504910672534362390528 / 1927983533652930855481078826533672813447199742802734375

Pour n=20, le résultat ne tient pas dans une ligne.
Pour avoir des valeurs réelles ( et des calculs plus rapides ), il suffit de remplacer l'instruction " P[0]:= 1 ; " par " P[0]:= 1.0 ; ".
On obtient alors :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
 
> w:=1: a:=-1: b:=1: Poly(6);
 
                           .2722314861 10^(-9)
 
> w:=1: a:=-1: b:=1: Poly(10);
 
                          .8026132078 10^(-28)
 
> w:=1: a:=-1: b:=1: Poly(20);
 
                          .2544999345 10^(-96)

C'est mieux ...

[Ankhotep]

Bonjour Prof,

Merci beaucoup pour votre réponse détaillée ! Utiliser la relation des trois termes semble en effet une bonne idée. J'aurais dû y penser moi-même...