Méthodes de dérivation numérique

Rappels d'analyse

Je pars du principe que vous connaissez la définition lycéenne de la dérivée, celle que l'on apprend en Première:
Soit f une fonction continue dans un intervalle donné. Soit un point a appartenant à I. On dit que la fonction f est dérivable en a si et seulement si la limite lim (1/h)*( f(a+h) -f(a)), quand h tend vers 0, existe; cette limite est notée f'(a) dans vos livres de lycée.
On la note plutôt f⁽¹⁾(a) pour tenir compte de l'existence possible de dérivées d'ordre supérieur!
A cette définition de la dérivée, il me faut ajouter un élément supplémentaire non connu des lycéens, mais fondamental en calcul numérique: le développement limité d'une fonction ou encore développement en série de Taylor. Je ne démontrerai rien ici, je me contenterai d'énoncer l'existence!
Donc, pour une fonction f analytique, dans un domaine de convergence de rayon R autour d'un point a, pour tout x tel que a-R <= x <= a+R, je peux écrire:
f(x) = f(a) + (x - a)*f⁽¹⁾(a) + (x - a)²/2! * f⁽²⁾(a) + ... + (x-a)n/n! * f⁽ⁿ⁾(a) + ...
Je vous laisse vous rafraichir la mémoire dans votre bouquin d'analyse préféré...
Comme vous le savez sans doute, le calcul de la dérivée d'une fonction en un point est un exercice fréquent en physique et en chimie: calcul de vitesse, d'extrema de courbe, etc..
Dans la suite on va donc chercher le moyen de calculer la dérivée première f⁽¹⁾ en un point.

Une méthode simple par différences finies

Tout d'abord, remarquons que la définition mathématique d'une dérivée est absolument inopérante à la fois pour les physiciens et pour les informaticiens. Pour la simple raison que faire tendre h vers 0 est impossible, à la fois physiquement et encore plus informatiquement.
Vous savez sans doute que si vous essayez de produire un nombre trop petit dans un programme, votre ordinateur protestera sévèrement et vous jettera! En effet, un nombre est codé avec un certain nombre de bit (16, 32, 64 ou 128), un nombre toujours limité. Aussi, vous ne pouvez pas choisir un nombre aussi petit que nécessaire!
Il va donc falloir trouver autre chose que la pure définition mathématique de la dérivée. Nous allons utiliser un outil très commun en analyse numérique: les différences finies. C'est pour ça que j'ai introduit le développement d'une fonction en série de Taylor. Pas de panique, si on admet l'existence de ce développement, alors le calcul est de niveau Première ou Terminale....
Je cherche donc la valeur de la dérivée d'une fonction f en un point a. Je vais supposer que f est continue et dérivable, et aussi qu'en a, sa dérivée à gauche est égale à sa dérivée à droite. Cette dernière condition n'est pas obligatoire, mais je n'ai pas envie de vous parler de différence à gauche et de différence à droite... Vous creuserez ce cas à titre d'exercice!
Donc, je vais développer ma fonction f(x) sur un intervalle [-h,h] centré en a, en appliquant la formule de Taylor:
f(a + h) = f(a) + h*f⁽¹⁾(a) + h²/2! * f⁽²⁾(a) + h³/3! * f⁽³⁾(a) + h⁴/4! * f⁽⁴⁾(a) + O(h⁵) (1)
f(a - h) = f(a) - h*f⁽¹⁾(a) + h²/2! * f⁽²⁾(a) - h^>3/3! * f⁽³⁾(a) + h⁴/4! * f⁽⁴⁾(a) + O(h⁵) (2)
O(h⁵) désigne tous les termes d'ordre 5 et supérieurs, que je vais négliger sans vergogne, car je fais une erreur de l'ordre h⁵/5! * f⁽⁵⁾(a), en choisissant bien sur h très petit devant 1, par exemple 0.001 ou 0.0001.
Si je retranche (2) à (1), j'obtiens:
f(a + h) - f(a - h) = 2*h*f⁽¹⁾(a) + h³/3! * f⁽³⁾(a)
et donc:
f⁽¹⁾(a) # (1/2h)*(f(a + h) - f(a - h))
si je néglige les termes d'ordre 3 et les suivants. Le signe # signifie "approximativement égal".
Notons que si j'additionne (1) et (2), j'obtiens de la même manière une approximation de la dérivée seconde de f en a:
f⁽²⁾(a) # (1/h²)*(f(a + h) - 2*f(a) + f(a - h))
Vous comprenez sans doute maintenant pourquoi on appelle cette méthode "différences finies"... Il s'agit en l'occurrence ici de différences centrées.
Voyons ce que cela donne sur un exemple dont on connait déjà le résultat. Par exemple, je veux calculer la dérivée de la fonction sin(x) en a=0. Et je vais essayer avec plusieurs valeurs de h, par exemple 0.01 et 0.001
Nous savons bien sur que la dérivée de sin(x) est cos(x), ce qui nous permet de calculer la valeur théorique de la dérivée de sin(x) en a=0, soit cos(0) = 1.
Le programme de calcul se décompose en plusieurs routines:

