"# Calcul de dérivées et équation différentielles\n",
"\n",
"- Date : 2023-2024\n",
"- Durée approximative : 2 x 1h15\n",
"\n",
"## Introduction\n",
"\n",
"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 utilisation 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 que nous allons appliquer au calcul des équations variationnelles (ou linéarisées) des équations différentielles.\n",
"\n",
"On notera $\\|\\cdot\\|$ 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": [
"# activate local project\n",
"using Pkg\n",
"Pkg.activate(\".\")\n",
"\n",
"# load packages\n",
"using DualNumbers\n",
"using DifferentialEquations\n",
"using ForwardDiff\n",
"using LinearAlgebra\n",
"using Plots\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()` (ou `eps(Float64)` pour les flottants codés sur 64 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",
" f'(x) \\cdot v = g'(0) \\approx \\frac{\\mathrm{Im}(g(x+ih))}{h}.\n",
"$$\n",
"\n",
"On peut prouver que l'approximation ci-dessus est au moins d'ordre 2 et aussi démontrer que tout pas inférieur à $h_*$ est optimal, avec \n",
"$$\n",
" h_{*} \\approx \\sqrt{eps_{mach}}.\n",
"$$\n",
"\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",
" error(\"Choose a valid method in \", methods)\n",
" end\n",
" if method == :forward\n",
" return f(x) # TO UPDATE\n",
" elseif method == :central\n",
" return f(x) # TO UPDATE\n",
" elseif method == :complex\n",
" return f(x) # TO UPDATE\n",
" elseif method == :dual \n",
" return f(x) # TO UPDATE\n",
" elseif method == :forward_ad\n",
" if x isa Number\n",
" return ForwardDiff.derivative(f, x)*v\n",
" else\n",
" return ForwardDiff.jacobian(f, x)*v\n",
" end\n",
" end\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercice 1\n",
"\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": [
"# function to print derivative values and errors\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."
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
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()` (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:
## 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
f'(x) \cdot v = g'(0) \approx \frac{\mathrm{Im}(g(x+ih))}{h}.
$$
On peut prouver que l'approximation ci-dessus est au moins d'ordre 2 et aussi démontrer que tout pas inférieur à $h_*$ est optimal, avec
$$
h_{*} \approx \sqrt{eps_{mach}}.
$$
**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
%% 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: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
# 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: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.
%% Cell type:code id: tags:
``` julia
# affichage des erreurs en fonction de h
function print_errors(steps,errors,h_star,title)
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