Cours ENPC - Pratique du calcul scientifique

Examen final

• Ce notebook est à soumettre sur Educnet avant 16h (ou 16h40 pour les quelques étudiants bénéficiant de temps supplémentaire).

• L’examen comporte trois exercices indépendants. Dans chaque exercice les cellules peuvent éventuellement dependre des cellules précèdentes.

• Afin de faciliter l'évaluation de votre code, ne pas changer les signatures des fonctions à implémenter.

• La cellulle ci-dessous importe les packages utilisés dans ce notebook et définit une macro qui a pour but de faciliter les tests unitaires figurant dans le sujet. Il est ainsi nécessaire d'exécuter le code dans cette cellule avant de poursuivre le notebook.

using Polynomials
using Plots
using LaTeXStrings

macro mark(bool_expr)
    return :( print($bool_expr ? "✔️" : "❌"))
end

@mark (macro with 1 method)

1. Exercice sur l'équation de Laplace

On se propose d'implémenter une méthode numérique pour résoudre approximativement l'équation de Laplace avec conditions au bord homogène de Dirichlet:

$$ u\in C^2([0,1]),\quad\left\{\begin{aligned} -u''(x) & = \varphi(x) & \forall\, x\in(0,1),\\ u(0) & = u(1) = 0. \end{aligned}\right.$$ Pour cela, on approxime $\varphi$ avec un polynome interpolateur $\widehat \varphi$, puis on résout exactement l'équation de Laplace avec membre de droite $\widehat \varphi$ au lieu de $\varphi$.

1. Écrire une fonction approx(n) implémentant cette approche. La fonction devra renvoyer une approximation polynomiale de la solution basée sur une interpolation de degré $n$ du membre de droite à des points équidistants entre 0 et 1 compris. On prendra comme membre de droite la fonction $$\varphi(x) = \exp\Bigl(\sin(2\pi x)\Bigr) \Bigl(\sin(2\pi x)-\cos(2\pi x)^2 \Bigr),$$ auquel cas la solution analytique est donnée par $$ u(x)=(2\pi)^{-2}\Bigl(\exp\bigl(\sin(2\pi x)\bigr)-1\Bigr).$$

Indications :

• Utiliser la fonction fit de la bibliothèque Polynomials.jl pour obtenir le polynôme interpolateur:

   p = fit(x,y)
  

  où x sont les noeuds d'interpolation et y sont les valeurs de la fonction à interpoler.

• Pour calculer la solution analytique avec membre de droite polynomial, on pourra remarquer que toutes les solutions sont des polynômes, et que, sans condition au bord, la solution est unique modulo $\mathbf{P}_1$.

• On pourra utiliser la fonction integrate de la bibliothèque Polynomials.jl qui calcule une primitive d'un polynôme:

  P = integrate(p)
  

• Utiliser le format BigFloat pour limiter les erreurs d'arrondi. En particulier, la fonction LinRange{BigFloat}(a, b, N) permet de créer un vecteur de N nombres au format BigFloat également distribués entre a et b.

φ(x) = exp(sin(2π*x))*(sin(2π*x)-cos(2π*x)^2)
u(x) = (exp(sin(2π*x))-1)/(2π)^2

function approx(n)
    X = LinRange{BigFloat}(0, 1, n + 1)
    Y = φ.(X)
    p = fit(X, Y)
    uh = - integrate(integrate(p))
    uh + Polynomial([-uh(0), uh(0) - uh(1)])
end

@mark typeof(approx(3)) == Polynomial{BigFloat, :x}
@mark degree(approx(3)) == 5
@mark approx(3)(0) ≈ 0
@mark approx(3)(1) ≈ 0

✔️✔️

✔️✔️

2. Faire une animation permettant de visualiser la solution exacte et la solution approchée pour des valeurs de n allant de 2 à 20, en utilisant la fonction animate de la bibliothèque Plots.

xs = LinRange(0, 1, 200)
anim = @animate for n in 2:20
    un = approx(n)
    plot(xs, u.(xs), label="Exact solution")
    plot!(xs, un.(xs), label="Approximate solution n=$n")
end
gif(anim, "interpolation.gif", fps=1, show_msg=false)

3. Écrire une fonction estimate_error(n) qui approxime l'erreur, en norme $L^\infty$, entre la solution approchée par l'approche ci-dessus et la solution exacte.

function estimate_error(n)
    un = approx(n)
    x_fine = LinRange(0, 1, 1000)
    un_fine, u_fine = un.(x_fine), u.(x_fine)
    return maximum(abs.(u_fine - un_fine))
end

@mark estimate_error(2) > 0
@mark estimate_error(20) < 1e-3
@mark estimate_error(20) > 1e-8
@mark estimate_error(40) < 1e-6

✔️✔️

✔️✔️

4. Tracer un graphique de l'erreur en fonction de $n$ pour $n$ allant de 5 à 50. Utiliser l'échelle par défaut pour l'axe des abcisses et une échelle logarithmique pour l'axe des ordonnées.

ns = 5:50
errors = estimate_error.(ns)
plot(ns, errors, marker = :circle, label=L"$L^{\infty}$ Error")
plot!(yaxis=:log, lw=2)

2. Calcul des racines cubiques de 1

On souhaite calculer numériquement les solutions dans $\mathbb C$ de l'équation $$ f(z) = 0, \qquad f(z) := z^3 - 1. $$ Pour ce faire, on utilisera l'itération de Newton-Raphson $$ z_{k+1} = z_k - \frac{f(z_k)}{f'(z_k)}, \tag{NR} $$ où $f'\colon \mathbb C \to \mathbb C$ est la dérivée complexe de la fonction $f$, ici donnée par $f'(z) = 3z^2$.

1. Écrire une fonction newton_raphson(f, df, z₀; maxiter = 100, ε = 1e-12) renvoyant le résultat de l'itération (NR) appliquée à la fonction f et initialisée à z₀, ou nothing si une solution n'a pas été trouvée après maxiter itérations. L'argument df est la dérivée complexe de la fonction f, et on utilisera comme critère d'arrêt $$ |f(z_k)| ≤ \varepsilon. $$

function newton_raphson(f, df, z₀; maxiter = 100, ε = 1e-12)
    for i in 1:maxiter
        z₀ -= f(z₀) / df(z₀)
        (abs∘f)(z₀) ≤ ε && return z₀
    end
    return nothing
end

@mark newton_raphson(z -> z^2 - 2, z -> 2z, 1) ≈ √2
@mark newton_raphson(z -> z^2 - 2, z -> 2z, -1) ≈ -√2
@mark newton_raphson(z -> z^3 - 2, z -> 3z^2, 1) ≈ cbrt(2)
@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1 + im) ≈ im
@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1 - im) ≈ -im
@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1.5) == nothing
@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1) == nothing

