Il existe plusieurs façon de calculer une dérivée sur un calculateur :
- par différences finies;
- par différentiation complexe;
- en utilisation la différentiation automatique;
- en utilisant le calcul formel et un générateur de code.
Nous allons étudier ici quelques cas que nous allons appliquer au calcul des équations variationnelles (ou linéarisées) des équations différentielles.
On notera $\|\cdot\|$ la norme euclidienne usuelle et $\mathcal{N}(0,1)$ une variable aléatoire Gaussienne centrée réduite.
%% Cell type:code id: tags:
``` julia
# activate local project
usingPkg
Pkg.activate(".")
# load packages
usingDualNumbers
usingDifferentialEquations
usingForwardDiff
usingLinearAlgebra
usingPlots
usingPlots.Measures
usingPrintf
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Dérivées par différences finies avant
</div>
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()` (ou `eps(Float64)` pour les flottants codés sur 64 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:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Dérivées par différentiation complexe
</div>
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:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Dérivées par différentiation automatique via les nombres duaux
</div>
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 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:
<div style="width:95%;
margin:10px;
padding:8px;
background-color:#afa;
border:2px solid #bbffbb;
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Pas optimal
</div>
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
