{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "2e62d7a8",
"metadata": {},
"source": [
"# Cours ENPC - Pratique du calcul scientifique"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "e56242bb",
"metadata": {},
"source": [
"## Examen final\n",
"\n",
"- Ce notebook est à soumettre sur Educnet avant 16h\n",
" (ou 16h40 pour les quelques étudiants bénéficiant de temps supplémentaire).\n",
"\n",
"- L’examen comporte trois exercices indépendants. Dans chaque exercice les\n",
" cellules peuvent éventuellement dependre des cellules précèdentes.\n",
"\n",
"- Afin de faciliter l'évaluation de votre code,\n",
" ne pas changer les signatures des fonctions à implémenter.\n",
"\n",
"- La cellulle ci-dessous importe les packages utilisés dans ce notebook et\n",
" définit une macro qui a pour but de faciliter les tests unitaires figurant\n",
" dans le sujet. Il est ainsi nécessaire d'exécuter le code dans cette\n",
" cellule avant de poursuivre le notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5c96ae5c",
"metadata": {},
"outputs": [],
"source": [
"using Polynomials\n",
"using Plots\n",
"using LaTeXStrings\n",
"\n",
"macro mark(bool_expr)\n",
" return :( print($bool_expr ? \"✔️\" : \"❌\"))\n",
"end"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "990af1ce",
"metadata": {},
"source": [
"### 1. Exercice sur l'équation de Laplace\n",
"\n",
"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:\n",
"\n",
"$$ 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.$$\n",
"Pour cela, on approxime $\\varphi$ avec un polynome interpolateur $\\widehat \\varphi$,\n",
"puis on résout exactement l'équation de Laplace avec membre de droite $\\widehat \\varphi$ au lieu de $\\varphi$.\n",
"\n",
"1. Écrire une fonction `approx(n)` implémentant cette approche.\n",
" 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.\n",
" On prendra comme membre de droite la fonction\n",
" $$\\varphi(x) = \\exp\\Bigl(\\sin(2\\pi x)\\Bigr) \\Bigl(\\sin(2\\pi x)-\\cos(2\\pi x)^2 \\Bigr),$$\n",
" auquel cas la solution analytique est donnée par\n",
" $$ u(x)=(2\\pi)^{-2}\\Bigl(\\exp\\bigl(\\sin(2\\pi x)\\bigr)-1\\Bigr).$$"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "4a8bb0cc",
"metadata": {},
"source": [
" > **Indications :**\n",
" > - Utiliser la fonction `fit` de la bibliothèque `Polynomials.jl` pour obtenir le polynôme interpolateur:\n",
" > ```julia\n",
" > p = fit(x,y)\n",
" > ```\n",
" > où `x` sont les noeuds d'interpolation et `y` sont les valeurs de la fonction à interpoler.\n",
" >\n",
" > - 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$.\n",
" >\n",
" > - On pourra utiliser la fonction `integrate` de la bibliothèque `Polynomials.jl` qui calcule une primitive d'un polynôme:\n",
" > ```julia\n",
" > P = integrate(p)\n",
" > ```\n",
" > - Utiliser le format `BigFloat` pour limiter les erreurs d'arrondi.\n",
" > 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`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f8f96ff5",
"metadata": {},
"outputs": [],
"source": [
"φ(x) = exp(sin(2π*x))*(sin(2π*x)-cos(2π*x)^2) # \\varphi + TAB pour écrire φ\n",
"u(x) = (exp(sin(2π*x))-1)/(2π)^2\n",
"\n",
"function approx(n)\n",
"\n",
"end\n",
"\n",
"@mark typeof(approx(3)) == Polynomial{BigFloat, :x}\n",
"@mark degree(approx(3)) == 5\n",
"@mark approx(3)(0) ≈ 0\n",
"@mark approx(3)(1) ≈ 0"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "d33d12ab",
"metadata": {},
"source": [
"2. Faire une animation permettant de visualiser la solution exacte et la solution approchée pour des valeurs de `n` allant de 2 à 20,\n",
" en utilisant la fonction `animate` de la bibliothèque `Plots`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "48d67ca8",
"metadata": {},
"outputs": [],
"source": []
},
{
"attachments": {},
"cell_type": "markdown",
"id": "9cb7efb5",
"metadata": {},
"source": [
"3. Écrire une fonction `estimate_error(n)` qui approxime l'erreur,\n",
" en norme $L^\\infty$,\n",
" entre la solution approchée par l'approche ci-dessus et la solution exacte."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9216d67c",
"metadata": {},
"outputs": [],
"source": [
"function estimate_error(n)\n",
"\n",
"end\n",
"\n",
"@mark estimate_error(2) > 0\n",
"@mark estimate_error(20) < 1e-3\n",
"@mark estimate_error(20) > 1e-8\n",
"@mark estimate_error(40) < 1e-6"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "531f17be",
"metadata": {},
"source": [
"4. Tracer un graphique de l'erreur en fonction de $n$ pour $n$ allant de 5 à 50.\n",
" Utiliser l'échelle par défaut pour l'axe des abcisses et une échelle logarithmique pour l'axe des ordonnées."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "68ae4534",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": []
},
{
"attachments": {},
"cell_type": "markdown",
"id": "586a7cf5",
"metadata": {},
"source": [
"### 2. Calcul des racines cubiques de 1\n",
"\n",
"On souhaite calculer numériquement les solutions dans $\\mathbb C$ de l'équation\n",
"$$\n",
"f(z) = 0, \\qquad f(z) := z^3 - 1.\n",
"$$\n",
"Pour ce faire, on utilisera l'itération de Newton-Raphson\n",
"$$\n",
"z_{k+1} = z_k - \\frac{f(z_k)}{f'(z_k)},\n",
"\\tag{NR}\n",
"$$\n",
"\n",
"où $f'\\colon \\mathbb C \\to \\mathbb C$ est la dérivée complexe de la fonction $f$,\n",
"ici donnée par $f'(z) = 3z^2$.\n",
"\n",
"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₀`,\n",
"ou `nothing` si une solution n'a pas été trouvée après `maxiter` itérations.\n",
"L'argument `df` est la dérivée complexe de la fonction `f`,\n",
"et on utilisera comme critère d'arrêt\n",
"$$\n",
"|f(z_k)| ≤ \\varepsilon.\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8ff73752",
"metadata": {},
"outputs": [],
"source": [
"function newton_raphson(f, df, z₀; maxiter = 100, ε = 1e-12)\n",
"\n",
"end\n",
"\n",
"@mark newton_raphson(z -> z^2 - 2, z -> 2z, 1) ≈ √2\n",
"@mark newton_raphson(z -> z^2 - 2, z -> 2z, -1) ≈ -√2\n",
"@mark newton_raphson(z -> z^3 - 2, z -> 3z^2, 1) ≈ cbrt(2)\n",
"@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1 + im) ≈ im\n",
"@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1 - im) ≈ -im\n",
"@mark newton_raphson(z -> z^2 + 1, z -> 2z, 1) == nothing"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "7a36a398",
"metadata": {},
"source": [
"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.\n",
" En vue de calculer numériquement les bassins d'attraction des trois racines,\n",
" écrire une fonction `which_basin(z₀)` qui, pour un argument `z₀`,\n",
"\n",
" - renvoit 0 si la méthode initialisée en `z₀` ne converge pas;\n",
"\n",
" - renvoit 1 si la méthode initialisée en `z₀` converge vers $1$;\n",
"\n",
" - renvoit 2 si la méthode initialisée en `z₀` converge vers $\\exp(2i\\pi/3)$;\n",
"\n",
" - renvoit 3 si la méthode initialisée en `z₀` converge vers $\\exp(4i\\pi/3)$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0b7a9a6c",
"metadata": {},
"outputs": [],
"source": [
"f(z) = z^3 - 1\n",
"df(z) = 3z^2\n",
"\n",
"function which_root(z₀)\n",
"\n",
"end\n",
"\n",
"@mark which_root(1) == 1\n",
"@mark which_root(im) == 2\n",
"@mark which_root(-im) == 3"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "66c71f90",
"metadata": {},
"source": [
"3. Calculer numériquement et représenter graphiquement les bassins d'attraction des trois racines\n",
" dans le domaine $[-2, 2] \\times [-2, 2]$ du plan complexe.\n",
" La fonction `heatmap` avec argument `levels=4` pourra étre utilisée pour la représentation graphique.\n",
" Cette fonction s'utilise de la même manière que les fonctions `contour` et `contourf`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "afea2fb6",
"metadata": {},
"outputs": [],
"source": []
},
{
"attachments": {},
"cell_type": "markdown",
"id": "106cdaec",
"metadata": {},
"source": [
"4. On considère à présent une autre approche,\n",
" appelée *méthode de la sécante*, pour le calcul des racines cubiques complexes de 1.\n",
" Une itération de cette méthode est donnée par\n",
" $$\n",
" z_{k+2} = \\frac{z_{k} f(z_{k+1}) - z_{k+1}f(z_k)}{f(z_{k+1}) - f(z_{k})}.\n",
" $$\n",
" Écrire une fonction `secant(f, z₀, z₁; maxiter = 100, ε = 1e-12)` qui,\n",
" contrairement à la méthode `newton_raphson` ci-dessus,\n",
" 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₁`,\n",
" ou `nothing` si une solution n'a pas été trouvée après `maxiter` itérations.\n",
" On utilisera comme critère d'arrêt\n",
" $$\n",
" |f(z_k)| ≤ \\varepsilon.\n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "48ebd3da",
"metadata": {},
"outputs": [],
"source": [
"function secant(f, z₀, z₁; maxiter = 100, ε = 1e-12)\n",
"\n",
"end\n",
"\n",
"@mark secant(z -> z^2 - 2, 1., 2.)[1:2] ≈ [1. ; 2.]\n",
"@mark secant(z -> z^2 - 2, 1., 2.)[end] ≈ √2\n",
"@mark secant(z -> z^2 - 2, -1., -2.)[end] ≈ -√2\n",
"@mark secant(z -> z^2 + 1, 1. + im, 1. + 2im)[end] ≈ im\n",
"@mark secant(z -> z^2 + 1, 1. + im, 1. - 2im)[end] ≈ -im\n",
"@mark secant(z -> z^2 + 1, 1., 2.) == nothing"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "447a24c6",
"metadata": {},
"source": [
"5. Il est possible de démontrer que si cette itération converge vers une solution $r$, alors\n",
" $$\n",
" \\tag{Q}\n",
" \\lim_{k \\to \\infty} \\frac{|z_{k+1} - r|}{|z_k - r|^q} \\in \\mathbb R_{>0},\n",
" $$\n",
" \n",
" pour un réel positif $q$ qui est l'*ordre de convergence*.\n",
" Calculer empiriquement dans la fonction `estimate_q` la valeur de $q$\n",
" en appliquant à la fonction $z \\mapsto z^3 - 1$ la méthode de la sécante initialisée avec\n",
" $$\n",
" z_0 = -1, \\qquad z_1 = - i.\n",
" $$\n",
" 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.\n",
" C'est ce qui est fait dans le début de code qui vous est donné ci-dessous.\n",
" La valeur correcte de $q$ est donnée par\n",
" $$\n",
" q = \\frac{1 + \\sqrt{5}}{2} \\approx 1.618.\n",
" $$\n",
"\n",
" > *Remarque*. Soit $y_k := - \\log(|z_{k} - r|)$.\n",
" > L'équation (Q) implique que\n",
" > $$\n",
" > \\lim_{k \\to \\infty} y_{k+1} - q y_{k} = C \\in \\mathbb R.\n",
" > $$\n",
" > On en déduit que\n",
" > $$\n",
" > q = \\lim_{k\\to \\infty} \\frac{y_{k+1}}{y_k}\n",
" > $$\n",
" > Cette équation permet d'estimer $q$ à partir de $y_{k+1}$ et $y_k$\n",
" > pour $k$ suffisamment grand."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5c3c9419",
"metadata": {},
"outputs": [],
"source": [
"# Set number of significant digits in base 10\n",
"precision = 1_000\n",
"setprecision(precision, base=10)\n",
"\n",
"# Initial values for the secant iteration\n",
"z₀ = big(-1.)\n",
"z₁ = big(-1.0im)\n",
"\n",
"# Exact root towards which convergence should occur\n",
"root = exp(im * 4big(π)/3)\n",
"\n",
"function estimate_q()\n",
"\n",
"end\n",
"\n",
"@mark abs(estimate_q() - (1 + √5)/2) < 1e-2"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "30f68a98",
"metadata": {},
"source": [
"### 3. Résolution d'une équation différentielle\n",
"\n",
"On s'intéresse dans cet exercice au calcul de la trajectoire d'une balle de golf frappée par un joueur.\n",
"Pour simplifier, on supposera que la trajectoire de la balle est confinée dans un plan perpendiculaire à la surface de jeu,\n",
"qu'on supposera parfaitement horizontale.\n",
"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.\n",
"Le mouvement peut être régit par l'équation de Newton:\n",
"$$\n",
"m \\mathbf x''(t) = - \\mathbf F^d(\\lVert \\mathbf v \\rVert) \\frac{\\mathbf v}{\\lVert \\mathbf v \\rVert} - mg \\mathbf e_2, \\qquad \\mathbf v := \\mathbf x'.\n",
"\\tag{Golf}\n",
"$$\n",
"\n",
"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$,\n",
"et $\\mathbf F^d(\\lVert v \\rVert)$ est la force de traînée, qui par analyse dimensionnelle peut s'écrire\n",
"$$\n",
"\\mathbf F^d(\\lVert \\mathbf v \\rVert) = \\frac{1}{2} \\rho A C^d \\lVert \\mathbf v \\rVert^2.\n",
"$$\n",
"La signification des constantes apparaissant dans cette formule,\n",
"et les valeurs que nous prendrons dans cette exercice,\n",
"sont données ci-dessous.\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f190ce81",
"metadata": {},
"outputs": [],
"source": [
"# Mass of the golf ball [kg]\n",
"const m = .045\n",
"\n",
"# Radius of the golf ball [m]\n",
"const r = .021\n",
"\n",
"# Air density at 0 ⁰C [kg/m³]\n",
"const ρ = 1.293\n",
"\n",
"# Gravity acceleration [m/s²]\n",
"const g = 9.81\n",
"\n",
"# Cross-sectional areal [m²]\n",
"const A = π*r^2\n",
"\n",
"# Drag coefficient [dimensionless]\n",
"const Cᵈ = .2\n",
"\n",
"# Drag force depending on V := ‖𝐯‖\n",
"Fᵈ(V) = 1/2 * ρ*A*Cᵈ*V^2"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "303bef45",
"metadata": {},
"source": [
"On modélise le joueur par un point au sol, qu'on prend comme origine du repère orthonormé.\n",
"Pour modéliser le fait que la balle est frappée par le joueur au temps $t = 0$,\n",
"on prend la condition initiale suivante pour (Golf).\n",
"$$\n",
"\\mathbf x(0) = \\begin{pmatrix} 0 \\\\ 0 \\end{pmatrix},\n",
"\\qquad\n",
"\\mathbf x'(0) = V_0 \\begin{pmatrix} \\cos(\\theta) \\\\ \\sin(\\theta) \\end{pmatrix}.\n",
"\\tag{IC}\n",
"$$\n",
"\n",
"Ici $V_0$ est la vitesse initiale de la balle, et $\\theta$ est l'angle de loft."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "736c639e",
"metadata": {},
"source": [
"1. Soit $\\mathbf v = \\mathbf x(t)$.\n",
" 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:\n",
" $$\n",
" \\mathbf Z'(t) = \\mathbf f(t, \\mathbf Z), \\qquad \\mathbf Z(0) = \\mathbf Z_0.\n",
" \\tag{1st-order}\n",
" $$\n",
" \n",
" (Dans le cas qui nous occupe, la fonction $f$ est en fait indépendante du temps.)\n",
" Écrire la fonction $f$ sous forme d'une fonction Julia `f(t, Z)`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ac778bac",
"metadata": {},
"outputs": [],
"source": [
"function f(t, Z)\n",
"\n",
"end\n",
"\n",
"@mark f(1, [0; 0; 100; 0]) ≈ [100.0, 0.0, -39.8083771506977, -9.81]\n",
"@mark f(0, [100; 100; 100; 100]) ≈ [100.0, 100.0, -56.297546862579914, -66.10754686257991]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "abd956dd",
"metadata": {},
"source": [
"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)$.\n",
" Cette méthode est basée sur l'itération suivante:\n",
" $$\n",
" \\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),\n",
" $$\n",
" où\n",
" \\begin{align*}\n",
" \\mathbf k_1 &= \\ f(t_n, \\mathbf Z_n), \\\\\n",
" \\mathbf k_2 &= \\ f\\!\\left(t_n + \\frac{\\Delta}{2}, \\mathbf Z_n + \\frac{\\Delta}{2} \\mathbf k_1\\right), \\\\\n",
" \\mathbf k_3 &= \\ f\\!\\left(t_n + \\frac{\\Delta}{2}, \\mathbf Z_n + \\frac{\\Delta}{2} \\mathbf k_2\\right), \\\\\n",
" \\mathbf k_4 &= \\ f\\!\\left(t_n + \\Delta, \\mathbf Z_n + \\Delta\\mathbf k_3\\right).\n",
" \\end{align*}\n",
" La fonction devra renvoyer $\\mathbf Z_{n+1}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "614f059e",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"function rk4(tₙ, Zₙ, f, Δ)\n",
"\n",
"end\n",
"\n",
"@mark rk4(0., [0.], (t, Z) -> [1.], 1.) ≈ [1.0]\n",
"@mark rk4(1., [0.], (t, Z) -> [t], 1.) ≈ [3/2]\n",
"@mark rk4(0., [0.; 0.], (t, Z) -> [t^2; t^3], 1.) ≈ [1/3; 1/4]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "a29612a7",
"metadata": {},
"source": [
"3. Écrire une fonction `solve_ode(V₀, θ, Δ)` pour résoudre (1st-order) pour une vitesse initiale `V₀` et un angle initial `θ`,\n",
" en utilisant la méthode de Runge-Kutta d'ordre 4 avec pas de temps fixe `Δ`.\n",
" Votre fonction devra renvoyer un vecteur de temps `ts` et un vecteur de vecteurs `Zs` contenant la solution à ces temps.\n",
" On calculera la trajectoire jusqu'à ce que la balle soit retombée sur le sol.\n",
" 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;\n",
" il faudra donc que `Zs[end-1][2]` soit positif et `Zs[end][2]` soit négatif.\n",
"\n",
" > **Indication**. Pour construire `Zs`, il est recommandé d'utiliser la fonction `push!`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fcef09cd",
"metadata": {},
"outputs": [],
"source": [
"function solve_ode(V₀, θ, Δ)\n",
"\n",
"end\n",
"\n",
"@mark solve_ode(50, π/4, .01) |> length == 2\n",
"@mark solve_ode(50, π/4, .01)[2][end-1][2] ≥ 0\n",
"@mark solve_ode(50, π/4, .01)[1][1:10] ≈ 0:.01:.09\n",
"@mark solve_ode(50, π/4, .01)[2][end][2] ≤ 0\n",
"@mark solve_ode(50, π/4, .01)[2][end][1] ≈ 149.30535166172655"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "296613b2",
"metadata": {},
"source": [
"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`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8b08a2da",
"metadata": {},
"outputs": [],
"source": [
"function my_plot(V₀, θs, Δ)\n",
"\n",
"end\n",
"my_plot(50, [π/8, π/4, 3π/8], .01)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "9ef74fd4",
"metadata": {},
"source": [
"5. Écrire une fonction `distance(V₀, θ, Δ)` pour calculer la distance entre le joueur et le point d'impact de la bal sur le sol.\n",
" 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,\n",
" durant lequel la hauteur de la balle passe d'une valeur positive à une valeur négative.\n",
"\n",
" > **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$)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "404db73c",
"metadata": {},
"outputs": [],
"source": [
"function distance(V₀, θ, Δ)\n",
"\n",
"end\n",
"\n",
"@mark distance(50, π/8, .01) ≈ 124.20307897577761\n",
"@mark distance(50, π/4, .01) ≈ 149.29126957486164\n",
"@mark distance(50, 3π/8, .01) ≈ 102.42272578334352\n",
"@mark distance(50, π/2 - .001, .01) < 1."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "f3ea670e",
"metadata": {},
"source": [
"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$,\n",
" estimée avec $\\Delta = .01s$.\n",
" Estimer graphiquement les angles permettant d'atteindre une distance de 80m.\n",
"\n",
" > **Indication** : il peut être utile de passer `xticks=0.:5.:90.` en argument à la fonction `plot` pour simplifier cette estimation."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ca660bd3",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"V₀, Δ = 50, .01"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "e7cc9aca",
"metadata": {},
"source": [
"7. Pour une vitesse initiale $V_0$ donnée,\n",
" on veut calculer l'angle de loft $\\theta$ permettant d'atteindre un trou situé à une distance $L$ du joueur.\n",
" Pour ce faire,\n",
" on se propose de mettre en œuvre l'algorithme de Newton-Raphson sur une fonction scalaire prenant comme argument $\\theta$\n",
" et s'annulant lorsque l'estimation de la distance obtenue par la fonction `distance` est égale à $L$.\n",
" La méthode de Newton-Raphson nécessite de connaître la dérivée de la fonction dont on cherche une racine,\n",
" ce qui peut être fait par différentiation automatique.\n",
" 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,\n",
" sauf les fonctions `cos(x::D)`, `sin(x::D)` et `sqrt(x::D)` qui sont à définir."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c090b0e2",
"metadata": {},
"outputs": [],
"source": [
"import Base: +, -, *, /, ==, ≈, sqrt, cos, sin, inv, conj, abs, isless, AbstractFloat, convert, promote_rule, show\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) = D(x.f .- y.f)\n",
"-(x::D) = D(.-(x.f))\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) = D((x.f[1]/y.f[1], (y.f[1]*x.f[2] - x.f[1]*y.f[2])/y.f[1]^2))\n",
"==(x::D, y::D) = x.f[1] == y.f[1] && x.f[2] == y.f[2]\n",
"≈(x::D, y::D) = x.f[1] ≈ y.f[1] && x.f[2] ≈ y.f[2]\n",
"abs(x::D) = D((abs(x.f[1]), x.f[2]*sign(x.f[1])))\n",
"abs2(x::D) = D((abs2(x.f[1]), 2x.f[1]*x.f[2]))\n",
"isless(x::D, y::D) = isless(x.f[1], y.f[1])\n",
"isless(x::D, y::Real) = isless(x.f[1], y)\n",
"isless(x::Real, y::D) = isless(x, y.f[1])\n",
"D(x::D) = x\n",
"AbstractFloat(x::D) = x.f[1]\n",
"convert(::Type{D}, x::Real) = D((x,zero(x)))\n",
"convert(::Type{Real}, x::D) = x.f[1]\n",
"promote_rule(::Type{D}, ::Type{<:Real}) = D\n",
"show(io::IO,x::D) = print(io,x.f[1],x.f[2]<0 ? \" - \" : \" + \",abs(x.f[2]),\" ε\")\n",
"ε = D((0, 1))\n",
"\n",
"@mark sin(π/4 + ε) ≈ √2/2 * (1. + ε)\n",
"@mark cos(ε) == 1\n",
"@mark sin(ε) ≈ ε\n",
"@mark sqrt(1 + ε) ≈ 1 + .5ε"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "dbf3d618",
"metadata": {},
"source": [
"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`,\n",
" dans laquelle la dérivée de `f` au point courant est obtenue par exploitation des nombres duaux,\n",
" 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)| < δ$.\n",
"\n",
" > **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)))`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "09023cc7",
"metadata": {},
"outputs": [],
"source": [
"function newton_raphson_dual(f, x, maxiter=100; δ = 1e-12)\n",
"\n",
"end\n",
"\n",
"@mark newton_raphson_dual(x -> x^2 - 2, 1) ≈ √2\n",
"@mark newton_raphson_dual(x -> x^2 - 2, -1) ≈ -√2\n",
"@mark newton_raphson_dual(x -> x^3 - 2, 1) ≈ cbrt(2)\n",
"@mark newton_raphson_dual(x -> cos(x) - .5, 2) ≈ acos(.5)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "016265d4",
"metadata": {},
"source": [
"9. Écrire une fonction `find_angle(L, V₀, Δ, θ₀)` renvoyant un angle qui permet d'atteindre une distance $L$ pour une vitesse initiale $V_0$.\n",
" On supposera que les arguments de la fonction sont tels qu'un tel angle existe.\n",
" Calculer une estimation des deux angles `θ₁` et `θ₂` permettant d'atteindre une distance de 80m et tracer les trajectoires correspondantes."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "17815cc2",
"metadata": {
"lines_to_next_cell": 0
},
"outputs": [],
"source": [
"function find_angle(L, V₀, Δ, θ₀)\n",
"\n",
"end\n",
"\n",
"# Calculer ici les deux angles permettant d'atteindre une distance de 80m.\n",
"# puis tracer les trajectoires correspondantes."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "edee5b8c",
"metadata": {},
"source": [
"10. Ayant constaté graphiquement à la question 6 l'existence d'un angle maximisant la portée,\n",
" 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$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "58a8e229",
"metadata": {},
"outputs": [],
"source": [
"function bisection(f, a, b; δ = 1e-10)\n",
"\n",
"end\n",
"\n",
"@mark bisection(x -> x^2 - 2, 1, 2) ≈ √2\n",
"@mark bisection(x -> x^3 - 2, 1, 2) ≈ cbrt(2)\n",
"@mark bisection(x -> cos(x) - .5, 1, 2) ≈ acos(.5)\n",
"\n",
"V₀, Δ = 50, .01\n",
"θᵐᵃˣ = # calculer ici l'angle associé à la portée maximale\n",
"@mark abs(θᵐᵃˣ*180/π - 40.94) < .1"
]
}
],
"metadata": {
"jupytext": {
"encoding": "# -*- coding: utf-8 -*-"
},
"kernelspec": {
"display_name": "Julia 1.8.5",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}