✔️✔️✔️✔️

✔️✔️✔️

2. On appelle le bassin d'attraction d'une racine l'ensemble des points de départ $z₀$ tels que l'itération de Newton-Raphson converge vers cette racine. En vue de calculer numériquement les bassins d'attraction des trois racines, écrire une fonction which_root(z₀) qui, pour un argument z₀,

   • renvoit 0 si la méthode initialisée en z₀ ne converge pas;

   • renvoit 1 si la méthode initialisée en z₀ converge vers $1$;

   • renvoit 2 si la méthode initialisée en z₀ converge vers $\exp(2i\pi/3)$;

   • renvoit 3 si la méthode initialisée en z₀ converge vers $\exp(4i\pi/3)$.

f(z) = z^3 - 1
df(z) = 3z^2

function which_root(z₀)
    roots = exp.([0, im*2π/3, im*4π/3])
    r = newton_raphson(f, df, z₀)
    r == nothing && return 0
    return argmin(abs.(roots .- r))
end

@mark which_root(1) == 1
@mark which_root(im) == 2
@mark which_root(-im) == 3

✔️✔️

✔️

3. Calculer numériquement et représenter graphiquement les bassins d'attraction des trois racines dans le domaine $[-2, 2] \times [-2, 2]$ du plan complexe. La fonction heatmap avec argument levels=4 pourra étre utilisée pour la représentation graphique. Cette fonction s'utilise de la même manière que les fonctions contour et contourf.

