Calcul de Dérivées via les nombres duaux

Considérons le cas simple de la fonction \(f(x)=x^2\). La dérivée est trivialement donnée par \(f^\prime(x)=2x\). Choisissons quelques valeurs:

On voit clairement que cette approche n'est pas stable et l'erreur varie fortement en fonction des valeurs choisies. Cette erreur s'explique principalement par le développement de Taylor.

Cette approche pourrait être améliorée en utilisant des ordres plus élevés du développement de Taylor. Nous choisissons une autre méthode.

Revenons aux principes élémentaires de dérivation et aux 4 opérations arithmétiques relatives à la dérivation:

Définissons un nouveau type de nombre composé de 2 réels. Soit \(x\) et \(y\) deux de ces nouveaux nombres qui s'écrivent,

et définissons comment ajouter, multiplier, soustraire et diviser ses nombres,

Nous notons ce nombre \((\texttt{value}, \texttt{deriv})\).

1. La partie value se comporte comme les nombres réels avec les opérations arithmétiques usuelles.
2. La partie deriv se comporte avec les opérations arithmétiques de la dérivation \eqref{org8e45ca5}. Intéressant!

Un nombre complexe peut s'écrire \(z = x + iy\) et le nombre imaginaire \(i\) est défini tel que \(i^2=-1\). À partir de là, en utilisant les règles élémentaires de calculs et les opérations arithmétiques usuelles, il est facile d'obtenir par exemple,

ce que nous écririons:

De manière similaire, un nombre dual peut s'écrire \(z = x + \varepsilon y\) et \(\varepsilon\) est défini tel que \(\varepsilon^2=0\). Par conséquent, il est facile d'obtenir \eqref{org567bf64}, par exemple,

Reprenons l'exemple \(f(x)=x^2\) et définissons le nombre dual \(z=(a, 1.)\). Appliquons les opérations avec \(z\) au lieu de \(x\) dans \(f\), cela donne,

Donc nous avons évalué la fonction en la valeur \(a\) et nous avons aussi obtenu l'évaluation de la dérivée en cette valeur.

Pourquoi le 1. dans \((a, 1.)\) ? Si nous considérons la fonction \(f(x)=x\) alors en substituant \(z\) à la place de \(x\) dans \(f\), nous obtenons bien que la dérivée vaut 1.

Avec le même argument, une constante \(c\) se définit comme \((c, 0.)\) puisque sa dérivée est nulle.

Il faut donc distinguer ce qui représente une variable au sens mathématique comme \(x\) et une constante au sens mathétmatique comme \(c\). La dérivation se fait par rapport à cette variable \(x\).

Il apparaît assez clair que nous allons représenter un nombre dual avec struct. Le plus simple est:

typedef struct {
    char name[MAXCHAR];
    float value;
    float deriv;
} dual;

Comme attendu, la structure possède 2 champs: l'un pour représenter les valeurs value et l'autre pour représenter la dérivée deriv. Nous ajoutons un champ name pour associer un nom en espérant faire de plus jolis affichages.

Pour tester, nous écrivons le corps de notre programme et l'affichage que nous obtenons.

/* provide printf/scanf and other */
#include <stdio.h>

/* Length maximum for some "string" (array of char) */
#define MAXCHAR 1000

typedef struct {
    char name[MAXCHAR];
    float value;
    float deriv;
} dual;

int main () {
    printf("Hi");

    /* the value we are interrested in */
    float a = 2.1;

    /* declare a dual number */
    dual z;

    /* fill the fields                                 */
    /*                                                 */
    /* sprintf is similar to printf                    */
    /*  instead of output to the console               */
    /*  sprintf store the output to its first argument */
    /*     (here z.name)                               */
    sprintf(z.name, "%s", "LE nom bien");
    z.value = a;
    z.deriv = 1.;

    /* display the dual number */
    printf("\n");
    printf("\tName : %s\n", z.name);
    printf("\tValue: %f\n", z.value);
    printf("\tDeriv: %f\n", z.deriv);

    printf("Bye.");
    return 0;
  }

Hi
	Name : LE nom bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

Nous souhaitons avoir une fonction qui prend en argument un nombre réel (float) et une chaîne de caractère (char*) et retourne un nombre dual (return z; avec la variable z de type dual). Nous appelons cette fonction newvar et sa signature est donc:

dual newvar(float, char*);

Son implémentation est aussi directe puisque nous avions déjà écrit cela dans le programme précédent.

dual newvar(float val, char* avoile) {
    dual z;
    sprintf(z.name, "%s", avoile);
    z.value = val;
    z.deriv = 1.;
    return z;
}

WARNING: Les quantités val et z sont des variables au sens informatique, et respectivement de type float et dual. Cependant, la variable informatique val ne modélise pas la notion de variable au sens mathématique, c'est-à-dire qu'au nombre val n'est pas associé une dérivée. Nous construisons justement un code pour cela et le nouveau type dual sera une modélisation informatique plus proche de la notion mathématique de variable en considérant la dérivation.

Ainsi nous adaptons notre programme pour obtenir le même affichage que précédemment.

Hi
	Name : LE nome bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

Bien évidemment, nous pouvons aussi définir une fonction similaire à newvar qui crée une constante, disons newcst. Laissé en exercice, sinon regardez dans le fichier source.

Hi
	Name : LE nome bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

Nous souhaitons avoir une fonction qui prend en argument un nombre dual et ne retourne rien (void). Sa signature est,

void print(dual);

Et nous l'avons déjà implémentée.

void print(dual z) {
    printf("\n");
    printf("\tName : %s\n", z.name);
    printf("\tValue: %f\n", z.value);
    printf("\tDeriv: %f\n", z.deriv);
}

Hi
	Name : LE nome bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

Au lieu d'afficher à chaque fois le programme complet, nous allons maintenant—pour des raisons de lisibilité—afficher uniquement les modifications que nous avons apportées.

--- Supplementary/main-v2b.c    2024-11-01 11:31:04.501237581 +0000
+++ Supplementary/main-v3.c     2024-11-01 11:31:04.502237580 +0000
@@ -14,6 +14,7 @@
       /* defined after the main    */
 dual newvar(float, char*);
 dual newcst(float, char*);
+void print(dual);

 int main () {
     printf("Hi");
@@ -29,10 +30,7 @@
     z = newvar(a, "LE nome bien");

     /* display the dual number */
-    printf("\n");
-    printf("\tName : %s\n", z.name);
-    printf("\tValue: %f\n", z.value);
-    printf("\tDeriv: %f\n", z.deriv);
+    print(z);

     printf("Bye.");
     return 0;
@@ -56,3 +54,10 @@
     z.deriv = 0.;
     return z;
 }
+
+void print(dual z) {
+    printf("\n");
+    printf("\tName : %s\n", z.name);
+    printf("\tValue: %f\n", z.value);
+    printf("\tDeriv: %f\n", z.deriv);
+}

La signature de la fonction qui permet la multiplication entre deux nombres duaux est assez naturelle. Cette fonction a 2 arguments qui sont des nombres duaux (dual) et elle retourne un nombre dual (return z; avec z une variable de type dual).

dual mul(dual, dual);

Et l'implémentation est donnée par \eqref{org567bf64}.

dual mul(dual z1, dual z2) {
    dual z;
    sprintf(z.name, "(%s*%s)", z1.name, z2.name);

    z.value = z1.value * z2.value;
    z.deriv = (z1.deriv * z2.value) + (z1.value * z2.deriv);

    return z;
}

Pour tester cette nouvelle fonction, nous écrivons le programme principal ci-dessous, que nous exécutons et nous obtenons bien le résultat attendu.

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : (x*x)
	Value: 4.409999
	Deriv: 4.200000
Bye.

