{ "cells": [ { "cell_type": "markdown", "id": "bb5469c0", "metadata": {}, "source": [ "# Cours ENPC - Pratique du calcul scientifique" ] }, { "cell_type": "markdown", "id": "9b241907", "metadata": {}, "source": [ "## Examen final\n", "\n", "- Ce notebook est à soumettre sur Educnet avant 16h30.\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." ] }, { "cell_type": "code", "execution_count": null, "id": "f484a3f0", "metadata": {}, "outputs": [], "source": [ "using Polynomials\n", "using Plots\n", "using LaTeXStrings\n", "using LinearAlgebra\n", "\n", "macro mark(bool_expr)\n", " return :(print($bool_expr ? \"✔️\" : \"❌\"))\n", "end" ] }, { "cell_type": "markdown", "id": "ba6bb13e", "metadata": {}, "source": [ "### 1. Intégration numérique" ] }, { "cell_type": "markdown", "id": "ef8090dd", "metadata": {}, "source": [ "1. Écrire une fonction `composite_trapezoidal(u, a, b, n)` permettant d'approximer l'intégrale\n", " $$\n", " I := \\int_a^b u(x) \\, \\mathrm d x\n", " $$\n", " par la méthode trapézoidale composite avec `n` de points équidistants $a = x_1 < x_2 < \\dots < x_{n-1} < x_n = b$.\n", " On supposera que $n \\geq 2$." ] }, { "cell_type": "code", "execution_count": null, "id": "1b35f128", "metadata": {}, "outputs": [], "source": [ "function composite_trapezoidal(u, a, b, n)\n", "\n", "end\n", "\n", "@mark composite_trapezoidal(x -> 5, 1, 2, 100) ≈ 5\n", "@mark composite_trapezoidal(x -> x, 1, 2, 100) ≈ 3/2\n", "@mark composite_trapezoidal(x -> x, 1, 2, 2) ≈ 3/2\n", "@mark composite_trapezoidal(x -> x^2, 0, 1, 2) ≈ 1/2\n", "@mark composite_trapezoidal(x -> x^2, 1, 2, 2) ≈ 5/2" ] }, { "cell_type": "markdown", "id": "3e79ef6c", "metadata": {}, "source": [ "2. Écrire une fonction `composite_simpson(u, a, b, n)` permettant d'approximer l'intégrale $I$ par la méthode de Simpson composite\n", " basée sur des évaluations de `u` à un nombre **impair** `n` de points équidistants tels que $a = x_1 < x_2 < \\dots < x_{n-1} < x_n = b$.\n", " On supposera que $n$ est impair et $n \\geq 3$.\n", "\n", " > **Remarque**: `n` est ici le nombre de points auxquels la fonction `u` est évaluée,\n", " > et pas un nombre d'intervalles où la règle de Simpson est appliquée localement." ] }, { "cell_type": "code", "execution_count": null, "id": "7dda7429", "metadata": {}, "outputs": [], "source": [ "function composite_simpson(u, a, b, n)\n", " @assert n % 2 == 1 \"`n` must be impair\"\n", "\n", "end\n", "\n", "@mark composite_simpson(x -> 1 , 1, 2, 101) ≈ 1\n", "@mark composite_simpson(x -> x , 1, 2, 101) ≈ 3/2\n", "@mark composite_simpson(x -> x^2, 1, 2, 101) ≈ 7/3\n", "@mark composite_simpson(x -> x^3, 1, 2, 101) ≈ 15/4\n", "@mark composite_simpson(x -> x , 0, 1, 3) ≈ 1/2\n", "@mark composite_simpson(x -> x^2, 0, 1, 3) ≈ 1/3\n", "@mark composite_simpson(x -> x^3, 0, 1, 3) ≈ 1/4" ] }, { "cell_type": "markdown", "id": "0070aaa8", "metadata": {}, "source": [ "3. Écrire une fonction `calculate_sum(N)` qui permet de calculer la somme\n", " $$\n", " S(n) := \\sum_{n = 1}^{N} n^{-n}.\n", " $$\n", " Afficher la valeur de $S(N)$ pour $n$ égal à 5, 10, 15, et 20." ] }, { "cell_type": "code", "execution_count": null, "id": "cd830922", "metadata": {}, "outputs": [], "source": [ "function calculate_sum(N)\n", "\n", "end\n", "\n", "@mark abs(calculate_sum(20) - 1.2912859970626636) < 1e-6\n", "@mark abs(calculate_sum(20) - 1.2912859970626636) < 1e-9\n", "@mark abs(calculate_sum(20) - 1.2912859970626636) < 1e-12" ] }, { "cell_type": "markdown", "id": "d5fa0650", "metadata": {}, "source": [ "4. On peut montrer que\n", " $$\n", " \\int_0^1 x^{-x} \\, \\mathrm d x = \\sum_{n=1}^{\\infty} n^{-n}.\n", " $$\n", " Illustrer l'erreur des méthodes composites définies ci-dessus en fonction de `n`,\n", " sur un même graphe.\n", " Utiliser $S(20)$ comme valeur de référence pour l'intégrale,\n", " et employer l'échelle logarithmique pour les deux axes du graphe.\n", "\n", " > **Remarque**: La fonction à intégrer dans cet exercice est continue,\n", " > mais sa dérivée diverge en $x = 0$.\n", " > Ne vous inquiétez donc pas si le taux de convergence que vous observez ne correspond pas au taux théorique." ] }, { "cell_type": "code", "execution_count": null, "id": "f67d698d", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "49b58aff", "metadata": {}, "source": [ "5. (**Bonus**). Estimer, en approximant à l'aide de la fonction `fit` le logarithme de l'erreur par une fonction affine du logarithme du pas d'intégration,\n", "l'ordre de convergence de la méthode composite de Simpson pour le calcul de l'intégrale dans la question précédente." ] }, { "cell_type": "code", "execution_count": null, "id": "59b31df7", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "715b2394", "metadata": {}, "source": [ "### 2. Résolution d'un système linéaire" ] }, { "cell_type": "markdown", "id": "e89f16dd", "metadata": {}, "source": [ "L'objectif de cet exercice est de proposer un algorithme permettant de réaliser la décomposition LU d'une matrice réelle $\\mathsf{A}\\in\\mathbb{R}^{n×n}$,\n", "**non pas par élimination gaussienne** mais par identification des entrées de $\\mathsf A$ avec celles de $\\mathsf L \\mathsf U$.\n", "Il s'agit de trouver un matrice triangulaire inférieure $\\mathsf L$ formée de 1 sur la diagonale\n", "et une matrice triangulaire supérieure $\\mathsf U$ telles que :\n", "\n", "$$\n", "\\tag{LU}\n", "\\mathsf{A}=\\mathsf{L}\\mathsf{U}\n", "$$\n", "\n", "1. Écrire une fonction `my_lu(A)` qui prend comme argument une matrice `A` et qui renvoie les matrices `L` et `U`.\n", " Pour calculer ces matrices, s'appuyer sur une identification successive des éléments des deux membres de (LU),\n", " ligne par ligne de haut en bas, et de gauche à droite au sein de chaque ligne.\n", "\n", " > **Indication**: lorsqu'on suit l'ordre conseillé,\n", " > la comparaison de l'élément $(i, j)$ fournit une équation pour $\\ell_{ij}$ si $j < i$,\n", " > et une équation pour $u_{ij}$ si $j \\geq i$.\n", " > Notons qu'il est possible de parcourir les éléments dans d'autres ordres." ] }, { "cell_type": "code", "execution_count": null, "id": "77163f0d", "metadata": {}, "outputs": [], "source": [ "function my_lu(A)\n", "\n", "end\n", "\n", "@mark my_lu(diagm([1; 2; 3])) == (diagm([1; 1; 1]), diagm([1; 2; 3]))\n", "@mark my_lu([2 -1 0; -1 2 -1; 0 -1 2])[1] ≈ [1 0 0; -1/2 1 0; 0 -2/3 1]\n", "@mark my_lu([2 -1 0; -1 2 -1; 0 -1 2])[2] ≈ [2 -1 0; 0 3/2 -1; 0 0 4/3]\n", "@mark begin C = [1 2 3 4; 4 3 2 1; 1 2 1 2; 1 5 4 1]; my_lu(C)[1] ≈ lu(C, NoPivot()).L end\n", "@mark begin C = [1 2 3 4; 4 3 2 1; 1 2 1 2; 1 5 4 1]; my_lu(C)[2] ≈ lu(C, NoPivot()).U end\n", "@mark begin C = randn(100, 100); my_lu(C)[1] ≈ lu(C, NoPivot()).L end\n", "@mark begin C = randn(100, 100); my_lu(C)[2] ≈ lu(C, NoPivot()).U end" ] }, { "cell_type": "markdown", "id": "6104ef50", "metadata": {}, "source": [ "2. On suppose maintenant que la matrice réelle définie positive `A` est à largeur de bande `b` supposée beaucoup plus petite que `n`.\n", " Réécrire la fonction de décomposition LU en exploitant la largeur de bande." ] }, { "cell_type": "code", "execution_count": null, "id": "a2d83280", "metadata": {}, "outputs": [], "source": [ "function my_banded_lu(A, b)\n", "\n", "end\n", "\n", "@mark begin C = [1 2 3 4; 4 3 2 1; 1 2 1 2; 1 5 4 1]; my_lu(C)[1] ≈ lu(C, NoPivot()).L end\n", "@mark begin C = [1 2 3 4; 4 3 2 1; 1 2 1 2; 1 5 4 1]; my_lu(C)[2] ≈ lu(C, NoPivot()).U end\n", "@mark begin C = randn(100, 100); my_banded_lu(C, 100)[1] ≈ lu(C, NoPivot()).L end\n", "@mark begin C = randn(100, 100); my_banded_lu(C, 100)[2] ≈ lu(C, NoPivot()).U end" ] }, { "cell_type": "markdown", "id": "e1a83dbd", "metadata": {}, "source": [ "3. Construire une fonction `generate_banded(n, b)` permettant de générer une matrice carrée aléatoire de taille `n` à largeur de bande donnée `b`." ] }, { "cell_type": "code", "execution_count": null, "id": "e3ebaadc", "metadata": {}, "outputs": [], "source": [ "function generate_banded(n, b)\n", "\n", "end\n", "\n", "@mark generate_banded(10, 2)[1, 5] == 0\n", "@mark generate_banded(10, 2)[2, 5] == 0\n", "@mark generate_banded(10, 2)[3, 5] != 0\n", "@mark generate_banded(10, 2)[4, 5] != 0\n", "@mark generate_banded(10, 2)[5, 5] != 0\n", "@mark generate_banded(10, 2)[6, 5] != 0\n", "@mark generate_banded(10, 2)[7, 5] != 0\n", "@mark generate_banded(10, 2)[8, 5] == 0\n", "@mark generate_banded(10, 2)[9, 5] == 0" ] }, { "cell_type": "markdown", "id": "d307d048", "metadata": {}, "source": [ "4. En utilisant `generate_banded`, tester votre implémentation de `my_banded_lu`,\n", " pour `n = 100` et des valeurs de `b` égales à 2, 3 et 10.\n", " Utiliser la fonction `lu` de la bibliothèque `LinearAlgebra`,\n", " avec l'argument `NoPivot()`, comme fonction de référence.\n", " Vous pouvez également utiliser la macro `@mark` pour vos tests." ] }, { "cell_type": "code", "execution_count": null, "id": "647e022b", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "a505537e", "metadata": {}, "source": [ "### 3. Résolution d'une équation différentielle\n", "\n", "Cet exercice vise à calculer la trajectoire d'une petite fusée\n", "et à dimensionner son chargement en carburant en vue d'atteindre un certain objectif.\n", "On fera plusieurs hypothèses simplificatrices :\n", "\n", "- On néglige la rotation de la terre.\n", "\n", "- On néglige la variation de l'accélération de gravité $g$ et de la densité de l'air $\\rho$ avec l'altitude.\n", "\n", "- On suppose que le coefficient de traînée $C^d$ est indépendant du nombre de Reynolds de l'écoulement.\n", "\n", "- La mouvement de la fusée est confiné à l'axe verticale. Son altitude et sa vitesse au lancement sont $z = 0$ et $v = 0$.\n", "\n", "- La fusée est approximée par un cylindre de rayon $r$ (sa hauteur n'a pas d'importance pour cet exercice).\n", "\n", "- Le carburant est éjecté à une vitesse constante par rapport à la fusée, notée $V_e$.\n", "\n", "- Le carburant est consommé à un taux $\\beta(\\mu)$,\n", " dépendant uniquement de la masse $\\mu$ de carburant restant.\n", "\n", "On note $z(t)$ l'altitude de la fusée, $v(t)$ sa vitesse et $\\mu(t)$ la masse de carburant restant.\n", "Sous les hypothèses que nous avons faites,\n", "le mouvement de la fusée peut-être modélisé par le système d'équations différentielles suivant:\n", "$$\n", "\\tag{Rocket}\n", "\\left\\{\n", "\\begin{aligned}\n", "z'(t) &= v(t), \\\\\n", "m(t) v'(t) &= \\beta\\bigl(\\mu(t)\\bigr) V_e - F^d\\bigl(v(t)\\bigr) - m(t) g, \\\\\n", "\\mu'(t) &= - \\beta\\bigl(\\mu(t)\\bigr).\n", "\\end{aligned}\n", "\\right.\n", "\\qquad\n", "\\left\\{\n", "\\begin{aligned}\n", "z(0) &= 0, \\\\\n", "v(0) &= 0, \\\\\n", "\\mu(0) &= \\mu_0.\n", "\\end{aligned}\n", "\\right.\n", "$$\n", "\n", "Ici, $\\mu_0$ est la masse de carburant au lancement,\n", "et $m(t) = m_r + \\mu(t)$ est la masse totale de la fusée,\n", "qui comporte une partie fixe $m_r$, correspondant à la structure et à la cargaison,\n", "et une partie variable $\\mu(t)$ correspondant au carburant.\n", "L'expression de la force de trainée $F^d$ est donnée dans la cellule ci-dessous.\n", "On fera varier au cours de l'exercice les paramètres $\\mu_0$ et $C^d$,\n", "et ceux-ci ne sont donc pas définis ci-dessous." ] }, { "cell_type": "code", "execution_count": null, "id": "fb062f97", "metadata": {}, "outputs": [], "source": [ "\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", "# Radius of the rocket [m]\n", "const r = .1\n", "\n", "# Cross-sectional area [m²]\n", "const A = π*r^2\n", "\n", "# Effective exhaust velocity [m/s]\n", "const Vₑ = 1000\n", "\n", "# Mass of the rocket without fuel [kg]\n", "const mᵣ = 5\n", "\n", "# Burn rate in the limit where μ → ∞ [kg/s]\n", "const β₊ = 1\n", "\n", "# Burn rate function [kg/s]\n", "β(μ) = μ > 0. ? β₊ * tanh(μ) : 0.\n", "\n", "# Drag force depending on v\n", "Fᵈ(v, Cᵈ) = 1/2 * ρ*A*Cᵈ*v^2" ] }, { "cell_type": "markdown", "id": "aff9b384", "metadata": { "lines_to_next_cell": 2 }, "source": [ "1. Le problème (Rocket) décrit peut être réécrit comme un problème aux valeurs initiales du premier ordre pour le vecteur $\\mathbf X := (z, v, \\mu)^T$ de la forme suivante:\n", " $$\n", " \\mathbf X'(t) = \\mathbf f\\bigl(t, \\mathbf X(t), C^d \\bigr), \\qquad \\mathbf X(0) = \\mathbf X_0.\n", " \\tag{1st-order}\n", " $$\n", " \n", " Écrire la fonction $f$ sous forme d'une fonction Julia `f(t, X, Cᵈ)`." ] }, { "cell_type": "code", "execution_count": null, "id": "fdee7209", "metadata": {}, "outputs": [], "source": [ "function f(t, X, Cᵈ)\n", "\n", "end\n", "\n", "@mark f(0, [0, 0, 0], 0) == [0; -g; 0]\n", "@mark f(1, [0, 0, 0], 0) == [0; -g; 0]\n", "@mark f(1, [0, 0, 0], 5) == [0; -g; 0.]\n", "@mark f(0, [0, 0, 1], 0) ≈ [0.; Vₑ * β(1) / (mᵣ + 1) - g; - β(1)]\n", "@mark f(0, [0, 100, 5], 0) ≈ [100.; Vₑ * β(5) / (mᵣ + 5) - g; - β(5)]\n", "@mark f(1, [5, 5, 5], 1) ≈ [5; Vₑ * β(5) / (mᵣ + 5) - Fᵈ(5, 1) / (mᵣ + 5) - g; - β(5)]" ] }, { "cell_type": "markdown", "id": "5fc0c556", "metadata": {}, "source": [ "2. Écrire une fonction `rkx(tₙ, Xₙ, f, Δ)` implémentant un pas de temps de taille $\\Delta$ de la méthode suivante de Runge-Kutta suivante pour une équation différentielle générique de la forme $X' = h(t, X)$:\n", " $$\n", " \\mathbf X_{n+1} = \\mathbf X_n + \\frac{\\Delta}{9}\\left(2\\mathbf k_1 + 3\\mathbf k_2 + 4\\mathbf k_3 \\right),\n", " $$\n", " où\n", " \\begin{align*}\n", " \\mathbf k_1 &= \\ h(t_n, \\mathbf X_n), \\\\\n", " \\mathbf k_2 &= \\ h\\!\\left(t_n + \\frac{\\Delta}{2}, \\mathbf X_n + \\frac{\\Delta}{2} \\mathbf k_1\\right), \\\\\n", " \\mathbf k_3 &= \\ h\\!\\left(t_n + \\frac{3\\Delta}{4}, \\mathbf X_n + \\frac{3\\Delta}{4} \\mathbf k_2\\right).\n", " \\end{align*}\n", " La fonction devra renvoyer $\\mathbf X_{n+1}$." ] }, { "cell_type": "code", "execution_count": null, "id": "7fdeb373", "metadata": {}, "outputs": [], "source": [ "function rkx(tₙ, Xₙ, h, Δ)\n", "\n", "end\n", "\n", "@mark rkx(0., [0.], (t, X) -> [1.], 1.) ≈ [1]\n", "@mark rkx(1., [0.], (t, X) -> [t], 1.) ≈ [3/2]\n", "@mark rkx(0., [0.; 0.; 0.], (t, X) -> [1, t, t^2], 1.) ≈ [1; 1/2; 1/3]" ] }, { "cell_type": "markdown", "id": "e3bb48d9", "metadata": {}, "source": [ "3. Écrire une fonction `solve_ode(Δ, Cᵈ, μ₀)` permettant de résoudre (1st-order) pour les paramètres $C^d$ et $μ_0$ donnés en arguments,\n", " en utilisant la méthode de Runge-Kutta donnée ci-dessus avec pas de temps fixe `Δ`.\n", " Votre fonction devra renvoyer un vecteur de temps `ts` et un vecteur de vecteurs `Xs` contenant la solution à ces temps.\n", " On calculera la trajectoire de la fusée jusqu'à la fin de son ascension.\n", " Plus précisément, on demande d'interrompre l'intégration numérique dès que la valeur de $v$ sera devenue strictement négative;\n", " il faudra donc que `Xs[end-1][2]` soit positif et `Xs[end][2]` soit strictement négatif." ] }, { "cell_type": "code", "execution_count": null, "id": "7f4bb488", "metadata": {}, "outputs": [], "source": [ "function solve_ode(Δ, Cᵈ, μ₀)\n", "\n", "end\n", "\n", "@mark solve_ode(.01, 0, 5) |> length == 2\n", "@mark solve_ode(.01, 0, 5)[2][end-1][2] ≥ 0\n", "@mark solve_ode(.01, 0, 5)[2][end][2] ≤ 0\n", "@mark solve_ode(.01, 0, 5)[1][1:10] ≈ 0:.01:.09" ] }, { "cell_type": "markdown", "id": "143894c7", "metadata": {}, "source": [ "4. Écrire une fonction `plot_altitude(Δ, Cᵈ, μ₀)` permettant d'illustrer sur un même graphe l'altitude de la fusée en fonction du temps,\n", "pour **une** valeur de `Cᵈ` donnée et **plusieurs** valeurs de $\\mu_0$ dans le vecteur `μs`." ] }, { "cell_type": "code", "execution_count": null, "id": "6cae530a", "metadata": {}, "outputs": [], "source": [ "function plot_altitude(Δ, Cᵈ, μs)\n", " p = plot(title=\"Altitude of the rocket\")\n", "\n", " return p\n", "end\n", "\n", "Δ, Cᵈ, μs = .01, .75, [1; 2; 3; 4; 5]\n", "plot_altitude(Δ, Cᵈ, μs)" ] }, { "cell_type": "markdown", "id": "2be940c3", "metadata": {}, "source": [ "5. Écrire une fonction `plot_velocity(Δ, Cᵈ, μ₀)` permettant d'illustrer sur un même graphe la vitesse de la fusée en fonction du temps,\n", "pour **une** valeur de $Cᵈ$ donnée et **plusieurs** valeurs de $\\mu_0$ dans le vecteur `μs`." ] }, { "cell_type": "code", "execution_count": null, "id": "a2677b77", "metadata": {}, "outputs": [], "source": [ "function plot_velocity(Δ, Cᵈ, μs)\n", " p = plot(title=\"Velocity of the rocket\")\n", "\n", " return p\n", "end\n", "\n", "Δ, Cᵈ, μ₀ = .01, .75, [1; 2; 3; 4; 5]\n", "plot_velocity(Δ, Cᵈ, μ₀)" ] }, { "cell_type": "markdown", "id": "07fd582f", "metadata": {}, "source": [ "6. On suppose ici que $C^d = 0$.\n", " Dans ce cas, une équation fondamentale de l'astronautique connue sous le nom d'équation de Tsiolkovski\n", " donne la vitesse de la fusée en fonction de sa masse:\n", " $$\n", " v(t) = V_e \\log\\left(\\frac{m(0)}{m(t)}\\right) - g t.\n", " $$\n", " Il suffit donc de connaître la masse de la fusée à un instant donné pour en connaître sa vitesse.\n", " Or, la troisième équation du système (Rocket) peut être résolue analytiquement:\n", " $$\n", " \\mu(t) = \\sinh^{-1} \\Bigl( \\exp\\bigl( \\log(\\sinh(\\mu_0)) - t \\bigr) \\Bigr).\n", " $$\n", " Il est donc possible d'obenir une expression analytique de la fonction vitesse $v(t)$.\n", " Écrire une fonction `error_velocity(Δ, μ₀)` qui renvoit l'erreur en norme maximum entre cette fonction et la composante vitesse de la solution numérique,\n", " définie comme\n", " $$\n", " e(\\Delta) := \\max_{i} \\bigl\\lvert v(t_i) - \\widehat v_i \\bigr\\rvert,\n", " $$\n", " où $\\widehat v_i$ est l'approximation numérique de la vitesse au temps $i \\Delta$." ] }, { "cell_type": "code", "execution_count": null, "id": "17541a07", "metadata": {}, "outputs": [], "source": [ "μ_exact(t, μ₀) = asinh(exp(log(sinh(μ₀)) - t))\n", "v_exact(t, μ₀) = Vₑ * log((mᵣ + μ₀) / (mᵣ + μ_exact(t, μ₀))) - g*t\n", "\n", "function error_velocity(Δ, μ₀)\n", "\n", "end\n", "\n", "@mark error_velocity(.1, 5) < 1e-2\n", "@mark error_velocity(.1, 5) > 1e-3\n", "@mark error_velocity(.01, 5) < 1e-5\n", "@mark error_velocity(.01, 5) > 1e-6\n", "@mark error_velocity(.001, 5) < 1e-8\n", "@mark error_velocity(.001, 5) > 1e-9" ] }, { "cell_type": "markdown", "id": "d0ddc38f", "metadata": {}, "source": [ "7. Tracer un graphique de l'erreur en fonction de $Δ$ pour les valeurs de $Δ$ données ci-dessous,\n", " et pour $μ₀ = 5$.\n", " Utiliser l'échelle logarithmique pour les deux axes." ] }, { "cell_type": "code", "execution_count": null, "id": "1127d708", "metadata": {}, "outputs": [], "source": [ "Δs, μ₀ = 2.0 .^ (-10:-2), 5" ] }, { "cell_type": "markdown", "id": "6209ec35", "metadata": {}, "source": [ "8. Écrire une fonction `max_altitude(Δ, Cᵈ, μ₀)` renvoyant l'altitude maximale de la fusée pour une quantité de carburant $μ₀$ donnée,\n", " approximée en utilisant un pas de temps `Δ` pour les paramètres `Cᵈ` et `μ₀`." ] }, { "cell_type": "code", "execution_count": null, "id": "9ee15d90", "metadata": {}, "outputs": [], "source": [ "function max_altitude(Δ, Cᵈ, μ₀)\n", "\n", "end\n", "\n", "@mark max_altitude(.0001, .75, 1) |> floor == 452\n", "@mark max_altitude(.0001, .75, 2) |> floor == 760\n", "@mark max_altitude(.0001, .75, 3) |> floor == 997" ] }, { "cell_type": "markdown", "id": "917b20b9", "metadata": {}, "source": [ "9. Faire un plot de l'altitude atteinte par la fusée pour $Cᵈ = .75$ et les valeurs de $\\mu_0$ données ci-dessous,\n", " estimée avec $\\Delta = .01$ [s].\n", " Estimer graphiquement la valeur minimale de $\\mu_0$ permettant d'atteindre des altitudes de 1000 [m], 2000 [m] et 3000 [m].\n", "\n", " > **Indication** : il peut être utile de passer `xticks=0:1:20` en argument à la fonction `plot` pour simplifier cette estimation." ] }, { "cell_type": "code", "execution_count": null, "id": "84e1f28c", "metadata": {}, "outputs": [], "source": [ "Cᵈ, Δ = .75, .01\n", "μs = LinRange(0, 30, 200)" ] }, { "cell_type": "markdown", "id": "48521184", "metadata": {}, "source": [ "10. Pour un coefficient de trainée donné,\n", " on veut calculer la masse de carburant minimale nécessaire en vue d'atteindre une altitude $H$ donnée,\n", " en vue d'y déposer du matériel météorologique.\n", " Pour ce faire,\n", " on se propose de mettre en œuvre l'algorithme de Newton-Raphson sur une fonction scalaire prenant comme argument $\\mu_0$\n", " et s'annulant lorsque l'estimation de l'altitude obtenue par la fonction `max_altitude` est égale à $H$.\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", " qui peut être obtenue 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 la fonction `tanh(x::D)` qui est à définir." ] }, { "cell_type": "code", "execution_count": null, "id": "2d2f30ba", "metadata": {}, "outputs": [], "source": [ "import Base: +, -, *, /, ==, ≈, cos, sin, tanh, 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", "tanh(x::D) = D((tanh(x.f[1]), 1 / cosh(x.f[1])^2 * x.f[2]))\n", "@mark tanh(ε) ≈ ε\n", "@mark tanh(1000 + ε) == 1\n", "@mark tanh(-1000 + ε) == -1\n", "@mark tanh(log(2) + ε) ≈ 3/5 + 16/25*ε" ] }, { "cell_type": "markdown", "id": "c14ec12e", "metadata": {}, "source": [ "11. 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": "1dac06c0", "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 -> tanh(x) - .5, 1) ≈ atanh(.5)" ] }, { "cell_type": "markdown", "id": "40014b04", "metadata": {}, "source": [ "12. Écrire une fonction `find_fuel(H, Δ, Cᵈ, μ₀)` renvoyant la masse de carburant minimale requise pour atteindre une altitude $H$.\n", " Calculer une estimation des quantités de carburant permettant d'atteindre des altidudes de 1000 [m], 2000 [m] et 3000 [m] pour `Cᵈ = .75`,\n", " et tracer les courbes d'altitude correspondantes." ] }, { "cell_type": "code", "execution_count": null, "id": "d6c109c5", "metadata": {}, "outputs": [], "source": [ "function find_fuel(H, Δ, Cᵈ, μ₀)\n", "\n", "end\n", "\n", "Δ, Cᵈ = .01, .75\n", "# Find values of μ₀ and plot here\n", "\n", "@mark find_fuel(1000, Δ, Cᵈ, 5) |> floor == 3\n", "@mark find_fuel(2000, Δ, Cᵈ, 5) |> floor == 7\n", "@mark find_fuel(3000, Δ, Cᵈ, 5) |> floor == 12" ] } ], "metadata": { "jupytext": { "encoding": "#@ -*- coding: utf-8 -*-" }, "kernelspec": { "display_name": "Julia 1.8.5", "language": "julia", "name": "julia-1.8" } }, "nbformat": 4, "nbformat_minor": 5 }