# Exact roots
X = LinRange(-2, 2, 2000)
Y = LinRange(-2, 2, 2000)
Z = [x + y*im for y in Y, x in X]
R = which_root.(Z)

using Plots
heatmap(X, Y, R, c=:viridis, levels=4, fg=:white, bg=:transparent)

4. On considère à présent une autre approche, appelée méthode de la sécante, pour le calcul des racines cubiques complexes de 1. Une itération de cette méthode est donnée par $$ z_{k+2} = \frac{z_{k} f(z_{k+1}) - z_{k+1}f(z_k)}{f(z_{k+1}) - f(z_{k})}. $$ Écrire une fonction secant(f, z₀, z₁; maxiter = 100, ε = 1e-12) qui, contrairement à la méthode newton_raphson ci-dessus, devra renvoyer toutes les itérations obtenues lorsque la méthode de la sécante est appliquée à la fonction f et initialisée à z₀ et z₁, ou nothing si une solution n'a pas été trouvée après maxiter itérations. On utilisera comme critère d'arrêt $$ |f(z_k)| ≤ \varepsilon. $$

function secant(f, z₀, z₁; maxiter = 100, ε = 1e-12)
    zs = [z₀, z₁]
    for i in 1:maxiter
        z₀, z₁ = zs[end-1:end]
        push!(zs, (z₀*f(z₁) - z₁*f(z₀)) / (f(z₁) - f(z₀)))
        (abs∘f)(zs[end]) ≤ ε && return zs
    end
    return nothing
end

@mark secant(z -> z^2 - 2, 1., 2.)[1:2] ≈ [1. ; 2.]
@mark secant(z -> z^2 - 2, 1., 2.)[end] ≈ √2
@mark secant(z -> z^2 - 2, -1., -2.)[end] ≈ -√2
@mark secant(z -> z^2 + 1, 1. + im, 1. + 2im)[end] ≈ im
@mark secant(z -> z^2 + 1, 1. + im, 1. - 2im)[end] ≈ -im
@mark secant(z -> z^2 + 1, 1.,  2.) == nothing

✔️✔️

✔️✔️

✔️✔️

5. Il est possible de démontrer que si cette itération converge vers une solution $r$, alors $$ \tag{Q} \lim_{k \to \infty} \frac{|z_{k+1} - r|}{|z_k - r|^q} \in \mathbb R_{>0}, $$ pour un réel positif $q$ qui est l'ordre de convergence. Calculer empiriquement dans la fonction estimate_q la valeur de $q$ en appliquant à la fonction $z \mapsto z^3 - 1$ la méthode de la sécante initialisée avec $$ z_0 = -1, \qquad z_1 = - i. $$ Il sera utile d'utiliser le format Complex{BigFloat} et d'appeler la fonction secant avec un petit ε afin d'éviter que l'erreur n'atteigne trop vite l'epsilon machine. C'est ce qui est fait dans le début de code qui vous est donné ci-dessous. La valeur correcte de $q$ est donnée par $$ q = \frac{1 + \sqrt{5}}{2} \approx 1.618. $$

   Remarque. Soit $y_k := - \log(|z_{k} - r|)$. L'équation (Q) implique que $$ \lim_{k \to \infty} y_{k+1} - q y_{k} = C \in \mathbb R. $$ On en déduit que $$ q = \lim_{k\to \infty} \frac{y_{k+1}}{y_k} $$ Cette équation permet d'estimer $q$ à partir de $y_{k+1}$ et $y_k$ pour $k$ suffisamment grand.

# Set number of significant digits in base 10
precision = 1_000
setprecision(precision, base=10)

# Initial values for the secant iteration
z₀ = big(-1.)
z₁ = big(-1.0im)

# Exact root towards which convergence should occur
root = exp(im * 4big(π)/3)

function estimate_q()
    zs = secant(f, z₀, z₁, ε = 1e-50)
    ys = @. log(abs(zs - root))
    return (ys[2:end] ./ ys[1:end-1])[end]
end

@mark abs(estimate_q() - (1 + √5)/2) < 1e-2

✔️

3. Résolution d'une équation différentielle

