{
"cells": [
{
"cell_type": "markdown",
"id": "657ee43e-abd8-4efe-b5e0-6a632114529a",
"metadata": {
"tags": []
},
"source": [
"# Cours ENPC - Pratique du calcul scientifique"
]
},
{
"cell_type": "markdown",
"id": "a30eaeba",
"metadata": {},
"source": [
"## Résolution de systèmes non linéaires avec introduction à la différentiation automatique"
]
},
{
"cell_type": "markdown",
"id": "d3623bef",
"metadata": {},
"source": [
"--------------"
]
},
{
"cell_type": "markdown",
"id": "df935bca",
"metadata": {},
"source": [
"### Exercice sur Newton-Raphson en dimension 2\n",
"\n",
"On considère le système non linéaire suivant\n",
"$$\n",
"\\left \\{\n",
" \\begin{aligned}\n",
" &y = (x-1)^2 \\\\\n",
" &x^2 + y^2 = 4\n",
" \\end{aligned}\n",
"\\right.\n",
"$$\n",
"\n",
"1. En exploitant des tracés graphiques appropriés pour visualiser grossièrement la (les) zone(s) contenant une (des) solution(s), implémenter un algorithme de Newton-Raphson pour calculer une approximation de cette (ces) solution(s) en choississant un (des) point(s) initial(aux) judicieux.\n",
"\n",
"1. Modifier l'algorithme pour recueillir l'ensemble des itérations et visualiser l'évolution de l'erreur en fonction de l'itération pour diverses valeurs initiales en utilisant une échelle appropriée en y. On pourra estimer qu'une solution donnée est fournie par la valeur finale de l'algorithme."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0800dc52",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "8782e8fd",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "32580e74",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "a0f9ef06",
"metadata": {},
"source": [
"--------------------\n",
"--------------------\n",
"--------------------"
]
},
{
"cell_type": "markdown",
"id": "f9e89b43",
"metadata": {},
"source": [
"### Exercice sur la méthode babylonienne\n",
"\n",
"Soit un paramètre réel $a>0$ et la suite définie par\n",
"\n",
"$$\n",
"\\tag{1}\n",
"x_0>0 \\qquad ; \\qquad ∀k∈\\mathbb{N},\\quad x_{k+1}=\\frac{1}{2}\\left(x_k+\\frac{a}{x_k}\\right)\n",
"$$\n",
"\n",
"1. En écrivant $x_{k+1}-\\sqrt{a}$ en fonction de $x_k-\\sqrt{a}$ puis $x_{k+1}-x_k$, montrer que $(x_k)$ converge quadratiquement vers $x_*=\\sqrt{a}$ quel que soit $x_0>0$.\n",
"\n",
"1. Montrer que la formulation par récurrence (1) n'est autre que l'algorithme de Newton-Raphson appliqué à une fonction à identifier s'annulant en $x_*=\\sqrt{a}$.\n",
"\n",
"1. Construire une fonction `Babylonian` prenant comme arguments `a` ainsi qu'un entier `n` valant `10` par défaut et qui renvoie le vecteur $[x_0,x_1,\\ldots,x_n]$ en initialisant la suite à $x_0=\\frac{1+a}{2}$.\n",
"\n",
"1. Tester sur quelques valeurs.\n",
"\n",
"1. Tracer l'erreur $|x_k-x_*|$ en fonction du rang $k$ pour $a=2$.\n",
"\n",
" L'idée de la suite de l'exercice est d'appliquer la fonction `Babylonian` définie plus haut à un argument `a` non pas de type `Float64` mais d'un nouveau type permettant d'estimer à la fois la valeur de $\\sqrt{a}$ mais également de la dérivée de $a\\mapsto\\sqrt{a}$ soit $\\frac{1}{2\\sqrt{a}}$. Pour cela, on introduit de nouveaux nombres qu'on appelle **nombres duaux**. Ceux-ci sont définis à l'instar des nombres complexes à partir de la définition d'un nombre particulier noté $ɛ$ de sorte qu'un nombre dual s'écrit $x=a+bɛ$ avec $a$ et $b$ réels. En quelque sorte $ɛ$ joue un rôle analogue au $i$ complexe à la différence que l'on pose ici $ɛ^2=0$. L'objectif de tels nombres est d'être en mesure de stocker à la fois la valeur d'une fonction mais également sa dérivée en posant\n",
"\n",
" \n",
" $$\n",
" \\tag{2}\n",
" f(a+bɛ)=f(a)+f'(a)bɛ\n",
" $$\n",
"\n",
" ce qui entraîne que la dérivée en $a$ de $f$ peut être obtenue en récupérant la composante sur $ɛ$ de $f(a+ɛ)$ (i.e. en prenant $b=1$).\n",
"\n",
" En pratique, il est donc nécessaire de redéfinir le comportement des fonctions usuelles en cohérence avec (2). Toutefois dans l'application actuelle, seules les opérations `+`, `-`, `*` et `/` seront nécessaires et donc devront être surchargées pour pouvoir prendre comme arguments des nombres duaux. En outre, il sera également nécessaire d'implémenter les fonctions `convert` pour convertir un réel en nombre dual et `promote_rule` pour exprimer le fait qu'en présence d'une opération impliquant deux nombres dont l'un est dual, les deux doivent d'abord être exprimés sous forme de nombres duaux avant de lancer l'opération. À noter que la surcharge des opérateurs et fonctions de base n'est possible que si ceux-ci sont explicitement importés à l'aide par exemple de `import Base: +, -, ...`. Il est également possible de définir la fonction `Base.show` de manière à ce que l'affichage d'un nombre dual ait explicitement la forme `a+bɛ`.\n",
"\n",
" La surcharge des opérateurs s'exprime mathématiquement par\n",
" $$\n",
" \\begin{align*}\n",
" (a+bɛ)+(c+dɛ)&=(a+c)+(b+d)ɛ\\\\\n",
" (a+bɛ)-(c+dɛ)&=(a-c)+(b-d)ɛ\\\\\n",
" (a+bɛ)*(c+dɛ)&=ac+(bc+ad)ɛ\\\\\n",
" (a+bɛ)/(c+dɛ)&=\\frac{a}{c}+\\frac{bc-ad}{c^2}ɛ\n",
" \\end{align*}\n",
" $$\n",
"\n",
" Alternativement à cette dernière opération, on peut aussi définir $\\mathrm{inv}(a+bɛ)=\\frac{1}{a}-\\frac{b}{a^2}ɛ$ puis $u/v=u*\\mathrm{inv}(v)$.\n",
"\n",
"1. Étudier le `struct D` défini pour représenter un nombre dual ainsi que les lignes de code associées. Compléter les parties de code manquantes.\n",
"\n",
"1. Définir une instance du nombre `ɛ` (\\varepsilon puis touche TAB pour afficher ε), en d'autres termes le nombre `0+1ɛ` et effectuer quelques opérations pour vérifier les implémentations, par exemple\n",
"\n",
" ```julia\n",
" (1+2ɛ)*(3+4ɛ)\n",
" 1/(1+ɛ)\n",
" (1+2ɛ)/(2-ɛ)\n",
" ```\n",
"\n",
"1. Exploiter la structure de nombre dual pour estimer la dérivée de la fonction racine à partir de la méthode babylonienne (en exploitant directement la fonction `Babylonian` sans la réécrire) sur quelques exemples ($a=0.1$, $a=2$, $a=25$) et vérifier les résultats avec la solution analytique.\n",
"\n",
"1. Superposer sur un graphe la dérivée de la racine obtenue par la méthode babylonienne sur nombre dual et l'expression analytique.\n",
"\n",
"1. Proposer une méthode analogue pour calculer la racine $p^\\textrm{ème}$ d'un nombre $a$ i.e. $\\sqrt[p]{a}$. Vérifier que la dérivée de la racine $p^\\textrm{ème}$ peut également être obtenue par exploitation des nombres duaux sans ligne de code supplémentaire."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "05eb72f1",
"metadata": {},
"outputs": [],
"source": [
"function Babylonian(a; n = 10)\n",
" # Votre code ici\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e4f5e336",
"metadata": {},
"outputs": [],
"source": [
"for a in (0.1, 2, 25) \n",
" # Votre code de vérification de l'algorithme ici\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "095a9ed7",
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"\n",
"# Votre code de tracé ici des erreurs ici"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dbe4bc43",
"metadata": {},
"outputs": [],
"source": [
"import Base: +, -, *, /, inv, conj, convert, promote_rule\n",
"using LinearAlgebra\n",
"\n",
"struct D <: Number\n",
" f::Tuple{Float64,Float64}\n",
"end\n",
"+(x::D, y::D) = D(x.f .+ y.f)\n",
"-(x::D, y::D) = # à compléter\n",
"*(x::D, y::D) = D((x.f[1]*y.f[1], (x.f[2]*y.f[1] + x.f[1]*y.f[2])))\n",
"/(x::D, y::D) = # à compléter\n",
"convert(::Type{D}, x::Real) = D((x,zero(x)))\n",
"promote_rule(::Type{D}, ::Type{<:Real}) = D\n",
"Base.show(io::IO,x::D) = print(io,x.f[1],x.f[2]<0 ? \" - \" : \" + \",abs(x.f[2]),\" ε\")\n",
"ε = D((0,1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "437e8442",
"metadata": {},
"outputs": [],
"source": [
"@show (1+2ɛ)*(3+4ɛ) ; @assert (1+2ɛ)*(3+4ɛ) == 3+10ɛ \"erreur\"\n",
"# Compléter avec d'autres vérifications\n",
";"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15733be0",
"metadata": {},
"outputs": [],
"source": [
"for a in (0.1, 2, 25) \n",
" # Votre code de vérification de l'algorithme avec nombres duaux ici\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "24c36846",
"metadata": {},
"outputs": [],
"source": [
"# Votre code avec les tracés ici"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69dfb790",
"metadata": {},
"outputs": [],
"source": [
"function nthrt(a, p=2; x=1, n=10)\n",
" # Votre code sur l'extension à la racine pième ici \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fc33ab1f",
"metadata": {},
"outputs": [],
"source": [
"a = 2 ; p = 3\n",
"# Votre code pour comparer le résultat de `nthrt` à la solution analytique"
]
},