la routine de définition de la fonction:

double Fonction(double x)

{

return sin(x);

}

la routine de calcul de la dérivée:

double DeriveOrdre1P3(double a, double h, double(*fonction)(double x))

{

double x;

x = (fonction(a+h) - fonction(a-h))/(2*h);

return x;

}

le programme principal:

double x,a,h;

// saisie des données

cout << "Valeur de a: "; cin >> a;

cout << "Valeur de h: "; cin >> h;

x = DeriveOrdre1P3(a,h,Fonction);

cout << "Derivee de f en " << a << ": " << x << endl;

Voyons les résultats obtenus:

pour a = 0 et h = 0.01, j'obtiens 0.999983
pour a = 0 et h = 0.001, j'obtiens 1
pour a = 3.141592 et h = 0.001, j'obtiens -1

Plutôt corrects, ces résultats. Il se peut cependant que pour certaines fonctions un peu "raides", les résultats ne soient pas très précis en utilisant une différence centrée sur 3 points. On peut améliorer un peu les choses en calculant la dérivée sur un voisinage plus étendu, en utilisant une différence à 5 points. Dans l'algorithme ci-dessus, j'ai négligé les termes en O(3). En travaillant un peu le calcul à partir des séries de Taylor (1) et (2), on peut obtenir une équation en O(5), ce qui améliore la précision de deux ordres de grandeur.
Au final, j'obtiens:
f⁽¹⁾(a) # (1/12h)*(f(a - 2*h) - 8*f(a - h) + 8*f(a + h) - f(a + 2*h))
ce qui nous donne la routine:

double DeriveOrdre1P5(double a, double h, double(*fonction)(double x))

{

double x;

x = (fonction(a-2*h) - 8*fonction(a-h) + 8*fonction(a+h) - fonction(a+2*h))/(12*h);

return x;

}

Vous pourrez contrôler qu'on obtient les mêmes résultats pour le calcul précédent.
Notons cependant que ces méthodes souffrent de quelques défauts, en particulier des erreurs d'arrondis et de troncatures. Elles peuvent être très gênantes sur des calculs précis.

Une méthode plus précise de dérivation

Il existe des méthodes qui améliorent grandement la technique des différences finies, en affinant le calcul des ordres supérieurs. Elles utilisent des algorithmes particuliers comme ceux de Neville, qu'il n'est pas question de décrire ici. Comme d'habitude, les "Numerical Recipes" nous fournissent une implémentation d'un algorithme de Neville revu et adapté par Ridders, qui est très efficace.
Vous en trouverez le code dans le programme Derivee.cpp.

Un programme pour tester les différentes méthodes de dérivation

Le programme Derivee.cpp, téléchargeable, permet, sur une fonction déclarée dans la routine Fonction, de comparer la précision des différentes méthodes.