TangenteX.com
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(1)(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(1)(a) + (x - a)2/2!
* f(2)(a) + ... + (x-a)n/n! * f(n)(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(1) en un point.
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 ! 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(1)(a) + h²/2! * f(2)(a)
+ h3/3! * f(3)(a) + h4/4! * f(4)(a)
+ O(h5) (1)
f(a - h) = f(a) - h*f(1)(a) + h²/2! * f(2)(a)
- h>3/3! * f(3)(a) + h4/4! * f(4)(a)
+ O(h5) (2)
O(h5) 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 h5/5! * f(5)(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(1)(a) + h3/3! *
f(3)(a)
et donc:
f(1)(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(2)(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:
double Fonction(double x)
{
return sin(x);
}
double DeriveOrdre1P3(double a, double h,
double(*fonction)(double x))
{
double x;
x = (fonction(a+h) - fonction(a-h))/(2*h);
return x;
}
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:
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(1)(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.
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. 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.
Le programme Derivee.cpp, téléchargeable dans la bibliothèque de codes de TangenteX, permet, sur une fonction déclarée dans la routine Fonction, de comparer la précision des différentes méthodes.
Le package scipy de Python propose une routine de calcul de
dérivée d'ordre n en un point, construite sur l'algorithme des
différences centrées. Il s'agit de la routine derivative()
que vous trouverez dans scipy.misc.
Les principaux arguments sont : le nom de la fonction à dériver, le point auquel vous voulez la dériver, le pas de calcul et l'ordre de la dérivation. Sur ce dernier point, assurez-vous que votre fonction est bien de classe Cn, sinon vous aurez droit à un beau message d'erreur. Attention au pas de calcul dx, qui doit être assez petit, le "assez" dépendant de la fonction. Faites plusieurs essais...
Vous trouverez la documentation de la fonction derivative()
sur
cette page .
Commençons par un exemple très simple : imaginons que je veuille calculer la dérivée de la fonction f(x) = sin(x) sur l'intervalle de 0 à 2*PI. Cette fonction est continue et dérivable sur le domaine, donc pas de problème.
Comme nous l'avons vu ci-dessus, la fonction derivative()
calcule la dérivée de f(x) en un point. Il faut donc découper le
domaine en un grand nombre de points, calculer la dérivée en
chaque point et construire la courbe dérivée à l'aide de ces
calculs. Pas très pratique...
Heureusement, Python qui sait travailler avec des vecteurs nous
permet de réaliser ce calcul très simplement. Voyons comment dans
le script BODerivationSinus.py, téléchargeable dans la bibliothèque
de codes de TangenteX .
Pour commencer les habituels appels de librairie :
from scipy import sin, pi, arange
from scipy.misc import derivative
from matplotlib.pyplot import *
La fonction derivative()
est dans le package
scipy.misc.
Je construit maintenant le vecteur x, à l'aide de la fonction
arange()
. Puis je calcule le vecteur y = f(x):
x = arange(0,2*pi,0.01)
y = f(x)
Et maintenant l'objet de la manoeuvre : je vais dériver la
fonction f(x) pour tous les points du domaine et construire le
vecteur yprime qui contiendra ces dérivées. J'appelle la fonction
derivative()
, en lui indiquant que je souhaite une
dérivée première (le n = 1) et que le pas de dérivation est égal à
0.01 (le dx = 0.01). C'est bien sur le même que pour arange()
.
Cela donne :
yprime = derivative(f,x,n=1,dx=0.01)
Ensuite viennent les instructions de tracé de courbe :
Vous reconnaissez bien sur la courbe représentative de cos(x) en
la courbe de la fonction dérivée.
Le calcul des dérivées intervient en optimisation lorsqu’il s’agit de rechercher le minimum d’une fonction. Il est alors associé au calcul des racines d’une équation que nous avons vu plus haut. Voyons un exemple simple et très classique. Une usine fabrique des boites de conserve de 1 l en métal et voudrait économiser le plus possible sur ses coûts de fabrication. On va supposer que les coûts de procédé ont été optimisés et qu’il ne reste plus qu’à travailler sur les coûts de matières premières. L’objectif est donc de réduire au mieux la surface de métal nécessaire pour englober un volume d’un litre, en supposant l’épaisseur du métal constante.
Appelons h la hauteur de la boite et r son rayon. La surface du cylindre que constitue la boite s’exprime donc par \(S = 2 \pi r h + 2\pi r^2 \) et son volume \( \pi r^2h \), si j’exprime les grandeurs en décimètres. Rappelons que 1 litre = 1dm^3 et donc que \( \pi r^2 h = 1\). Je me sers de cette dernière égalité pour exprimer S en fonction de r uniquement, ce qui me donne \(S = \dfrac{2}{r} + 2\pi r^2 \).
Il me reste à calculer la dérivée de cette fonction, et à
chercher sa racine. Le script BODeriveeOpt.py, téléchargeable dans
la bibliothèque de codes
de TangenteX va le faire pour moi.
Après avoir défini la fonction de calcul de la surface du cylindre et le domaine de variation du rayon (attention à éliminer le rayon nul), je calcule la dérivée de la fonction Surface et je cherche le rayon ropt pour lequel elle s'annule.
Le programme trace les deux courbes S et S', ce qui nous permet
de lire le rayon optimum, qui est le même que celui qu'affiche le
programme, soit 0,54 dm.
Voici les deux courbes :