"Il existe plusieurs façon de calculer une dérivée sur un calculateur : \n",
"\n",
"- par différences finies;\n",
"- par différentiation complexe;\n",
"- en utilisant la différentiation automatique;\n",
"- en utilisant le calcul formel et un générateur de code.\n",
"\n",
"Nous allons étudier ici quelques cas.\n",
"\n",
"On notera $\\Vert{\\cdot}\\Vert$ la norme euclidienne usuelle et $\\mathcal{N}(0,1)$ une variable aléatoire Gaussienne centrée réduite."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load packages\n",
"using DualNumbers\n",
"using DifferentialEquations\n",
"using ForwardDiff\n",
"using LinearAlgebra\n",
"using Plots\n",
"using Plots.Measures\n",
"using Printf"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dérivées par différences finies avant\n",
"\n",
"Soit $f$ une fonction lisse de $\\mathbb{R}^{n}$ dans $\\mathbb{R}^{m}$, $x$ un point de $\\mathbb{R}^{n}$ et $v$ un vecteur de $\\mathbb{R}^{n}$. Si on note $g \\colon h \\mapsto f(x+hv)$, alors on a d'après la formule de Taylor-Young :\n",
"L'approximation ainsi obtenue de $ g'(0) = f'(x) \\cdot v \\in \\mathbb{R}^m$ est d'ordre 1 si $g^{(2)}(0) \\neq 0$ ou au moins d'ordre 2 sinon. \n",
"\n",
"**Remarque.** Sur machine, Les calculs se font en virgule flottante. On note epsilon machine, le plus petit nombre $\\mathrm{eps}_\\mathrm{mach}$ tel que $1+\\mathrm{eps}_\\mathrm{mach}\\ne 1$. Cette quantité dépend des machines et de l'encodage des données. Pour l'optenir en `Julia` il suffit de taper `eps()`. On peut indiquer la précision numérique : `eps(Float64)` ou `eps(Float64(1))` pour les flottants codés sur 64 bits et `eps(Float32(1))` sur 32 bits.\n",
"\n",
"Notons $\\mathrm{num}(g,\\, h)$ la valeur de $g(h)$ calculée numériquement et supposons que l'on puisse majorer l'erreur relative numérique par : \n",
"l'approximation ainsi obtenue de $f'(x) \\cdot v \\in \\mathbb{R}^{m}$ est d'ordre 2 si $g^{(3)}(0) \\neq 0$ ou au moins d'ordre 4 sinon. À noter que ce schéma nécessite plus d'évaluations de la fonction $f$. On peut montrer comme précédemment que le meilleur $h$ est de l'ordre \n",
"Les formules des schémas avant et centrée sont sensibles aux calculs de la différence $\\Delta f = f(x+h) - f(x)$ ou $\\Delta f = f(x+h) - f(x-h)$. Pour remédier à ce problème, les [différences finies à pas complexe](https://dl.acm.org/doi/10.1145/838250.838251) ont été introduites.\n",
"Si on suppose que la fonction $g$ est holomorphe, c'est-à-dire dérivable au sens complexe,\n",
"on peut considérer un pas complexe $ih$. Un développement limité de $g$ en $0$ s'écrit\n",
"**Remarque.** Utiliser en `Julia` la commande `imag` pour calculer la partie imaginaire d'un nombre complexe et la variable `im` pour représenter l'unité imaginaire."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dérivées par différentiation automatique via les nombres duaux\n",
"\n",
"Les nombres duaux s'écrivent sous la forme $a + b\\, \\varepsilon$ avec $(a,b)\\in \\mathbb{R}^2$ et $\\varepsilon^2 = 0$. Nous allons voir comment nous pouvons les utiliser pour calculer des dérivées.\n",
"\n",
"Soit deux fonctions $f, g \\colon \\mathbb{R} \\to \\mathbb{R}$ dérivables, de dérivées respectives $f'$ et $g'$. On pose\n",
"1. Compléter la fonction `derivative` ci-dessus avec les méthodes de différences finies avant, centrée, par différentiation complexe et par différentiation automatique via les nombres duaux.\n",
"2. Exécuter le code ci-dessous et vérifier les résultats obtenus."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Scalar case\n",
"\n",
"# check if the derivatives are correct\n",
"f(x) = cos(x)\n",
"x = π/4\n",
"v = 1.0\n",
"\n",
"# solution\n",
"sol = -sin(x)*v\n",
"\n",
"# print derivatives and errors for each method\n",
"# print derivatives and errors for each method\n",
"print_derivatives(f, x, v, sol)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pas optimal\n",
"\n",
"On se propose de tester pour la fonction $\\cos$ aux points $x_0=\\pi/3$, $x_1 = 10^6\\times\\pi/3$ et la fonction $\\cos+10^{-8} \\mathcal{N}(0, \\, 1)$ au point $x_0=\\pi/3$ l'erreur entre les différences finies et la dérivée au point considéré en fonction de $h$. On prendra $h=10^{-i}$ pour $i= \\{1,\\ldots,16\\}$ et on tracera ces erreurs dans une échelle logarithmique (en `Julia`, avec le package `Plots` on utilise l'option `scale=:log10`).\n",
"\n",
"## Exercice 2\n",
"\n",
"- Visualiser les différentes erreurs en fonction de $h$ pour les différentes méthodes de calcul de dérivées. Commentaires.\n",
"- Modifier la précision de $x_0$ et $x_1$ en `Float32`. Commentaires."
Il existe plusieurs façon de calculer une dérivée sur un calculateur :
- par différences finies;
- par différentiation complexe;
- en utilisant la différentiation automatique;
- en utilisant le calcul formel et un générateur de code.
Nous allons étudier ici quelques cas.
On notera $\Vert{\cdot}\Vert$ la norme euclidienne usuelle et $\mathcal{N}(0,1)$ une variable aléatoire Gaussienne centrée réduite.
%% Cell type:code id: tags:
``` julia
# load packages
usingDualNumbers
usingDifferentialEquations
usingForwardDiff
usingLinearAlgebra
usingPlots
usingPlots.Measures
usingPrintf
```
%% Cell type:markdown id: tags:
## Dérivées par différences finies avant
Soit $f$ une fonction lisse de $\mathbb{R}^{n}$ dans $\mathbb{R}^{m}$, $x$ un point de $\mathbb{R}^{n}$ et $v$ un vecteur de $\mathbb{R}^{n}$. Si on note $g \colon h \mapsto f(x+hv)$, alors on a d'après la formule de Taylor-Young :
L'approximation ainsi obtenue de $ g'(0) = f'(x) \cdot v \in \mathbb{R}^m$ est d'ordre 1 si $g^{(2)}(0) \neq 0$ ou au moins d'ordre 2 sinon.
**Remarque.** Sur machine, Les calculs se font en virgule flottante. On note epsilon machine, le plus petit nombre $\mathrm{eps}_\mathrm{mach}$ tel que $1+\mathrm{eps}_\mathrm{mach}\ne 1$. Cette quantité dépend des machines et de l'encodage des données. Pour l'optenir en `Julia` il suffit de taper `eps()`. On peut indiquer la précision numérique : `eps(Float64)` ou `eps(Float64(1))` pour les flottants codés sur 64 bits et `eps(Float32(1))` sur 32 bits.
Notons $\mathrm{num}(g,\, h)$ la valeur de $g(h)$ calculée numériquement et supposons que l'on puisse majorer l'erreur relative numérique par :
l'approximation ainsi obtenue de $f'(x) \cdot v \in \mathbb{R}^{m}$ est d'ordre 2 si $g^{(3)}(0) \neq 0$ ou au moins d'ordre 4 sinon. À noter que ce schéma nécessite plus d'évaluations de la fonction $f$. On peut montrer comme précédemment que le meilleur $h$ est de l'ordre
$$
h_*\approx \sqrt[3]{\mathrm{eps}_\mathrm{mach}}.
$$
%% Cell type:markdown id: tags:
## Dérivées par différentiation complexe
Les formules des schémas avant et centrée sont sensibles aux calculs de la différence $\Delta f = f(x+h) - f(x)$ ou $\Delta f = f(x+h) - f(x-h)$. Pour remédier à ce problème, les [différences finies à pas complexe](https://dl.acm.org/doi/10.1145/838250.838251) ont été introduites.
Si on suppose que la fonction $g$ est holomorphe, c'est-à-dire dérivable au sens complexe,
on peut considérer un pas complexe $ih$. Un développement limité de $g$ en $0$ s'écrit
**Remarque.** Utiliser en `Julia` la commande `imag` pour calculer la partie imaginaire d'un nombre complexe et la variable `im` pour représenter l'unité imaginaire.
%% Cell type:markdown id: tags:
## Dérivées par différentiation automatique via les nombres duaux
Les nombres duaux s'écrivent sous la forme $a + b\,\varepsilon$ avec $(a,b)\in \mathbb{R}^2$ et $\varepsilon^2 = 0$. Nous allons voir comment nous pouvons les utiliser pour calculer des dérivées.
Soit deux fonctions $f, g \colon \mathbb{R} \to \mathbb{R}$ dérivables, de dérivées respectives $f'$ et $g'$. On pose
## La méthode principale pour le calcul de dérivées
La fonction `derivative` ci-dessous calcule la dérivée directionnelle
$$
f'(x) \cdot v.
$$
%% Cell type:code id: tags:
``` julia
## TO COMPLETE
# compute directional derivative
function derivative(f,x,v;method=_method(),h=_step(x,v,method))
ifmethod∉methods
error("Choose a valid method in ",methods)
end
ifmethod==:forward
returnf(x)# TO UPDATE
elseifmethod==:central
returnf(x)# TO UPDATE
elseifmethod==:complex
returnf(x)# TO UPDATE
elseifmethod==:dual
returnf(x)# TO UPDATE
elseifmethod==:forward_ad
ifxisaNumber
returnForwardDiff.derivative(f,x)*v
else
returnForwardDiff.jacobian(f,x)*v
end
end
end;
```
%% Cell type:code id: tags:
``` julia
# function to print derivative values and errors
function print_derivatives(f,x,v,sol)
println("Hand derivative: ",sol,"\n")
formethod∈methods
dfv=derivative(f,x,v,method=method)
println("Method: ",method)
println(" derivative: ",dfv)
@printf(" error: %e\n",norm(dfv-sol))
ifmethod∈(:forward,:central,:complex)
step=_step(x,v,method)
println(" step: ",step)
@printf(" error/step: %e\n",norm(dfv-sol)/step)
end
println()
end
end;
```
%% Cell type:markdown id: tags:
## Exercice 1
1. Compléter la fonction `derivative` ci-dessus avec les méthodes de différences finies avant, centrée, par différentiation complexe et par différentiation automatique via les nombres duaux.
2. Exécuter le code ci-dessous et vérifier les résultats obtenus.
%% Cell type:code id: tags:
``` julia
# Scalar case
# check if the derivatives are correct
f(x)=cos(x)
x=π/4
v=1.0
# solution
sol=-sin(x)*v
# print derivatives and errors for each method
print_derivatives(f,x,v,sol)
```
%% Cell type:code id: tags:
``` julia
# Vectorial case
# check if the derivatives are correct
f(x)=[0.5*(x[1]^2+x[2]^2);x[1]*x[2]]
x=[1.0,2.0]
v=[1.0,-1.0]
# solution
sol=[x[1]*v[1]+x[2]*v[2],x[1]*v[2]+x[2]*v[1]]
# print derivatives and errors for each method
print_derivatives(f,x,v,sol)
```
%% Cell type:markdown id: tags:
## Pas optimal
On se propose de tester pour la fonction $\cos$ aux points $x_0=\pi/3$, $x_1 = 10^6\times\pi/3$ et la fonction $\cos+10^{-8} \mathcal{N}(0, \, 1)$ au point $x_0=\pi/3$ l'erreur entre les différences finies et la dérivée au point considéré en fonction de $h$. On prendra $h=10^{-i}$ pour $i= \{1,\ldots,16\}$ et on tracera ces erreurs dans une échelle logarithmique (en `Julia`, avec le package `Plots` on utilise l'option `scale=:log10`).
## Exercice 2
- Visualiser les différentes erreurs en fonction de $h$ pour les différentes méthodes de calcul de dérivées. Commentaires.
- Modifier la précision de $x_0$ et $x_1$ en `Float32`. Commentaires.
%% Cell type:code id: tags:
``` julia
# affichage des erreurs en fonction de h
function plot_errors(steps,errors,h_star,title)
steps_save=steps
ymax=10^10
# supprimer les erreurs nulles
non_nul_element=findall(!iszero,errors)
errors=errors[non_nul_element]
steps=steps[non_nul_element]
# Courbe des erreurs pour les differents steps en bleu