On s'intéresse dans cet exercice au calcul de la trajectoire d'une balle de golf frappée par un joueur. Pour simplifier, on supposera que la trajectoire de la balle est confinée dans un plan perpendiculaire à la surface de jeu, qu'on supposera parfaitement horizontale. La position de la balle à un instant donné peut alors être décrite par les coordonnées horizontale ($x_1$) et verticale ($x_2$) dans le plan de la trajectoire. Le mouvement peut être régit par l'équation de Newton: $$ m \mathbf x''(t) = - \mathbf F_x(\lVert \mathbf v \rVert) \frac{\mathbf v}{\lVert \mathbf v \rVert} - mg \mathbf e_2, \qquad \mathbf v := \mathbf x'. \tag{Golf} $$ Ici, $m$ est la masse de la balle de golf, $g$ est l'accélération de gravité, $\mathbf e_2$ est le vecteur unitaire $(0, 1)^T$, et $\mathbf F_x(\lVert v \rVert)$ est la force de traînée, qui par analyse dimensionnelle peut s'écrire $$ \mathbf F^d(\lVert \mathbf v \rVert) = \frac{1}{2} \rho A C^d \lVert \mathbf v \rVert^2. $$ La signification des constantes apparaissant dans cette formule, et les valeurs que nous prendrons dans cette exercice, sont données ci-dessous. Noter qu'on suppose, pour simplifier, que le coefficient de traînée $C^d$ est une constante indépendante du nombre de Reynolds de l'écoulement. On néglige en outre les effets de la rotation propre de la balle à l'origine de l'effet Magnus observé en pratique.

# Mass of the golf ball [kg]
const m = .045

# Radius of the golf ball [m]
const r = .021

# Air density at 0 ⁰C [kg/m³]
const ρ = 1.293

# Gravity acceleration [m/s²]
const g = 9.81

# Cross-sectional areal [m²]
const A = π*r^2

# Drag coefficient [dimensionless]
const Cᵈ = .2

# Drag force depending on V := ‖𝐯‖
Fᵈ(V) = 1/2 * ρ*A*Cᵈ*V^2

Fᵈ (generic function with 1 method)

On modélise le joueur par un point au sol, qu'on prend comme origine du repère orthonormé. Pour modéliser le fait que la balle est frappée par le joueur au temps $t = 0$, on prend la condition initiale suivante pour (Golf). $$ \mathbf x(0) = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \qquad \mathbf x'(0) = V_0 \begin{pmatrix} \cos(\theta) \\ \sin(\theta) \end{pmatrix}. \tag{IC} $$ Ici $V_0$ est la vitesse initiale de la balle, et $\theta$ est l'angle de loft.

1. Soit $\mathbf v = \mathbf x'$. L'équation (Golf) munie de la condition initiale (IC) peut être réécrite comme un problème aux valeurs initiales du premier ordre pour le vecteur $\mathbf Z := (x_1, x_2, v_1, v_2)^T$ de la forme suivante: $$ \mathbf Z'(t) = \mathbf f\bigl(t, \mathbf Z(t)\bigr), \qquad \mathbf Z(0) = \mathbf Z_0. \tag{1st-order} $$ (Dans le cas qui nous occupe, la fonction $f$ est en fait indépendante du temps.) Écrire la fonction $f$ sous forme d'une fonction Julia f(t, Z).

function f(t, Z)
   x₁, x₂, v₁, v₂ = Z
   u = √(v₁^2 + v₂^2)
   return [v₁, v₂, - 1/m * Fᵈ(u)/u*v₁, - 1/m*  Fᵈ(u)/u*v₂ - g]
end

@mark f(1, [0; 0; 100; 0]) ≈ [100.0, 0.0, -39.8083771506977, -9.81]
@mark f(0, [100; 100; 100; 100]) ≈ [100.0, 100.0, -56.297546862579914, -66.10754686257991]

✔️✔️

