Systèmes de Vandermonde

Les systèmes d'équations linéaires de Vandermonde (en l'honneur du mathématicien français Alexandre Vandermonde (1735-1796)) apparaissent naturellement lorsqu'on veut trouver un polynôme de degré N-1 qui prend des valeurs données q_i en N points donnés distincts x_i (polynôme d'interpolation de Lagrange). On est alors amené à résoudre le système d'équations linéaires

$$w_1+w_2x_i+w_3x_i^2+...+w_Nx_i^{N-1}=q_i \quad 1\leq i \leq N$$

par rapport aux inconnues w_j où

$$P(x)=\sum _{j=1}^{N}w_jx^{j-1}$$

est le polynôme cherché. Sous des habillages variés, des systèmes de ce type se présentent dans diverses questions d'algèbre linéaire.

Remarquons que souvent il n'est pas nécessaire de connaître la décomposition du polynôme de Lagrange dans la base naturelle, mais que sa décomposition dans une autre base suffit. C'est ainsi que si on se contente de la décomposition de P(x) dans la base contituée des polynômes $P_j(x)= \prod_{n=1 \atop {n \neq j}}^N {{x-x_n} \over {x_j-x_n}}$ qui valent 1 en x_j et aux autres points x_n il n'y a pas de calcul à faire car dans ce cas les coefficients de la décomposition sont les valeurs q_i. De même si on choisit de décomposer P(x) dans la base de Newton 1,x-x₁,(x-x₁)(x-x₂),...,(x-x₁)(x-x₂)...(x-x_N-1) le calcul se ramène à l'algorithme simple des différences divisées. Cependant il existe des problèmes concrets dont la solution nécessite la résolution effective d'un système de Vandermonde, par exemple quand on veut connaître explicitement les coefficients du polynôme d'interpolation de Lagrange. Dans le sujet qui suit on travaille en fait avec la matrice transposée de la précédente.

1.Algorithme

Soit N un entier $\geq 2$. Etant donné x=(x₁,x₂,...,x_N) un N-uplet de réels deux à deux distincts on note B la matrice

$$B=\left( \begin{array} {cccc} 1 & 1 & \cdots & 1 \\ x_1 &... ..._1^{N-1} & x_2^{N-1} & \cdots & x_N^{N-1} \end{array} \right)$$

On cherche à résoudre le système

BW=Q

où Q est la matrice colonne constituée des seconds membres $q_1,q_2,\cdots,q_N$ du système et où W est la matrice colonne constituée des inconnues $w_1,w_2,\cdots,w_N$ du système.

Pour tout entier j vérifiant $1 \leq j \leq N$ on pose

$$P_j(x)= \prod_{n=1 \atop {n \neq j}}^N {{x-x_n} \over {x_j-x_n}}$$

P_j est donc un polynôme en x de degré N-1 qui peut s'écrire

$$P_j(x)=\sum_{k=1}^N {A_{j,k}x^{k-1}}$$

Montrer que la matrice (A_j,k)_j,k est l'inverse de la matrice B. En conclure que pour tout $1 \leq j \leq N$

$$w_j=\sum_{k=1}^N {A_{j,k}q_k}.$$

Nous allons dans la suite mettre en place une méthode de résolution qui calcule les coefficients des polynômes P_j, donc qui calcule l'inverse de la matrice B. Pour calculer les coefficients de P_j on sera amené à calculer les coefficients de

$$N_j(x)= \prod_{n=1 \atop {n \neq j}}^N {x-x_n}$$

et aussi le dénominateur intervenant dans la formule qui définit P_j, c'est à dire le nombre N_j(x_j).

Posons

P(x)=(x-x₁)(x-x₂)...(x-x_N)

P(x) est donc un polynôme de degré N qui s'écrit sous la forme:

P(x)=x^N+c_Nx^N-1+...+c₂x+c₁

Montrons tout d'abord comment, si on connaît les coefficients c_j, on peut calculer les coefficients du polynôme

$$N_j(x)= \prod_{n=1 \atop {n \neq j}}^N {x-x_n}$$

Pour cela posons

N_j(x)=b_Nx^N-1+...+b₂x+b₁

Etablir que

