Approximation de Padé

Au voisinage d'un point le développement de Taylor d'une fonction donnée fournit un polynôme qui réalise une approximation de la fonction. Mais on peut également chercher une fraction rationnelle qui satisfait le même type de conditions que la partie polynômiale de Taylor : l'égalité des dérivées jusqu'à un certain ordre. On peut alors comparer les 2 types de développements suivant le comportement de la fonction. Cette technique est à la base des approximations dites de Padé [15], donnant pour une fonction des fractions rationnelles l'approchant avec une précision croissante quand on augmente le degré du numérateur et du dénominateur.

1.Algorithme

Soit f une fraction rationnelle, $f(x)={{\textstyle P(x)}\over{\textstyle Q(x)}}$ où P et Q sont des polynômes de degré n :

$$P(x)=\sum_{i=0}^{n}a_ix^i\quad Q(x)=\sum_{i=0}^{n}b_ix^i$$

On écrit Qf=P : dériver cette relation à l'ordre k pour $k=0,\dots,2n$.

Ecrire les relations obtenues pour x=0.

En déduire que si les dérivées de f en x=0 sont connues on peut calculer les coefficients a_i et b_i en résolvant un système de dimension 2n+1 ; calculer ses coefficients et son second membre. Quel choix arbitraire supplémentaire reste-t-il à faire ?

Pour une fonction f quelconque dont les dérivées en x=0 sont connues on réalisera une approximation en la remplaçant par la fraction rationnelle obtenue grâce aux coefficients a_i et b_i.

2.Programmation

On dispose d'un programme qui calcule les coefficients a_i et b_i et compare l'approximation d'une fonction connue avec ses valeurs exactes d'une part et avec les valeurs obtenues grâce à sa série de Taylor.

Le programme utilise pour cela les types

Polynome = ARRAY[0..PolyMax] OF REAL;

destiné à contenir les coefficients a_i d'un polynôme $P(x)=\sum_{i=0}^{n}a_ix^i$ et

Subdivision = ARRAY[0..SubdivMax] OF REAL;

pour les valeurs des dérivées de f en 0.

Ecrire la procédure et la fonction suivantes dont seul l'en-tête a été indiqué :

PROCEDURE TableauCoeff(D:Subdivision;N:INTEGER;VAR P,Q:Polynome);

qui calcule à partir des dérivées de f données dans le tableau D les coefficients des polynômes P et Q ; on pourra écrire une fonction annexe pour calculer les coefficients de la matrice du système.

ATTENTION : Déclarer les entiers utilisés de type LongInt afin d'éviter des problèmes de débordement dans les calculs des coefficients

FUNCTION ApproxTaylor(D:Subdivision;N:INTEGER;T:Real):Real;

qui calcule la valeur approchée $$\sum_{i=0}^{n}D[i]{{\textstyle T^i}\over{\textstyle i!}}$$

Pour résoudre le système linéaire on utilisera la procédure déjà incluse dans le programme

Gauss(N:INTEGER;A:Matrice;B:Vecteur;VAR X:Vecteur;VAR Err:BOOLEAN);

qui résout le système AX=B de dimension $N\times N$ et met le booléen Err à vrai si le système n'a pas de solution.

$\bullet$ Partie théorique

La formule de Leibnitz donne à l'ordre k :

d'où l'égalité :

$$\sum_{j=0}^kC_k^j \sum_{i=j}^{n}i(i-1)\dots(i-j+1)b_ix^{i-j}f^{(k-j)}(x)$$

$$=\cases{ \sum_{i=k}^{n}i(i-1)\dots(i-k+1)a_ix^{i-k}\qquad k\le n\cr \ \cr 0\qquad k\gt n\cr}$$

Pour x=0 on obtient :

Ces relations forment un système linéaire de 2n+1 équations aux 2n+2 inconnues $a_0,\dots,a_n,b_0\dots,b_n$. Il reste un coefficient à choisir car la fraction ${{\textstyle P}\over{\textstyle Q}}$ n'est définie qu'a un coefficient multiplicatif près au numérateur et au dénominateur. On peut prendre par exemple a₀=1, b₀=1, a_n=1 ou b_n=1. Le choix b₀=1 permet d'être sûr que le dénominateur Q ne risque pas de s'annuler. Le système se met alors sous la forme :

a₀ = f(0)

$\bullet$ Programmation

FUNCTION D(K,P:INTEGER):LongInt;
VAR J:LongInt;
    Q:INTEGER;
BEGIN
  J:=1;
  IF K>0 THEN FOR Q:=K DOWNTO K-P+1 DO J:=J*Q;
  D:=J;
END;

PROCEDURE TableauCoeff(Der:Subdivision;N:INTEGER;VAR A,B:Polynome);
VAR M:Matrice;
    X,Y:Vecteur;
    I,J:INTEGER;
BEGIN
{Affectation coefficients matrice dimension 2n+1}
{A: numerateur   B: denominateur}
  FOR I:=1 TO 2*N+1 DO Y[I]:=0;
  Y[N+1]:=D(N,N);
  FOR I:=1 TO N+1 DO
  BEGIN
    FOR J:=1 TO I DO M[I,J]:=Der[I-J]*D(I-1,J-1);
    FOR J:=I+1 TO 2*N+1 DO M[I,J]:=0;
    M[I,N+1+I]:=-D(I-1,I-1);
  END;
  FOR I:=N+2 TO 2*N+1 DO
  BEGIN
    FOR J:=1 TO N+1 DO M[I,J]:=Der[I-J]*D(I-1,J-1);
    FOR J:=N+2 TO 2*N+1 DO M[I,J]:=0;
  END;
  {Solution du systeme}
  Gauss(2*N+1,M,Y,X,Err);
  IF ERR THEN
  BEGIN
   writeln('erreur !');
   pause;
  END;
  Polynome_Nul(B);
  FOR I:=1 TO N+1 DO {Coefficients B}
  BEGIN
    B[I-1]:=X[I];
  END;
  Polynome_Nul(A);
  A[N]:=1; {Normalisation}
  FOR I:=N+2 TO 2*N+1 DO {Coefficients A}
  BEGIN
    A[I-N-2]:=X[I];
  END;
END;

{Approximation de Taylor en 0}
FUNCTION ApproxTaylor(d:Subdivision;N:INTEGER;T:Real):Real;
VAR X,Z,P,Q:Real;
    I,J:INTEGER;
BEGIN
  Z:=D[0];
  P:=1;
  FOR I:=1 TO N DO
  BEGIN
    P:=P*T/I;
    Z:=Z+P*D[I];
  END;
  ApproxTaylor:=Z;
END;

$\bullet$ Remarques d'exécution

Pour la fonction exponentielle on constate que pour de petites valeurs de n (2, 3, 4) l'approximation de Padé est de meilleure qualité que celle de Taylor. Lorsque n est plus grand (jusqu'à 9, valeur maximale compatible avec une résolution sans problème du système) les 2 approximations semblent comparables. Dans tous les cas l'approximation n'est acceptable que dans un intervalle assez restreint autour de l'origine.

On remarque aussi que les coefficients du numérateur et du dénominateur sont les mêmes au signe près, à cause de $e^{-x}={{\textstyle 1}\over{\textstyle e^x}}$.

Les approximations par fractions rationnelles qui ont été introduites par Padé au début du siècle n'ont cessé de se développer et elles sont utilisées actuellement dans de nombreux domaines où on doit calculer des valeurs approchées de fonctions (calculatrices, bibliothèques numériques), y compris même pour des fonctions matricielles [6].