2. Écrire une fonction rk4(tₙ, Zₙ, f, Δ) implémentant un pas de temps de taille $\Delta$ de la méthode de Runge-Kutta d'ordre 4 pour une équation différentielle générique de la forme $Z' = f(t, Z)$. Cette méthode est basée sur l'itération suivante: $$ \mathbf Z_{n+1} = \mathbf Z_n + \frac{\Delta}{6}\left(\mathbf k_1 + 2\mathbf k_2 + 2\mathbf k_3 + \mathbf k_4 \right), $$ où \begin{align*} \mathbf k_1 &= \ f(t_n, \mathbf Z_n), \\ \mathbf k_2 &= \ f\!\left(t_n + \frac{\Delta}{2}, \mathbf Z_n + \frac{\Delta}{2} \mathbf k_1\right), \\ \mathbf k_3 &= \ f\!\left(t_n + \frac{\Delta}{2}, \mathbf Z_n + \frac{\Delta}{2} \mathbf k_2\right), \\ \mathbf k_4 &= \ f\!\left(t_n + \Delta, \mathbf Z_n + \Delta\mathbf k_3\right). \end{align*} La fonction devra renvoyer $\mathbf Z_{n+1}$.

function rk4(tₙ, Zₙ, f, Δ)
    k₁ = f(tₙ,       Zₙ           )
    k₂ = f(tₙ + Δ/2, Zₙ + Δ/2 * k₁)
    k₃ = f(tₙ + Δ/2, Zₙ + Δ/2 * k₂)
    k₄ = f(tₙ + Δ  , Zₙ + Δ   * k₃)
    return Zₙ + Δ/6 * (k₁ + 2k₂ + 2k₃ + k₄)
end

@mark rk4(0., [0.], (t, Z) -> [1.], 1.) ≈ [1.0]
@mark rk4(1., [0.], (t, Z) -> [t], 1.)  ≈ [3/2]
@mark rk4(0., [0.; 0.], (t, Z) -> [t^2; t^3], 1.) ≈ [1/3; 1/4]

✔️✔️

✔️

3. Écrire une fonction solve_ode(V₀, θ, Δ) pour résoudre (1st-order) pour une vitesse initiale V₀ et un angle initial θ, en utilisant la méthode de Runge-Kutta d'ordre 4 avec pas de temps fixe Δ. Votre fonction devra renvoyer un vecteur de temps ts et un vecteur de vecteurs Zs contenant la solution à ces temps. On calculera la trajectoire jusqu'à ce que la balle soit retombée sur le sol. Plus précisément, on demande d'interrompre l'intégration numérique dès que la valeur de la coordonnée $x_2$ sera devenue négative; il faudra donc que Zs[end-1][2] soit positif et Zs[end][2] soit négatif.

   Indication. Pour construire Zs, il est recommandé d'utiliser la fonction push!.

function solve_ode(V₀, θ, Δ)
    Z₀ = [0.; 0.; V₀*cos(θ); V₀*sin(θ)]
    ts = [0.]
    Zs = [Z₀]
    while Zs[end][2] ≥ 0.
        push!(Zs, rk4(ts[end], Zs[end], f, Δ))
        push!(ts, ts[end] + Δ)
    end
    return ts, Zs
end

@mark solve_ode(50, π/4, .01) |> length == 2
@mark solve_ode(50, π/4, .01)[2][end-1][2] ≥ 0
@mark solve_ode(50, π/4, .01)[1][1:10] ≈ 0:.01:.09
@mark solve_ode(50, π/4, .01)[2][end][2] ≤ 0
@mark solve_ode(50, π/4, .01)[2][end][1] ≈ 149.30535166172655

✔️✔️

✔️✔️✔️

4. Écrire une fonction my_plot(Δ, V₀, θs) permettant d'illustrer sur un même plot la trajectoire de la balle dans le plan $x_1 \, x_2$ pour une valeur de $V_0$ donnée et plusieurs valeurs de l'angle de loft contenues dans le vecteur θs.

function my_plot(V₀, θs, Δ)
    p = plot(title="Trajectories of golf ball")
    for θ in θs
        ts, Zs = solve_ode(V₀, θ, Δ)
        xs = [z[1] for z in Zs]
        ys = [z[2] for z in Zs]
        plot!(xs, ys, label="\\theta = $(round(θ*180/π,digits=1))")
    end
    return p
end
my_plot(50, [π/8, π/4, 3π/8], .01)

5. Écrire une fonction distance(V₀, θ, Δ) pour calculer la distance entre le joueur et le point d'impact de la bal sur le sol. Pour ce faire, résoudre l'équation différentielle avec un pas de temps Δ et trouver le point d'impact par interpolation linéaire sur le dernier pas de temps, durant lequel la hauteur de la balle passe d'une valeur positive à une valeur négative.

   Indication. Il n'est pas nécessaire d'utiliser de bibliothèque externe car il suffit de trouver l'abscisse de l'intersection avec l'axe $y=0$ de la droite passant par les deux derniers points de la trajectoire $(x_1, y_1)$ et $(x_2,y_2)$ (on rappelle que par construction $y_1≥0$ et $y_2<0$).

function distance(V₀, θ, Δ)
    ts, Zs = solve_ode(V₀, θ, Δ)
    x₁, y₁ = Zs[end-1][1:2]
    x₂, y₂ = Zs[end][1:2]
    return x₁ + y₁ * (x₂ - x₁) / (y₁ - y₂)
end

@mark distance(50, π/8, .01) ≈ 124.20307897577761
@mark distance(50, π/4, .01) ≈ 149.29126957486164
@mark distance(50, 3π/8, .01) ≈ 102.42272578334352
@mark distance(50, π/2 - .001, .01) < 1.

✔️✔️✔️✔️

6. Faire un plot de la distance parcourue pour par la balle en fonction de l'angle $\theta \in (0, \pi/2)$ pour $V_0 = 50m/s$, estimée avec $\Delta = .01s$. Estimer graphiquement les angles permettant d'atteindre une distance de 80m.

   Indication : il peut être utile de passer xticks=0.:5.:90. en argument à la fonction plot pour simplifier cette estimation.

V₀, Δ = 50, .01
θs = LinRange(0, π/2, 100)
plot(θs * 180/π, distance.(V₀, θs, Δ), label="Distance")
plot!(θs * 180/π, zero(θs) .+ 80, color=:red, label="80m", xticks=0.:5.:90.)
plot!(xlabel="\\theta (en degrés)", xlims=(0, 90))

7. Pour une vitesse initiale $V_0$ donnée, on veut calculer l'angle de loft $\theta$ permettant d'atteindre un trou situé à une distance $L$ du joueur. Pour ce faire, on se propose de mettre en œuvre l'algorithme de Newton-Raphson sur une fonction scalaire prenant comme argument $\theta$ et s'annulant lorsque l'estimation de la distance obtenue par la fonction distance est égale à $L$. La méthode de Newton-Raphson nécessite de connaître la dérivée de la fonction dont on cherche une racine, ce qui peut être fait par différentiation automatique. On donne ci-dessous la structure D de nombre dual avec presque toutes les fonctions suffisantes pour entreprendre la résolution de l'équation différentielle, sauf les fonctions cos(x::D), sin(x::D) et sqrt(x::D) qui sont à définir.

import Base: +, -, *, /, ==, ≈, sqrt, cos, sin, inv, conj, abs, isless, AbstractFloat, convert, promote_rule, show
struct D <: Number
    f::Tuple{Float64,Float64}
end
+(x::D, y::D)                           = D(x.f .+ y.f)
-(x::D, y::D)                           = D(x.f .- y.f)
-(x::D)                                 = D(.-(x.f))
*(x::D, y::D)                           = D((x.f[1]*y.f[1], (x.f[2]*y.f[1] + x.f[1]*y.f[2])))
/(x::D, y::D)                           = D((x.f[1]/y.f[1], (y.f[1]*x.f[2] - x.f[1]*y.f[2])/y.f[1]^2))
==(x::D, y::D)                          = x.f[1] == y.f[1] && x.f[2] == y.f[2]
≈(x::D, y::D)                           = x.f[1] ≈ y.f[1] && x.f[2] ≈ y.f[2]
sqrt(x::D)                              = D((sqrt(x.f[1]), 1/(2sqrt(x.f[1])) * x.f[2]))
cos(x::D)                               = D((cos(x.f[1]), -sin(x.f[1]) * x.f[2]))
sin(x::D)                               = D((sin(x.f[1]), cos(x.f[1]) * x.f[2]))
abs(x::D)                               = D((abs(x.f[1]), x.f[2]*sign(x.f[1])))
abs2(x::D)                              = D((abs2(x.f[1]), 2x.f[1]*x.f[2]))
isless(x::D, y::D)                      = isless(x.f[1], y.f[1])
isless(x::D, y::Real)                = isless(x.f[1], y)
isless(x::Real, y::D)                = isless(x, y.f[1])
D(x::D)                                 = x
AbstractFloat(x::D)                     = x.f[1]
convert(::Type{D}, x::Real)             = D((x,zero(x)))
convert(::Type{Real}, x::D)          = x.f[1]
promote_rule(::Type{D}, ::Type{<:Real}) = D
show(io::IO,x::D)                       = print(io,x.f[1],x.f[2]<0 ? " - " : " + ",abs(x.f[2])," ε")
ε = D((0, 1))