$$\left\{ \begin{array} {c} b_N=1\\ b_{k-1}=c_k+x_jb_k \end{array} \right.$$

Connaissant les coefficients de N_j il est alors facile de calculer le dénominateur N_j(x_j) intervenant dans la définition de P_j.

En effet posons t_N=b_N=1 et définissons pour tout $k \leq N$

t_k-1=x_jt_k+b_k-1

Montrer que t₁=N_j(x_j).

Il reste maintenant à calculer les coefficients c_j de P.

Pour tout entier k vérifiant $1 \leq k \leq N$ on définit

Q_k(x)=(x-x₁)(x-x₂)...(x-x_k)

et on écrit Q_k sous la forme

$$Q_k(x)=x^k+\alpha _{k,k}x^{k-1}+\alpha _{k,k-1}x^{k-2}+...+\alpha _{k,1}$$

Etablir que

$$\alpha _{1,1}=-x_1$$

et pour k=2,3,...,N

$$\alpha _{k,k}=\alpha _{k-1,k-1}-x_k$$

$$\alpha _{k,j}=\alpha _{k-1,j-1}-x_k\alpha _{k-1,j}~~~j=k-1,...,2$$

2.Programmation

On dispose d'un programme qui résout les systèmes du type que l'on vient d'étudier. Une procédure a été laissée en blanc, celle qui correspond au calcul des coefficients de P(x).

Le programme déclare une constante VraiDim qui est la dimension du système (le N des paragraphes précédents). Pour connaître la matrice B, il suffit de connaître (x₁,x₂,...,x_Vraidim) si bien qu'on définit un type

VDM_Mat=ARRAY[1..VraiDim] OF REAL;

De même le polynôme P devra être stocké en rangeant les coefficients c₁,c₂,...,c_N dans un tableau. On définit donc un type

VDM_Poly=ARRAY[1..VraiDim] OF REAL;

Ecrire la procédure

PROCEDURE PolyNoyau(X:VDM_Mat;VAR Noyau:VDM_Poly);

qui étant donné une matrice de Vandermonde définie par les nombres

(x₁,x₂,...,x_VraiDim)

rangés dans le tableau X (x_j est rangé dans X[j]) fait ressortir dans la variable Noyau les coefficients c₁,c₂,...,c_VraiDim (c_j dans Noyau[j]) du polynôme P(x)=(x-x₁)...(x-x_VraiDim) écrit sous la forme

P(x)=x^VraiDim+c_VraiDimx^VraiDim-1+...+c₁.

$\bullet$ Partie théorique

En effectuant le produit de la matrice A=(A_j,k)_j,k par la matrice B on constate que

AB=(P_j(x_k))_j,k

ce qui prouve que A est l'inverse de B.

On peut alors écrire que W=AQ. On obtient ainsi les formules

$$w_j=\sum_{k=1}^N {A_{j,k}q_k}.$$

On vérifie immédiatement sur l'expression de N_j(x) que le coefficient du terme de plus haut degré est 1. En remarquant que P(x)=N_j(x)(x-x_j) on établit la formule

b_k-1=c_k+x_jb_k.

Le calcul proposé alors pour N_j(x_j) n'est autre que l'algorithme de Horner.

Il est facile de voir sur l'expression de Q₁(x)=x-x₁ que $\alpha_{1,1}=1$.

Les autres formules découlent de

Q_k(x)=Q_k-1(x)(x-x_k).

$\bullet$ Programmation

PROCEDURE PolyNoyau(X:VDM_Mat;VAR Noyau:VDM_Poly);

VAR
  I,J    : INTEGER;
  Auxil  : REAL;

BEGIN
  Noyau[1]:=-X[1];
  FOR I:=2 TO VraiDim
    DO BEGIN
         Auxil:=-X[I];
         Noyau[I]:=Noyau[I-1]+Auxil;
         FOR J:=I-1 DOWNTO 2
           DO Noyau[J]:=Auxil*Noyau[J]+Noyau[J-1];
         Noyau[1]:=Auxil*Noyau[1];
       END;
END;

On peut voir que le nombre d'opérations à faire dans cet algorithme est en O(N²), la partie la plus coûteuse étant le calcul des coefficients c_k.