{
"cells": [
{
"cell_type": "markdown",
"id": "5f336a96",
"metadata": {},
"source": [
"# Cours ENPC - Pratique du calcul scientifique"
]
},
{
"cell_type": "markdown",
"id": "0bcbc4e2",
"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": "d36048ba",
"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": "0984b0f1",
"metadata": {},
"source": [
"### 1. Intégration numérique"
]
},
{
"cell_type": "markdown",
"id": "1ef3da43",
"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": "9731f8f5",
"metadata": {},
"outputs": [],
"source": [
"function composite_trapezoidal(u, a, b, n)\n",
" x = LinRange(a, b, n)\n",
" Δ = x[2] - x[1]\n",
" ux = u.(x)\n",
" return Δ * (ux[1]/2 + sum(ux[2:end-1]) + ux[end]/2)\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": "f5b51329",
"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": "f444434f",
"metadata": {},
"outputs": [],
"source": [
"function composite_simpson(u, a, b, n)\n",
" @assert n % 2 == 1 \"`n` must be impair\"\n",
" x = LinRange(a, b, n)\n",
" Δ = x[2] - x[1]\n",
" ux = u.(x)\n",
" return Δ/3 * sum([ux[1]; ux[end]; 4ux[2:2:end-1]; 2ux[3:2:end-2]])\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": "54399873",
"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": "12e110cc",
"metadata": {},
"outputs": [],
"source": [
"function calculate_sum(N)\n",
" sum(n^(-n) for n in N:-1.:1)\n",
"end\n",
"\n",
"println(calculate_sum(5))\n",
"println(calculate_sum(10))\n",
"println(calculate_sum(15))\n",
"println(calculate_sum(20))\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": "d98b951d",
"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": "de10d830",
"metadata": {},
"outputs": [],
"source": [
"ns = 3:2:400\n",
"u = x -> x^-x\n",
"I_exact = calculate_sum(20)\n",
"I_trap = composite_trapezoidal.(u, 0, 1, ns)\n",
"I_simp = composite_simpson.(u, 0, 1, ns)\n",
"plot(ns, abs.(I_trap .- I_exact), label=\"Trapezoidal\")\n",
"plot!(ns, abs.(I_simp .- I_exact), label=\"Simpson\")\n",
"plot!(xaxis=:log, yaxis=:log, lw=2)"
]
},
{
"cell_type": "markdown",
"id": "c8fe409e",
"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": "c379f11e",
"metadata": {},
"outputs": [],
"source": [
"I_simp = composite_simpson.(u, 0, 1, ns)\n",
"ns = 3:2:400\n",
"log_Δ = @. log(1/ns)\n",
"log_e = @. log(abs(I_simp - I_exact))\n",
"fit(log_Δ, log_e, 1).coeffs[2]"
]
},
{
"cell_type": "markdown",
"id": "53eb10eb",
"metadata": {},
"source": [
"### 2. Résolution d'un système linéaire"
]
},
{
"cell_type": "markdown",
"id": "3694282d",
"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": "9832e8c6",
"metadata": {},
"outputs": [],
"source": [
"function my_lu(A)\n",
" n, n = size(A)\n",
" L, U = zeros(n, n), zeros(n, n)\n",
" for i in 1:n, j in 1:i\n",
" U[j, i] = A[j, i] - L[j, 1:j-1]'U[1:j-1, i]\n",
" L[i, j] = (A[i, j] - L[i, 1:j-1]'U[1:j-1, j]) / U[j, j]\n",
" end\n",
" return L, U\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": "ca181ea0",
"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": "17f2dfb1",
"metadata": {},
"outputs": [],
"source": [
"function my_banded_lu(A, b)\n",
" n, n = size(A)\n",
" L, U = zeros(n, n), zeros(n, n)\n",
" for i in 1:n, j in max(1, i-b):i\n",
" U[j, i] = A[j, i] - L[j, max(j-b, 1):j-1]'U[max(i-b, 1):j-1, i]\n",
" L[i, j] = (A[i, j] - L[i, max(i-b, 1):j-1]'U[max(j-b, 1):j-1, j]) / U[j, j]\n",
" end\n",
" return L, U\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": "574914e9",
"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": "6119aed6",
"metadata": {},
"outputs": [],
"source": [
"function generate_banded(n, b)\n",
" A = zeros(n, n)\n",
" for i in 1:n\n",
" for j in max(1, i-b):min(n, i+b)\n",
" A[i, j] = randn()\n",
" end\n",
" end\n",
" return A\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": "0f0fda5a",
"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": "2a3b546d",
"metadata": {},
"outputs": [],
"source": [
"# Write your tests here"
]
},
{
"cell_type": "markdown",
"id": "f4776586",
"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": "c46b4d80",
"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": "c23a869b",
"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": "db8b14a3",
"metadata": {},
"outputs": [],
"source": [
"function f(t, X, Cᵈ)\n",
" z, v, μ = X\n",
" m = mᵣ + μ\n",
" return [v; Vₑ * β(μ) / m - Fᵈ(v, Cᵈ) / m - g; - β(μ)]\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": "f108a690",
"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": "1b205d2b",
"metadata": {},
"outputs": [],
"source": [
"function rkx(tₙ, Xₙ, h, Δ)\n",
" k₁ = h(tₙ, Xₙ )\n",
" k₂ = h(tₙ + Δ/2, Xₙ + Δ/2 * k₁)\n",
" k₃ = h(tₙ + 3Δ/4, Xₙ + 3Δ/4 * k₂)\n",
" return Xₙ + Δ/9 * (2k₁ + 3k₂ + 4k₃)\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": "19dffae4",
"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": "703452a1",
"metadata": {},
"outputs": [],
"source": [
"function solve_ode(Δ, Cᵈ, μ₀)\n",
" X₀ = [0.; 0.; μ₀]\n",
" ts = [0.]\n",
" Xs = [X₀]\n",
" h(t, X) = f(t, X, Cᵈ)\n",
" while Xs[end][2] ≥ 0.\n",
" push!(Xs, rkx(ts[end], Xs[end], h, Δ))\n",
" push!(ts, ts[end] + Δ)\n",
" end\n",
" return ts, Xs\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": "4676cdc5",
"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": "8c4f4a0d",
"metadata": {},
"outputs": [],
"source": [
"function plot_altitude(Δ, Cᵈ, μs)\n",
" p = plot(title=\"Altitude of the rocket\")\n",
" for μ₀ ∈ μs\n",
" ts, Xs = solve_ode(Δ, Cᵈ, μ₀)\n",
" plot!(ts, [x[1] for x in Xs], label=\"μ₀ = $μ₀\")\n",
" xlabel!(\"t [s]\")\n",
" ylabel!(\"z [m]\")\n",
" end\n",
" p\n",
"end\n",
"\n",
"Δ, Cᵈ, μs = .01, .75, [1; 2; 3; 4; 5]\n",
"plot_altitude(Δ, Cᵈ, μs)"
]
},
{
"cell_type": "markdown",
"id": "2009e4af",
"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": "87915e13",
"metadata": {},
"outputs": [],
"source": [
"function plot_velocity(Δ, Cᵈ, μs)\n",
" p = plot(title=\"Velocity of the rocket\")\n",
" for μ₀ ∈ μs\n",
" ts, Xs = solve_ode(Δ, Cᵈ, μ₀)\n",
" plot!(ts, [x[2] for x in Xs], label=\"μ₀ = $μ₀\")\n",
" xlabel!(\"t [s]\")\n",
" ylabel!(\"v [m/s]\")\n",
" end\n",
" p\n",
"end\n",
"\n",
"Δ, Cᵈ, μ₀ = .01, .75, [1; 2; 3; 4; 5]\n",
"plot_velocity(Δ, Cᵈ, μ₀)"
]
},
{
"cell_type": "markdown",
"id": "02e63691",
"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": "ef6aa79f",
"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",
" ts, Xs = solve_ode(Δ, 0, μ₀)\n",
" v = v_exact.(ts, μ₀)\n",
" v_approx = [x[2] for x in Xs]\n",
" return maximum(abs.(v_approx - v))\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": "87228159",
"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": "bc8621b1",
"metadata": {},
"outputs": [],
"source": [
"Δs, μ₀ = 2.0 .^ (-10:-2), 5\n",
"errors = error_velocity.(Δs, μ₀)\n",
"plot(Δs, errors, marker = :circle, label=L\"$L^{\\infty}$ Error\")\n",
"plot!(xaxis=:log, yaxis=:log, lw=2)"
]
},
{
"cell_type": "markdown",
"id": "7375479d",
"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": "c98a2ad5",
"metadata": {},
"outputs": [],
"source": [
"function max_altitude(Δ, Cᵈ, μ₀)\n",
" ts, Xs = solve_ode(Δ, Cᵈ, μ₀)\n",
" return max(Xs[end-1][1], Xs[end][1])\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": "7b01ef45",
"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": "dff90c65",
"metadata": {},
"outputs": [],
"source": [
"Cᵈ, Δ = .75, .01\n",
"μs = LinRange(0, 30, 200)\n",
"plot(μs, max_altitude.(Δ, Cᵈ, μs), label=\"Altitude\")\n",
"plot!(μs, zero(μs) .+ 1000, color=:red, label=\"1000 [m]\", xticks=0:1:20)\n",
"plot!(μs, zero(μs) .+ 2000, color=:red, label=\"2000 [m]\", xticks=0:1:20)\n",
"plot!(μs, zero(μs) .+ 3000, color=:red, label=\"3000 [m]\", xticks=0:1:20)\n",
"plot!(xlabel=\"\\\\mu_0 [kg]\", xlims=(0, 20))"
]
},
{
"cell_type": "markdown",
"id": "0d56604e",
"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": "65a0402f",
"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": "a146ee5c",
"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": "81f63dab",
"metadata": {},
"outputs": [],
"source": [
"function newton_raphson_dual(f, x, maxiter=100; δ = 1e-12)\n",
" for i in 1:maxiter\n",
" y = f(x + ε)\n",
" x -= y.f[1]/y.f[2]\n",
" abs(f(x)) < δ && return x\n",
" end\n",
" error(\"Failed to converge!\")\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": "8e379f06",
"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": "6d18c371",
"metadata": {},
"outputs": [],
"source": [
"function find_fuel(H, Δ, Cᵈ, μ₀)\n",
" obj = m -> max_altitude(Δ, Cᵈ, m) - H\n",
" newton_raphson_dual(obj, μ₀)\n",
"end\n",
"\n",
"Δ, Cᵈ = .01, .75\n",
"μ₁ = find_fuel(1000, Δ, Cᵈ, 5)\n",
"μ₂ = find_fuel(2000, Δ, Cᵈ, 5)\n",
"μ₃ = find_fuel(3000, Δ, Cᵈ, 5)\n",
"plot_altitude(Δ, Cᵈ, [μ₁; μ₂; μ₃])\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 release",
"language": "julia",
"name": "julia-1.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}