@mark sin(π/4 + ε) ≈ √2/2 * (1. + ε)
@mark cos(ε) == 1
@mark sin(ε) ≈ ε
@mark sqrt(1 + ε) ≈ 1 + .5ε

✔️✔️✔️✔️

8. Ecrire une fonction newton_raphson_dual(f, x; maxiter = 100, δ = 1e-12) de résolution par Newton-Raphson d'une équation scalaire f(x) = 0, dans laquelle la dérivée de f au point courant est obtenue par exploitation des nombres duaux, avec un nombre d'itérations maximal maxiter ($100$ par défaut) et une tolérance δ ($10^{-12}$ par défaut) pour un critère d'arrêt $|f(x)| < δ$.

   Indication : à chaque itération de la résolution, les valeurs de f et de sa dérivée en x peuvent être extraites du calcul de y = f(x + D((0,1))).

function newton_raphson_dual(f, x, maxiter=100; δ = 1e-12)
    for i in 1:maxiter
        y = f(x + ε)
        x -= y.f[1]/y.f[2]
        abs(f(x)) < δ && return x
    end
    error("Failed to converge!")
end

@mark newton_raphson_dual(x -> x^2 - 2, 1) ≈ √2
@mark newton_raphson_dual(x -> x^2 - 2, -1) ≈ -√2
@mark newton_raphson_dual(x -> x^3 - 2, 1) ≈ cbrt(2)
@mark newton_raphson_dual(x -> cos(x) - .5, 2) ≈ acos(.5)

✔️✔️

✔️✔️

9. Écrire une fonction find_angle(L, V₀, Δ, θ₀) renvoyant un angle qui permet d'atteindre une distance $L$ pour une vitesse initiale $V_0$. On supposera que les arguments de la fonction sont tels qu'un tel angle existe. Calculer une estimation des deux angles θ₁ et θ₂ permettant d'atteindre une distance de 80m et tracer les trajectoires correspondantes.

function find_angle(L, V₀, Δ, θ₀)
    obj = θ -> distance(V₀, θ, Δ) - L
    newton_raphson_dual(obj, θ₀)
end

V₀, Δ = 50, .01

θ₁ = find_angle(80, V₀, Δ, 25*π/180)
θ₂ = find_angle(80, V₀, Δ, 60*π/180)

my_plot(50, [θ₁, θ₂], .01)

10. Ayant constaté graphiquement à la question 6 l'existence d'un angle maximisant la portée, estimer la valeur de cet angle en utilisant la méthode de dichotomie (aussi appelée méthode de la bissection) initialisée avec $a = \pi/20$ et $b = \pi/2$.

function bisection(f, a, b; δ = 1e-10)
    @assert f(a) * f(b) ≤ 0
    while abs(b - a) ≥ δ
        x = (a + b) / 2
        a, b = f(a) * f(x) ≤ 0 ? [a, x] : [x, b]
    end
    return (a + b) / 2
end

@mark bisection(x -> x^2 - 2, 1, 2) ≈ √2
@mark bisection(x -> x^3 - 2, 1, 2) ≈ cbrt(2)
@mark bisection(x -> cos(x) - .5, 1, 2) ≈ acos(.5)

V₀, Δ = 50, .01
θᵐᵃˣ = bisection(θ -> distance(V₀, θ + ε, Δ).f[2], π/20, π/2)
@show θᵐᵃˣ*180/π

my_plot(50, θᵐᵃˣ .* collect(0.5:0.1:1.5), .01)

✔️✔️✔️(θᵐᵃˣ * 180) / π = 

40.94384355985676