La programmation de la synthèse polynomiale

Nous allons maintenant passer à la généralisation de l'algorithme précédent pour un polynôme de degré quelconque. Toute la difficulté vient de la construction de la matrice principale et du second terme pour un degré quelconque.

Recherchons d'abord les termes du polynôme. Il suffit de trouver tous les termes de degré inférieurs ou égaux au degré du polynôme cherché. Pour le degré 0, il n'y a qu'un seul terme de degré 0 donc fonction de 1 :
 

1

Pour le degré 1 nous avons le terme du degré 0 et deux termes de degré 1, un fonction de x et l'autre de y
 

1 x y

Pour le degré 2 nous avons le terme du degré 0, les deux termes de degré 1 et 3 termes fonction de x et y, un fonction de x2, un fonction de xy et un fonction de y2
 

1 x y x2 xy y2

En fait pour trouver les termes de degré n, il suffit de prendre pour premier terme xny0=xn, puis d'enlever 1 à l'exposant de x tout en ajoutant 1 a l'exposant de y. Le terme suivant sera x(n-1)y, et ceci jusqu'au terme en yn. Pour trouver les termes de la matrice principale, il suffit alors de créer un tableau dont la première ligne et la première colonne contiendront la sommes des termes du polynôme. Par exemple pour le degré 2 :
 
         
         
         
         
         

Puis de compléter chaque case vide en multipliant le terme en haut de la colonne par le terme a gauche de la ligne :
 
n
 

Evidemment, il ne s'agit pas d'une vraie multiplication mais d'une vue de l'esprit qui nous permettra de trouver l'algorithme. Pour obtenir le second terme, il suffit multiplier toutes les sommes des termes par z :
 

 
Voyons maintenant la programmation. D'abord écrivons une petite routine qui calcule la puissance b d'un nombre a :
function Puiss(a:Longint;b:Integer):Extended;
var
i:Integer;
Valr:Extended;
begin
if b>1 then
begin
Valr:=1.0*a;
for i:=1 to b-1 do
Valr:=1.0*Valr*a;
end

else if b=1 then Valr:=1.0*a
else if b=0 then Valr:=1.0;
Result:=Valr;
end;

 

Dans la procédure principale, on commence par trouver la dimension de la matrice principale :
Dim:=1;
for i:=1 to Degre do Dim:=Dim+i+1;
et la dimension du tableau de sommes :
Dims:=2*Degre+1;
En effet, on peut remarquer que le degré maximum des sommes est double de celui du polynôme. Puis calculons les sommes :
for i:=1 to dims do
for j:=1 to dims do
sums^[(i-1)*dims+j]:=0;
for j:=1 to dims do
begin
m:=1;
k:=j-1; l:=0;
while l<j do
begin
for o:=1 to n do
Sums^[(m-1)*Dims+j]:=Sums^[(m-1)*Dims+j]+Puiss(xi^[o],k)*Puiss(yi^[o],l);
Inc(m);
k:=k-1; l:=l+1;
end;
end;
Pour le degré deux, ce tableau sera le suivant :
 
n 0 0 0 0
0 0 0
0 0
0
 

Pour pouvoir remplir la matrice principale, il faut savoir à quelle position dans la matrice va correspondre quelle somme. Pour ce faire, imitons le processus que nous avons utilisé précédemment pour remplir la matrice principale et cherchons les puissances de la ligne supérieure de celle ci. Le tableau Rang1[i] contiendra la puissance de x de la colonne i et le tableau Rang2[i] contiendra la puissance de y de la colonne i.

Rang1^[1]:=0; Rang2^[1]:=0;
m:=1;
for i:=1 to Degre do
begin
k:=i; l:=0;
while l<i+1 do
begin
Rang1^[m+1]:=k; Rang2^[m+1]:=l;
Inc(m);
k:=k-1; l:=l+1;
end;
end;
Pour le degré 2 cela donnera :
Pour Rang1 :
 
0 1 0 2 1 0

et pour Rang2 :
 
0 0 1 0 1 2

Il est alors facile de trouver les sommes correspondant à chaque coefficient de la matrice principale. Il suffit de trouver les puissances de x et y correspondant à ces coefficients grâce aux tableaux des rangs, puis de faire la correspondance avec le tableau des sommes :

for i:=1 to Dim do
for j:=1 to Dim do
begin
k:=Rang1^[j]+Rang1^[i];
l:=Rang2^[j]+Rang2^[i];
Mat^[(i-1)*Dim+j]:=Sums^[(l+1-1)*Dims+k+l+1];
end;
Il ne reste plus alors qu'à compléter le second terme de la même manière que les sommes pour pouvoir résoudre le système par la méthode de Gauss.
for i:=1 to Dim do yg^[i]:=0;
for o:=1 to n do yg^[1]:=yg^[1]+zi^[o];
m:=1;
for i:=1 to degre do
begin
k:=i; l:=0;
while l<i+1 do
begin
for o:=1 to n do
yg^[m+1]:=yg^[m+1]+puiss(xi^[o],k)*puiss(yi^[o],l)*zi^[o];
inc(m);
k:=k-1; l:=l+1;
end;
end;
L'algorithme de résolution est ici légèrement différent de celui du degré deux car il permute certaines colonnes pour éviter les problèmes dus aux pivots nuls et pour diminuer les effets de troncature.

Voici la procédure que vous pouvez télécharger ici : Synthèse degré n (1K).

Remarque :
Ne connaissant pas à priori le degré du polynôme demandé par l'utilisateur, les divers tableaux devront être dynamiques. Or en Pascal, les tableaux dynamiques à deux dimensions n'existent pas, il faut donc utiliser des tableaux à une dimension dont on calcule la position. Ainsi, si le tableau T a pour dimensions Dim sur Dim, l'élément T[i,j] devient T^[(i-1)*Dim+j].

Page précédente