
{
"cell_type": "markdown",
"id": "06b9ad64",
"metadata": {},
"source": [
"--------------------\n",
"--------------------\n",
"--------------------"
]
},
{
"cell_type": "markdown",
"id": "b0e42447",
"metadata": {},
"source": [
"### Extension de la différentiation automatique au second ordre et à plusieurs variables : application à une méthode des moindres carrés non linéaire\n",
"\n",
"L'objectif de cet exercice est d'étendre le concept de différentiation automatique au cas d'une fonction scalaire de plusieurs variables au second ordre, autrement dit de permettre le calcul du gradient et de la hessienne en $a$ de la fonction deux fois différentiable $f: Ω⊂\\mathbb{R}^N → \\mathbb{R}$.\n",
"\n",
"En partant de son développement de Taylor à l'ordre 2\n",
"\n",
"\n",
"$$\n",
"\\tag{3}\n",
"x₀, ɛ ∈ \\mathbb{R}^N,\\quad\n",
"f(x₀+ɛ)=f(x₀)+∇f(x₀)ᵀɛ+\\frac{1}{2}ɛᵀ∇²f(x₀)ɛ+ℴ(ɛ²)\n",
"\\quad\\textrm{avec}\\quad\n",
"f(x₀)∈\\mathbb{R},\\,∇f(x₀)∈\\mathbb{R}^N,\\,∇²f(x₀)∈\\mathbb{R}^{N×N}\n",
"$$\n",
"\n",
"on a l'idée d'introduire une nouvelle classe de nombres scalaires formés par un triplet constitué d'un scalaire $a$, d'un vecteur $b$ et d'une matrice carrée symétrique $c$ et d'écrire\n",
"\n",
"\n",
"$$\n",
"\\tag{4}\n",
"x=a+bᵀɛ+\\frac{1}{2}\\mathrm{Tr}(cη)\n",
"$$\n",
"\n",
"où $ɛ$ et $η$ sont ici respectivement un vecteur de $\\mathbb{R}^N$ et une matrice de $\\mathbb{R}^{N×N}$ telles que $ɛɛᵀ=η$, $ηɛ=0$ et $η²=0$. On peut ainsi faire l'économie de $η$ dans (4) et réécrire tout nombre $x$ sous la forme\n",
"\n",
"\n",
"$$\n",
"\\tag{5}\n",
"x=a+bᵀɛ+\\frac{1}{2}ɛᵀcɛ\n",
"$$\n",
"\n",
"Une telle famille de nombres est représentée ci-dessous par un `struct DD{N} <: Number` (où `N` permet de paramétriser ce nouveau type par la dimension $N$ et `<: Number` indique que le type dérive du type abstrait `Number` donnant un sens aux opérations usuelles) contenant les données membres `val`, `grad` et `hess` désignant respectivement $a$, $b$ et $c$ dans la décomposition (5). A noter que 3 constructeurs sont définis pour ce type\n",
"\n",
"- `DD(val, grad, hess)` définit un nombre `DD` à partir de ses 3 données membres `val`, `grad` et `hess`. Il n'est pas utile de préciser `N` car la dimension est déduite de `length(grad)`.\n",
"\n",
"- `DD(val, grad)` définit un nombre `DD` à partir de `val` et `grad` et fixe `hess` à la matrice nulle de dimension cohérente avec celle de `grad`.\n",
"\n",
"- `DD{N}(val)` définit un nombre `DD` à partir de `val` uniquement. Il est alors nécessaire de préciser dans la construction la valeur de `N` qui ne peut se déduire de `val`. `grad` et `hess` sont alors initialisés à des valeurs nulles de dimension cohérente avec `N`.\n",
"\n",
"\n",
"Il est naturel de surcharger les opérateurs de base de la façon suivante\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"\\left(a+bᵀɛ+\\frac{1}{2}ɛᵀcɛ\\right)+\\left(a'+b'ᵀɛ+\\frac{1}{2}ɛᵀc'ɛ\\right)&=(a+a')+(b+b')ᵀɛ+\\frac{1}{2}ɛᵀ(c+c')ɛ\\\\\n",
"\\left(a+bᵀɛ+\\frac{1}{2}ɛᵀcɛ\\right)-\\left(a'+b'ᵀɛ+\\frac{1}{2}ɛᵀc'ɛ\\right)&=(a-a')+(b-b')ᵀɛ+\\frac{1}{2}ɛᵀ(c-c')ɛ\\\\\n",
"\\left(a+bᵀɛ+\\frac{1}{2}ɛᵀcɛ\\right)*\\left(a'+b'ᵀɛ+\\frac{1}{2}ɛᵀc'ɛ\\right)&=(aa')+(ab'+a'b)ᵀɛ+\\frac{1}{2}ɛᵀ(ac'+a'c+bb'ᵀ+b'bᵀ)ɛ\n",
"\\end{align*}\n",
"$$\n",
"\n",
"1. Compléter la surcharge des opérateurs `-` et `*` dans le code. \n",
" Remarque : la syntaxe `where {N}` rappelle que `N` est un paramètre qui sera remplacé à la compilation par la valeur voulue. \n",
"\n",
"1. On choisit de construire l'opérateur de division à partir de la définition préalable de la function d'inversion `inv`. Établir à la main l'expression des données membres `val` et `grad` et fixe `hess` de l'inverse d'un nombre de type `DD{N}` par identification en écrivant que la mutiplication d'un nombre par son inverse correspond à `val=1`, `grad=0` et `hess=0`. Compléter le code de la fonction `inv(x::DD{N}) where {N}`.\n",
"\n",
"1. On se place en dimension $N=2$. On peut donc introduire les deux nombres\n",
"\n",
" ```julia\n",
" εˣ = DD(0,[1,0])\n",
" εʸ = DD(0,[0,1])\n",
" ```\n",
"\n",
" Tester les développements des opérateurs sur des cas simples `1+εˣ`, `1+εˣ+εʸ` et sur des cas plus complexes `(1+εˣ)*(2+3εʸ)` et `εˣ/(1+εˣ+εʸ)` en comparant avec les valeurs calculées à la main. Par exemple\n",
"\n",
" $$\n",
" ε^x ε^y=\\frac{1}{2} \\begin{pmatrix}ε^x \\\\ ε^y\\end{pmatrix}^T \\begin{pmatrix}0 & 1 \\\\ 1 & 0\\end{pmatrix}\\begin{pmatrix}ε^x \\\\ ε^y\\end{pmatrix}\n",
" ⟹\n",
" \\left(0, \\begin{pmatrix}1 \\\\ 0\\end{pmatrix}, \\begin{pmatrix}0 & 0 \\\\ 0 & 0\\end{pmatrix}\\right)\n",
" *\n",
" \\left(0, \\begin{pmatrix}0 \\\\ 1\\end{pmatrix}, \\begin{pmatrix}0 & 0 \\\\ 0 & 0\\end{pmatrix}\\right)\n",
" =\n",
" \\left(0, \\begin{pmatrix}0 \\\\ 0\\end{pmatrix}, \\begin{pmatrix}0 & 1 \\\\ 1 & 0\\end{pmatrix}\\right)\n",
" $$\n",
"\n",
"1. Obtenir le gradient et la hessienne de $f(x,y)=\\frac{(x+y)^2}{2}$ par différentiation automatique au point $(1,1)$.\n",
"\n",
"1. Soit un entier arbitraire $n$, $a$, $b$ et $c$ respectivement un scalaire de $[0,1]$, un vecteur de $[0,1]^n$ et une matrice symétrique de $[0,1]^{n×n}$ aléatoires. Implémenter la fonction $f(x)=a+bᵀx+\\frac{1}{2}xᵀcx$ pour $x∈\\mathbb{R}^n$. \n",
"\n",
" En exploitant judicieusement le vecteur de nombres\n",
"\n",
" ```julia\n",
" ε = [DD(0,(1:n .== i)) for i in 1:n]\n",
" ```\n",
"\n",
" et un tirage aléatoire d'un vecteur $x$, vérifier par différentiation automatique que $∇f(x)=b+cx$ et $∇²f(x)=c$.\n",
"\n",
"1. On souhaite maintenant réaliser un programme permettant de minimiser une fonction $J$ arbitraire dépendant de plusieurs variables\n",
"\n",
" $$\n",
" \\min_{(p_1,\\ldots,p_n)\\in\\mathbb{R}^n} J(p_1,\\ldots,p_n)\n",
" $$\n",
"\n",
" On suppose que le problème revient à chercher un vecteur de paramètres annulant son gradient par la méthode de Newton-Raphson. Autrement dit à chaque étape il est nécessaire d'évaluer le gradient ainsi que la hessienne de $J$\n",
" \n",
" $$\n",
" p^{k+1}=p^{k}-∇²J(p^{k})^{-1}∇J(p^{k})\n",
" $$\n",
" \n",
" Écrire une fonction générique `minimizer` prenant comme entrées un vecteur de paramètres initiaux `p⁰`, une fonction `J`, un nombre maximal d'itérations `maxiter` (par défaut `100`) et un paramètre d'arrêt `ϵ` (par défaut `1.e-15`). On considérera un algorithme de Newton-Raphson avec un critère d'arrêt $\\lVert ∇J(p^{k}\\rVert<ϵ$.\n",
"\n",
"1. Application n°1.\n",
"\n",
" On dispose de $n$ points $(x_i, y_i)$ d'une fonction inconnue $y = f(x)$.\n",
" \n",
" ```julia\n",
" x = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1.0]\n",
" y = [0.6761488864859304; 0.6345697680852508; 0.6396283580587062; 0.6132010027973919;\n",
" 0.5906142598705267; 0.5718728461471725; 0.5524549902830562; 0.538938885654085;\n",
" 0.5373495476994958; 0.514904589752926; 0.49243437874655027]\n",
" ```\n",
" \n",
" On souhaite approximer $f$ par une fonction de la forme\n",
"\n",
" $$\n",
" \\widetilde f(x) = \\frac{a}{b + x}\n",
" $$\n",
" en minimisant\n",
" $$\n",
" \\sum_{i=1}^{n} |\\widetilde f(x_i) - y_i|^2\n",
" $$\n",
"\n",
" Écrire un code formant la fonction à minimiser et appliquer `minimizer` pour obtenir la meilleure approximation $\\widetilde f$. Tracer sur le même graphe les points donnés et la fonction approximante.\n",
"\n",
"1. Application n°2.\n",
"\n",
" On se donne $n$ nouveaux points $(x_i, y_i)$ d'une fonction inconnue $y = f(x)$\n",
" \n",
" ```julia\n",
" x = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1.0]\n",
" y = [-0.9187980789440975; -0.6159791344678258; -0.25568734869121856;\n",
" -0.14269370171581808; 0.3094396057228459; 0.6318327173549161;\n",
" 0.8370437988106428; 1.0970402798788812; 1.6057799131867696;\n",
" 1.869090784869698; 2.075369730726694]\n",
" ``` \n",
" \n",
" et on souhaite approximer $f$ par une fonction affine\n",
"\n",
" $$\n",
" \\widetilde f(x) = ax+b\n",
" $$\n",
"\n",
" en minimisant la somme des distances euclidiennes entre les points et la droite définie par $\\widetilde f$. Etant donné que la distance entre un point $(x_i,y_i)$ et la ligne droite est donnée par\n",
"\n",
" $$\n",
" \\frac{\\lvert y_i - a x_i - b \\rvert}{\\sqrt{1+a^2}}\n",
" $$\n",
"\n",
" la fonction objectif à minimiser s'écrit\n",
"\n",
" $$\n",
" J(a, b) := \\sum_{i=1}^{n} \\frac{ \\left( y_i - a x_i - b \\right)^2 }{1+a^2}\n",
" $$\n",
"\n",
" Trouver les paramètres optimaux $a$ et $b$ à l'aide de `minimizer` et tracer la droite ainsi que les points."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "995fb215",
"metadata": {},
"outputs": [],
"source": [
"import Base: +, -, *, /, inv, conj, convert, promote_rule\n",
"using LinearAlgebra\n",
"\n",
"struct DD{N} <: Number\n",
" val::Float64\n",
" grad::Vector{Float64}\n",
" hess::Symmetric{Float64, Matrix{Float64}}\n",
" DD(val::Real, grad, hess) = new{length(grad)}(val, grad, Symmetric(hess))\n",
" DD(val::Real, grad) = new{length(grad)}(val, grad, Symmetric(zeros(length(grad),length(grad))))\n",
" DD{N}(val::Real) where {N} = new{N}(val, zeros(N), Symmetric(zeros(N,N)))\n",
"end\n",
"val(x::DD{N}) where {N} = x.val\n",
"grad(x::DD{N}) where {N} = x.grad\n",
"hess(x::DD{N}) where {N} = x.hess\n",
"conj(x::DD{N}) where {N} = DD(conj(x.val),conj(x.grad),conj(x.hess))\n",
"+(x::DD{N}, y::DD{N}) where {N} = DD(x.val+y.val,x.grad+y.grad,x.hess+y.hess)\n",
"-(x::DD{N}, y::DD{N}) where {N} = # à compléter ici\n",
"*(x::DD{N}, y::DD{N}) where {N} = # à compléter ici\n",
"inv(x::DD{N}) where {N} = # à compléter ici\n",
"/(x::DD{N}, y::DD{N}) where {N} = x*inv(y)\n",
"convert(::Type{DD{N}}, x::Real) where {N} = DD{N}(x)\n",
"promote_rule(::Type{DD{N}}, ::Type{<:Real}) where {N} = DD{N}\n",
"Base.show(io::IO,x::DD{N}) where {N} = print(io,x.val,\" + \",x.grad,\"ᵀ𝛆\",\" + ½ 𝛆ᵀ\",x.hess,\"𝛆\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aea3d848",
"metadata": {},
"outputs": [],
"source": [
"εˣ = DD(0,[1,0])\n",
"εʸ = DD(0,[0,1])\n",
"@show 1+εˣ\n",
"# Votre code de vérification sur plusieurs expressions ici\n",
";"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7f3ad593",
"metadata": {},
"outputs": [],
"source": [
"f(x,y) = (x+y)^2/2\n",
"@show # Votre code de calcul du gradient et de la hessienne de f en (1,1) ici"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "be44974c",
"metadata": {},
"outputs": [],
"source": [
"function test_diff(n)\n",
" # Votre code ici\n",
"end\n",
"\n",
"test_diff(1)\n",
"test_diff(10)\n",
"test_diff(100)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b81d8ee9",
"metadata": {},
"outputs": [],
"source": [
"function minimizer(p⁰, J::Function, ϵ = 1.e-15, maxiter = 100)\n",
" n = length(p⁰)\n",
" p = p⁰ + [DD(0,(1:n .== i)) for i in 1:n]\n",
" # Votre code ici\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4fd6ffc8",
"metadata": {},
"outputs": [],
"source": [
"x = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1.0]\n",
"y = [0.6761488864859304; 0.6345697680852508; 0.6396283580587062; 0.6132010027973919;\n",
" 0.5906142598705267; 0.5718728461471725; 0.5524549902830562; 0.538938885654085;\n",
" 0.5373495476994958; 0.514904589752926; 0.49243437874655027]\n",
"f(a,b) = x -> a / (b+x)\n",
"# Votre code ici"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9c580d22",
"metadata": {},
"outputs": [],
"source": [
"x = [0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.6; 0.7; 0.8; 0.9; 1.0]\n",
"y = [-0.9187980789440975; -0.6159791344678258; -0.25568734869121856;\n",
" -0.14269370171581808; 0.3094396057228459; 0.6318327173549161;\n",
" 0.8370437988106428; 1.0970402798788812; 1.6057799131867696;\n",
" 1.869090784869698; 2.075369730726694]\n",
"f(a,b) = x -> a*x+b\n",
"# Votre code ici"
]
}
],
"metadata": {
"citation-manager": {
"items": {}
},
"kernelspec": {
"display_name": "Julia 1.9.0",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}