Nous obtenons bien pour la fonction \(x^2\) évaluée en 2.1 la valeur \(2.1^2=\) 4.41 (à l'arrondi près) et sa dérivée \(2\times 2.1=\) 4.2.

--- Supplementary/main-v3.c     2024-11-01 11:31:04.502237580 +0000
+++ Supplementary/main-v4a.c    2024-11-01 11:31:04.503237579 +0000
@@ -16,6 +16,8 @@
 dual newcst(float, char*);
 void print(dual);

+dual mul(dual, dual);
+
 int main () {
     printf("Hi");

@@ -27,11 +29,14 @@


     /* fill the fields */
-    z = newvar(a, "LE nome bien");
+    z = newvar(a, "x");

     /* display the dual number */
     print(z);

+    dual zz = mul(z, z);
+    print(zz);
+
     printf("Bye.");
     return 0;
   }
@@ -61,3 +66,13 @@
     printf("\tValue: %f\n", z.value);
     printf("\tDeriv: %f\n", z.deriv);
 }
+
+dual mul(dual z1, dual z2) {
+    dual z;
+    sprintf(z.name, "(%s*%s)", z1.name, z2.name);
+
+    z.value = z1.value * z2.value;
+    z.deriv = (z1.deriv * z2.value) + (z1.value * z2.deriv);
+
+    return z;
+}

Sur le même principe, il est facile d'implémenter les opérations manquantes add, sub et div. Laissé en exercice, sinon regardez dans le fichier source.

Nous commençons pour définir une valeur variable (newvar). C'est une variable informatique arbitrairement nommée z et d'autre part nous lui attribuons le nom "x" dans notre modélisation des nombres duaux.

Nous définissons aussi une valeur constante (newcst). C'est une variable informatique arbitrairement nommée c et d'autre part nous lui attribuons le nom "c" dans notre modélisation des nombres duaux.

Puis nous définissons un programme principal qui affiche:

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : c
	Value: 1.200000
	Deriv: 0.000000

	Name : (x*c)
	Value: 2.520000
	Deriv: 1.200000

	Name : (x+x)
	Value: 4.200000
	Deriv: 2.000000

	Name : (c-x)
	Value: -0.900000
	Deriv: -1.000000

	Name : ((x+c)/(c-x))
	Value: -3.666667
	Deriv: 2.962964
Bye.

À partir de ces exemples, nous pouvons tester si notre implémentation est correcte. Il faut faire plusieurs tests pour s'assurer que tout fonctionne comme espéré.

Nous commençons à voir l'intérêt d'avoir un champ name dans notre modélisation des nombres duaux.

Il est assez clair que notre fonction exponential aura comme argument un nombre dual et retournera un nombre dual, donc sa signature est,

#define NEXP 20
dual exponential(dual);

Et son implémentation est simple mais mérite un peu d'attention. La somme \(\sum\) va être naturellement transformée en boucle for. Par ailleurs, cette somme commence par k=0 mais comme nous allons accumuler son résultat, nous allons initialiser l'accumulateur (variable informatique expx) à la valeur k=0 et faire commencer la boucle à k=1.

dual exponential(dual x) {

    /* temporary variables */
    dual expx  = newcst(1., "");
    dual xk    = newcst(1., "");
    dual kfac  = newcst(1., "");

    float k;
    for (k=1; k<=NEXP; k++) {
        dual kk = newcst(k, "");

        xk   = mul(xk, x);      /* x, x*x, x*x*x, etc. => eval x^k */
        kfac = mul(kfac, kk);   /* 1, 1*2, 1*2*3, etc. => eval k!  */

        dual term_k = div(xk, kfac);

        /* update the Sum with another term */
        expx = add(expx, term_k);

        /* to avoid segmentation fault (because of ugly C)   */
        /* instead of accumalating all the ( and ) and +,*,/ */
        sprintf(expx.name, "");
        sprintf(xk.name, "");
        sprintf(kfac.name, "");
        sprintf(kk.name, "");
    }

    /* pretty print the name */
    sprintf(expx.name, "~(e^(%s))", x.name);
    return expx;
}

Avec un papier et un crayon, et en écrivant les premiers termes, il est rapidement clair pourquoi expx, xk et kfrac sont initialisés avec la fonction newcst et non la fonction newvar.

Finalement, nous écrivons un programme principal pour tester.

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : c
	Value: 1.200000
	Deriv: 0.000000

	Name : ~(e^(x))
	Value: 8.166168
	Deriv: 8.166168

	Name : ~(e^((((c*x)*x)+c)))
	Value: 659.838440
	Deriv: 3325.548828
Bye.

Il est clair que nous obtenons bien le bon résultat, même pour une fonction non-triviale. Notre implémentation n'utilise que les 4 opérations arithmétiques.

Nous pourrions facilement étendre cela à toutes les fonctions (sinus, cosinus, tangente, etc.). Par exemple, nous pourrions calculer la fonction sinus en remarquant que le sinus est la partie imaginaire d'une fonction exponentielle complexe. Par conséquent, nous aurions besoin d'implémenter une structure représentant les nombres complexes puis d'implémenter les 4 opérations arithmétiques pour ce nouveau type complex. Et enfin modifier notre structure dual en remplaçant le type float par ce type complex.

Nous avons choisi la fonction exponentielle mais nous aurions aussi pu choisir la fonction racine carrée comme illustration. Elle se calcule par la méthode de Héron comme,

Pour illustrer notre propos, étendons la fonction exp de math.h. Ou pour le dire autrement, nous allons utiliser la fonction exp de math.h avec notre type dual.

En sachant que la fonction exp de la bibliothèque math.h a la signature:

double exp(double);

or double et float sont des types compatibles. Nous voulons une signature similaire pour le type dual ; en fait la même signature que la fonction exponential:

#include <math.h>
dual dexp(dual);

Nous choisissons arbitraiment le nom dexp, et il nous semble parlant: l'extension de exp au type dual.

Pour l'implémentation, il faut commencer par se rappeler la définition mathématique:

Ce qui se traduit naturellement par:

dual dexp(dual u) {
    dual z;
    sprintf(z.name, "(e^(%s))", u.name);

    z.value = exp(u.value);
    z.deriv = u.deriv * z.value;

    return z;
}

Il faut noter que la multiplication u.deriv * z.value se fait dans le type float. Finalement, nous comparons nos 2 versions.

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : c
	Value: 1.200000
	Deriv: 0.000000

# by +,*,-,/

	Name : ~(e^((((c*x)*x)+c)))
	Value: 659.838440
	Deriv: 3325.548828

# by math.h

	Name : (e^((((c*x)*x)+c)))
	Value: 659.841492
	Deriv: 3325.601074
Bye.

Sur ce même principe, nous pouvons étendre toutes les fonctions qui sont dans la bibliothèque math.h et dont nous connaissons les formules de dérivations, comme sinus, cosinus, tangente, etc.

Avec notre code, nous sommes maintenant capable de calculer la dérivée de n'importe quelque fonction. Même des fonctions qui sont impossibles à dériver à la main. Par exemple,

dual fweird(dual x) {
    dual y;
    sprintf(y.name, "%s", x.name);
    y.deriv = x.deriv;
    if (x.value < 0) {
        return x;
    } else if (((int)x.value) % 2 == 0) { /* even */
        y.value = x.value - 1.;
        return add(div(sub(x, newcst(1, "1")),
                       add(x, newcst(2, "2"))),
                   fweird(y));
    } else {                              /* odd   */
        y.value = x.value - 1.;
        return mul(x, fweird(y));
    }
}

Hi
	Name : x
	Value: 12.300000
	Deriv: 1.000000

	Name : (((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+x)))))))))))))
	Value: -10758.750000
	Deriv: 6566.896484
Bye.

Pour connaître la taille mémoire utilisée pour représenter un nombre entier en mémoire, le plus simple est d'utiliser sizeof(int) (retourne en octet/byte). En général, ceci est de 4 octets soit 32 bits. Seul les entiers dans l'intervalle \(\left[-2^{31} ~;~ 2^{31} - 1\right]\) sont représentables dans l'ordinateur. Le plus grand entier est donc 2147483647. Il est possible de représenter un nombre 2 fois plus grand en considérant le type unsigned c'est-à-dire uniquement les entiers positifs. Pour des nombres encore plus grands, il faut utiliser le type long et unsigned long.

Mais qu'en est-il du type float ?

Ceci dépend ! La virgule est flottante et donc il y a une précision variable.

#include <stdio.h>
#include <math.h>
#include <limits.h>
int main () {
    float x=1., f=1., ff;
    while (!isinf(f)) {
        ff = f;

        f = f * x;

        x = x + 1.;
        printf("%.0f!  \t= %f\n", x-1, f);
    }
    return 0;
}

1!  	= 1.000000
2!  	= 2.000000
3!  	= 6.000000
4!  	= 24.000000
5!  	= 120.000000
6!  	= 720.000000
7!  	= 5040.000000
8!  	= 40320.000000
9!  	= 362880.000000
10!  	= 3628800.000000
11!  	= 39916800.000000
12!  	= 479001600.000000
13!  	= 6227020800.000000
14!  	= 87178289152.000000
15!  	= 1307674279936.000000
16!  	= 20922788478976.000000
17!  	= 355687414628352.000000
18!  	= 6402373530419200.000000
19!  	= 121645096004222976.000000
20!  	= 2432902023163674624.000000
21!  	= 51090940837169725440.000000
22!  	= 1124000724806013026304.000000
23!  	= 25852017444594485559296.000000
24!  	= 620448454699064672387072.000000
25!  	= 15511211079246240657965056.000000
26!  	= 403291499589617303175561216.000000
27!  	= 10888870415132690890901946368.000000
28!  	= 304888371623715344945254498304.000000
29!  	= 8841763079319199907069674127360.000000
30!  	= 265252889961724357982831874408448.000000
31!  	= 8222839685527520666638122083155968.000000
32!  	= 263130869936880661332419906660990976.000000
33!  	= 8683318509846655538309012935952826368.000000
34!  	= 295232822996533287161359432338880069632.000000
35!  	= inf

À comparer avec la séquence exacte:

1!      = 1
2!      = 2
3!      = 6
4!      = 24
5!      = 120
6!      = 720
7!      = 5040
8!      = 40320
9!      = 362880
10!     = 3628800
11!     = 39916800
12!     = 479001600
13!     = 6227020800
14!     = 87178291200
15!     = 1307674368000
16!     = 20922789888000
17!     = 355687428096000
18!     = 6402373705728000
19!     = 121645100408832000
20!     = 2432902008176640000
21!     = 51090942171709440000
22!     = 1124000727777607680000
23!     = 25852016738884976640000
24!     = 620448401733239439360000
25!     = 15511210043330985984000000
26!     = 403291461126605635584000000
27!     = 10888869450418352160768000000
28!     = 304888344611713860501504000000
29!     = 8841761993739701954543616000000
30!     = 265252859812191058636308480000000
31!     = 8222838654177922817725562880000000
32!     = 263130836933693530167218012160000000
33!     = 8683317618811886495518194401280000000
34!     = 295232799039604140847618609643520000000
35!     = 10333147966386144929666651337523200000000
36!     = 371993326789901217467999448150835200000000
37!     = 13763753091226345046315979581580902400000000

Donc nous voyons que la représentation en nombre flottant fait un calcul faux pour factorielle 14 (14!), et ensuite les erreurs se cumulent. Ceci est attendu et la précision dépend de ce que l'on cherche à calculer. Pour les calculs standards, nous sommes très rarement dans ce cas. Ici, nous testons les limites.

Par ailleurs, nous voyons que nous ne pouvons pas utiliser plus de 36 termes dans notre calcul de l'exponentielle utilisant la série.