{
"cells": [
{
"cell_type": "markdown",
"id": "d7de331f",
"metadata": {},
"source": [
"# Cours ENPC - Pratique du calcul scientifique"
]
},
{
"cell_type": "markdown",
"id": "8bf20f40",
"metadata": {},
"source": [
"Before you turn in this assignment, make sure everything runs as expected. First, **restart the kernel** and then **run all cells**. Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\", as well as your name and group number below:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "79e538d9",
"metadata": {},
"outputs": [],
"source": [
"NAME = \"\"\n",
"GROUP = \"\""
]
},
{
"cell_type": "markdown",
"id": "0136f8ef",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "5d715064cf2c5af13022a609b86f1d64",
"grade": false,
"grade_id": "cell-9546f1bdb893cdcd",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Final exam\n",
"\n",
"- Ce notebook est à soumettre sur Educnet avant 11h45.\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 bibliothèques utilisées dans ce notebook. Si une ou plusieurs d'entre elles manquent sur votre machine, vous êtes invités à les installer au préalable dans le gestionnaire de bibliothèques d'une console."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a57b1a8f",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "be157bcca3a7f8baa4d2e61f8454feb8",
"grade": false,
"grade_id": "cell-4d638a879ba6f86e",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"using ForwardDiff\n",
"using LaTeXStrings\n",
"using LinearAlgebra\n",
"using Plots\n",
"using Polynomials\n",
"using Random\n",
"\n",
"Plots.default(titlefontsize=14,\n",
" xlabelfontsize=12,\n",
" ylabelfontsize=12,\n",
" legendfontsize=12,\n",
" xtickfontsize=12,\n",
" ytickfontsize=12)"
]
},
{
"cell_type": "markdown",
"id": "fe25b04b",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "76a2826e12cd247361c44db258355aaa",
"grade": false,
"grade_id": "cell-ae1c522a2bcd7f2a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## [Exercice 1] Une méthode de déflation pour calculer plusieurs valeurs propres dominantes"
]
},
{
"cell_type": "markdown",
"id": "a9d93a90",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "830e46d98ded01cf0f9a655c4681b435",
"grade": false,
"grade_id": "cell-25b349cbda605e93",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Le but de cet exercice est de mettre en œuvre une méthode pour approximer plusieurs valeurs propres dominantes d'une matrice hermitienne ${\\sf A} \\in \\mathbb{C}^{n\\times n}$.\n",
"Pour décrire la méthode,\n",
"notons $ |\\lambda_1|\\geq |\\lambda_2|\\geq \\dots \\geq |\\lambda_n|$ les valeurs propres de $\\sf A$,\n",
"et définissons la suite de matrices :\n",
"$$ {\\sf A}_1 := {\\sf A},\\qquad {\\sf A}_{k+1} = {\\sf A}_{k} - \\lambda_{k} v_{k} v_{k}^*,$$\n",
"où $v_{k}$ est un vecteur propre normalisé de $\\mathsf A$ associé à la valeur propre $\\lambda_k$.\n",
"Il est facile de vérifier par récurrence que $(\\lambda_k,v_k)$ est le couple propre dominant de ${\\sf A}_k$.\n",
"La méthode de déflation consiste simplement à définir la suite\n",
"$$ \\widetilde{\\sf A}_1 := {\\sf A},\\quad \\widetilde{\\sf A}_{k+1} = \\widetilde{\\sf A}_{k} - \\widetilde\\lambda_{k} \\widetilde v_{k} \\widetilde v_{k}^*,$$\n",
"où $(\\widetilde\\lambda_{k},\\widetilde v_{k})$ est une approximation numérique du couple propre dominant de $ \\widetilde{\\sf A}_{k} $ calculée par l'itération de la puissance.\n",
"\n",
"1. Implémenter une fonction `power_iter(A, x₀; ε=1e-8, maxiter=100000)` prenant comme arguments une matrice `A` de taille $n\\times n$,\n",
" un vecteur initial de $n$ éléments `x₀`, un seuil de tolérance `ɛ` et un nombre maximal d'itérations `maxiter`,\n",
" et retournant un tuple `(λ, x)` contenant le résultat de l'itération de la puissance appliquée à `A`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "94188e2a",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "15b6aa9d4f53b7e1535592b66d9b8ad9",
"grade": false,
"grade_id": "power_iter",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function power_iter(A, x₀; ε=1e-8, maxiter=100000)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8b28140c",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "540b842d0a2ec7b62c4eed146b4f5dbb",
"grade": true,
"grade_id": "power_iter_tests",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"@assert power_iter([1. 2.; 2. 1.], [1., 0.])[1] ≈ 3.\n",
"@assert power_iter([1. 0.; 0. .5], [1., 1.])[1] ≈ 1.\n",
"@assert [1, -1]'power_iter([1. 2.; 2. 1.], [1., 0.])[2] |> abs < 1e-6\n",
"@assert [1, 0]'power_iter([0. 0.; 0. 1.], [1., 1.])[2] |> abs < 1e-6\n",
"@assert [0, 1]'power_iter([1. 0.; 0. .5], [1., 1.])[2] |> abs < 1e-6"
]
},
{
"cell_type": "markdown",
"id": "ec8cdd87",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "918633d1c058f993a979f73371ca9ec5",
"grade": false,
"grade_id": "cell-c9618e0c7e895d1d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"2. Implémenter une fonction `deflation_method(A, x₀, nev; ε=1e-8, maxiter=100000)` prenant comme arguments la matrice `A`, un vecteur initial `x₀`, le nombre de valeurs propres désirées `nev`, la tolérance `ɛ` et le nombre d'itérations `maxiter` pour chaque appel à l'itération de la puissance.\n",
" La fonction devra retourner un tuple `(λs, vs)` où `λs` est un vecteur de `nev` valeurs propres approximatives,\n",
" classées de la plus grande à la plus petite en valeur absolue,\n",
" et `vs` est une matrice dont les colonnes sont constituées des vecteurs propres associés."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b9087f01",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "1be5362158ce6bebb352ac44b548ffc3",
"grade": false,
"grade_id": "deflation_method",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function deflation_method(A, x₀, nev; ε=1e-8, maxiter=100000)\n",
" n = size(A, 1)\n",
"\n",
" (nev > n) && error(\"$nev eigenvalues required for a $n×$n matrix\")\n",
"\n",
" λs = zeros(eltype(A), nev)\n",
" vs = zeros(eltype(A), n, nev)\n",
"\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"\n",
" return λs,vs\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "515eedfa",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "fe64ac9bd69e32f982ac46f020cabf46",
"grade": true,
"grade_id": "deflation_method_tests",
"locked": true,
"points": 2,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"# Tests automatiques\n",
"N = 200\n",
"seed = 2024\n",
"A = randn(Xoshiro(seed), N, N)\n",
"A = (A+A')/sqrt(N)\n",
"x₀ = randn(N)\n",
"\n",
"nev = 10\n",
"\n",
"@time λs, us = deflation_method(A,x₀,nev)\n",
"@assert all(λs .≈ sort(eigvals(A), by=abs, rev=true)[1:nev])"
]
},
{
"cell_type": "markdown",
"id": "d0540971",
"metadata": {
"deletable": false,
"editable": false,
"lines_to_next_cell": 0,
"nbgrader": {
"cell_type": "markdown",
"checksum": "da67a9144df9393c017629d42b917aa5",
"grade": false,
"grade_id": "cell-2b7096f389287011",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### [Exercice 2] Intégration composite de Gauss-Legendre\n",
"\n",
"1. Écrire une fonction `legendre(n)` qui retourne le polynôme de Legendre de degré $n$,\n",
" sous forme d'une structure `Polynomial` de la bibliothèque `Polynomials`.\n",
" Pour ce faire, vous pouvez utiliser la bibliothèque `Polynomials` et la formule de Rodrigues :\n",
" $$\n",
" L_n(x) = \\frac{1}{2^n n!} \\frac{\\mathrm{d}^n}{\\mathrm{d} x^n} \\left(x^2 - 1\\right)^n.\n",
" $$\n",
"\n",
" \n",
" \n",
" Indication (cliquer pour afficher)\n",
"
\n",
"\n",
" - La fonction `factorial(n)` permet de calculer la factiorielle de `n`.\n",
"\n",
" - La fonction `Polynomials.Polynomial` permet de créer un polynôme à partir de ses coefficients :\n",
" ```julia\n",
" p = Polynomial([1, 2, 3]) # p(x) = 1 + 2x + 3x²\n",
" ```\n",
" - La fonction `Polynomials.derivative` permet de calculer les dérivées d'un polynôme :\n",
" ```julia\n",
" dp = derivative(p) # dp(x) = 2 + 6x\n",
" ddp = derivative(p, 2) # ddp(x) = 6\n",
" ```\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "40c5b77b",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "903c8b6f23a488dcae3feb31be2d65c5",
"grade": false,
"grade_id": "legendre",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function legendre(n)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "124a4c21-3a99-4b2c-a400-84f7de3730d9",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "f6a32e444046df48f24ee946deab3ed5",
"grade": true,
"grade_id": "legendre_tests",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"X = Polynomial([0, 1])\n",
"for n in 1:5\n",
" @assert (n+1)*legendre(n+1) == (2n+1)*X*legendre(n)-n*legendre(n-1)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "4d1d4efb",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "e418daa42cd874ba9dc47a827efe9edd",
"grade": false,
"grade_id": "cell-84fc8003305ec45c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"2. Écrire une fonction `get_nodes_and_weights(n)` qui calcule,\n",
" sans utiliser d'autres bibliothèques logicielles que celles importées au début du notebook,\n",
" les nœuds $(x_i)_{i \\in \\{1, \\dots, n\\}}$ et poids $(w_i)_{i \\in \\{1, \\dots, n\\}}$ de la quadrature de Gauss-Legendre avec $n$ nœuds.\n",
" Pour rappel, les nœuds et poids doivent être tels que l'approximation\n",
" $$\n",
" \\int_{-1}^{1} f(x) \\, \\mathrm d x\n",
" \\approx \\sum_{i=1}^{n} w_i f(x_i)\n",
" $$\n",
" soit exacte pour tout polynôme $f$ de degré au plus $2n-1$.\n",
" \n",
" \n",
" Indication (cliquer pour afficher)\n",
"
\n",
"\n",
" - On rappelle que les nœuds d'intégration sont donnés par les racines du polynôme de Legendre de degré `n`.\n",
" Ces racines peuvent être calculées par la fonction `roots` de la biblothèque `Polynomials.jl`.\n",
"\n",
" - Pour construire les polynômes de Lagrange en vue de calculer les poids,\n",
" il peut être utile d'utiliser les fonctions `fromroots` et `integrate` de la biblothèque `Polynomials.jl`.\n",
"\n",
" ```julia\n",
" p = fromroots([1., 2.]) # Constructs (x - 1)(x - 2) = x² - 3x + 2\n",
" q = integrate(p) # q = x^3/3 - 3x^2/2 + 2x\n",
" ```\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6499aa1",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "9763785292b76c26d29e8a22212e1501",
"grade": false,
"grade_id": "get_nodes_and_weights",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function get_nodes_and_weights(n)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
" return nodes, weights\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e7d5b3e",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "1d1b51c4267f42c14895d2ef161b6e6f",
"grade": true,
"grade_id": "get_nodes_and_weights_tests",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"@assert get_nodes_and_weights(5) |> length == 2\n",
"@assert get_nodes_and_weights(5)[1] |> length == 5\n",
"@assert get_nodes_and_weights(5)[2] |> length == 5\n",
"@assert get_nodes_and_weights(1)[1] ≈ [0.]\n",
"@assert get_nodes_and_weights(1)[2] ≈ [2.0]\n",
"@assert get_nodes_and_weights(3)[1] .|> legendre(3) |> abs∘sum < 1e-10\n",
"@assert get_nodes_and_weights(5)[1] .|> legendre(5) |> abs∘sum < 1e-10\n",
"@assert get_nodes_and_weights(5)[2] |> sum ≈ 2"
]
},
{
"cell_type": "markdown",
"id": "780f44f1",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "1ce1d116cd4d114d7929202d4be26d53",
"grade": false,
"grade_id": "cell-7d8cb5ba773ea661",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"3. Écrire une fonction `composite_gauss_legendre(u, a, b, n, N)` qui renvoie une approximation de l'intégrale\n",
" $$\n",
" \\int_{a}^{b} u(x) \\, \\mathrm{d} x\n",
" $$\n",
" obtenue en partitionnant l'intervalle d'intégration $[a, b]$ en $N$ cellules de même taille,\n",
" et en appliquant la quadrature de Gauss-Legendre avec $n$ nœuds dans chaque cellule."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0e642215",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "f38b06e3ba4ff79a0556b69ccf451b09",
"grade": false,
"grade_id": "composite_gauss_legendre",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function composite_gauss_legendre(u, a, b, n, N)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad1d6e71",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "77e7dd11883ea7b4e3d5aceebdf7fe06",
"grade": true,
"grade_id": "composite_gauss_legendre_tests",
"locked": true,
"points": 2,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"_short(f, n, N) = composite_gauss_legendre(f, 0, 1, n, N)\n",
"for d in 1:9\n",
" @assert _short(x -> x^d, 5, 1) ≈ 1/(d+1)\n",
" @assert _short(x -> x^d, 5, 2) ≈ 1/(d+1)\n",
" @assert _short(x -> x^d, 5, 3) ≈ 1/(d+1)\n",
"end\n",
"@assert !(_short(x -> x^10, 2, 1) ≈ 1/11)\n",
"@assert !(_short(x -> x^10, 2, 2) ≈ 1/11)\n",
"@assert _short(x -> x^10, 5, 200) ≈ 1/11\n",
"@assert _short(x -> exp(x), 5, 200) ≈ ℯ - 1"
]
},
{
"cell_type": "markdown",
"id": "fe3fae17",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "f6eec1826356ec7b431ca4061db21a9b",
"grade": false,
"grade_id": "cell-8f08122a9c4bcdad",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"4. Considérons le cas particulier où $u(x) = \\cos(x)$, $a = -1$ et $b = 1$,\n",
" et définissons l'erreur d'intégration,\n",
" vue comme une fonction de $N$ où $n$ est un paramètre fixé,\n",
" par la formule\n",
" $$\n",
" E_{n}(N) = \\lvert \\widehat I_{n, N} - I_{\\rm exact} \\rvert.\n",
" $$\n",
" Dans cette équation,\n",
" $I_{\\rm exact}$ est la valeur exacte de l'intégrale\n",
" tandis que $\\widehat I_{n, N}$ est son approximation par la règle de Gauss-Legendre composite.\n",
" Il est demandé\n",
"\n",
" - d'estimer, pour chaque valeur de $n \\in \\{1, 2, 3\\}$,\n",
" l'ordre de convergence de la quadrature de Gauss-Legendre composite par rapport à $N$,\n",
" c'est-à-dire de trouver $\\beta = \\beta(n)$ tel que\n",
" $$\n",
" E_n(N) \\propto C N^{-\\beta}.\n",
" $$\n",
"\n",
" - d'illustrer sur un même graphique,\n",
" à l'aide de la fonction `Plots.scatter`,\n",
" les fonctions $E_1, E_2, E_3$,\n",
" pour des valeurs de $N$ variant de 1 à 40.\n",
" Utiliser l'échelle logarithmique pour les deux axes,\n",
" et inclure l'ordre de convergence `β` trouvé au point précédent dans la légende,\n",
" en passant par exemple à la fonction `scatter` l'argument `label=\"n=$n, β=$β\"`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "caa0f074",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "9a595624a443bf9c7c1457927440dee6",
"grade": true,
"grade_id": "plot_legendre_error",
"locked": false,
"points": 2,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Function to integrate\n",
"uᵢ(x) = cos(x)\n",
"\n",
"# Integration interval\n",
"a, b = -1, 1\n",
"\n",
"# Exact value of the integral\n",
"I_exact = 2sin(1)\n",
"\n",
"# YOUR CODE HERE\n",
"throw(ErrorException(\"No code provided\"))"
]
},
{
"cell_type": "markdown",
"id": "03d1fc13",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "b0a640f2b4088d9c10ccc9a3b2f21455",
"grade": false,
"grade_id": "cell-4a19561ec8ac613f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### [Exercice 3] Trajectoire d'une masse ponctuelle sur un rail"
]
},
{
"cell_type": "markdown",
"id": "0cfe0e47",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "5815253dc8c7408ffbeaf55031b7be0c",
"grade": false,
"grade_id": "cell-8a547a50ccf6374b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"On considère dans cet exercice un solide de masse $m$ modélisé comme un point susceptible de se déplacer sans frottement le long d'un arc paramétré $\\mathbf{M}(x)=\\left(x, y(x)\\right)^T$ de classe $\\mathcal{C}^2$ dans le plan $xy$. On introduit les vecteurs tangent $\\mathbf{t}(x)$ et normal $\\mathbf{n}(x)$ à l'arc par\n",
"\n",
"$$\n",
"\\mathbf{t}(x)=\\frac{\\mathbf{M}'(x)}{\\lVert \\mathbf{M}'(x) \\rVert}=\\frac{(1, y'(x))^T}{\\lVert \\mathbf{M}'(x) \\rVert}\n",
"\\quad ; \\quad\n",
"\\mathbf{n}(x)=\\frac{(-y'(x), 1)^T}{\\lVert \\mathbf{M}'(x) \\rVert}\n",
"\\quad \\textrm{avec} \\quad\n",
"\\lVert \\mathbf{M}'(x) \\rVert=\\sqrt{1+y'(x)^2}\n",
"$$\n",
"\n",
"Par convention on notera dans la suite $y'$ la dérivation de $y$ par rapport à $x$ et par $\\dot{y}$ la dérivation par rapport au temps, ce qui entraîne que $\\dot{y}=y' \\dot{x}$.\n",
"\n",
"Le solide est soumis à la gravité $\\mathbf{P}=-mg\\mathbf{e}_y$ et à la réaction du support $\\mathbf{R}=R_n\\mathbf{n}$ (pas de frottement). On montre alors par projection sur $\\mathbf{t}(x)$ du principe fondamental de la dynamique que le mouvement est régi par la résolution d'une équation différentielle d'ordre 2 sur $x(t)$\n",
"\n",
"$$\n",
"\\ddot{x}=\\frac{-g y'(x) - y'(x) y''(x) \\dot{x}^2}{1+y'(x)^2}\n",
"\\quad ; \\quad\n",
"x(0)=x_0\n",
"\\quad ; \\quad\n",
"\\dot{x}(0)=\\dot{x}_0\n",
"\\tag{EDOx}\n",
"$$\n",
"\n",
"Cette équation différentielle se ramène classiquement à une équation vectorielle d'ordre 1 d'inconnue $\\mathbf{u}=(x, \\dot{x})^T$ sous la forme suivante\n",
"\n",
"$$\n",
"\\dot{\\mathbf{u}}=\\mathbf{f}(\\mathbf{u})=\n",
"\\begin{pmatrix}\n",
"u_2\\\\\n",
"\\frac{-g y'(u_1) - y'(u_1) y''(u_1) u_2^2}{1+y'(u_1)^2}\n",
"\\end{pmatrix}\n",
"\\quad ; \\quad\n",
"\\mathbf{u}(0)=\\begin{pmatrix} x_0 \\\\ \\dot{x}_0 \\end{pmatrix}\n",
"\\tag{EDOu}\n",
"$$\n",
"On considère que la fonction $y(x)$ est donnée par\n",
"$$\n",
"y_{\\alpha}(x) = \\frac{e^{-\\alpha x} - e^{-\\alpha}}{e^{\\alpha} - e^{-\\alpha}},\n",
"$$\n",
"où $α$ est un paramètre strictement positif.\n",
"Noter que $y_{\\alpha}(-1) = 1$ et $y_{\\alpha}(1) = 0$.\n",
"La fonction $y_\\alpha$ et ses dérivées sont implémentées ci-dessous sous la forme de fonctions de fonctions, i.e. `y(α)` renvoie la fonction $y_\\alpha$ de $x$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e2971944",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "d0bdaa1162b6f2932542d36144177eb6",
"grade": false,
"grade_id": "cell-e9d4527fa7ace911",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"const g = 9.81 ;\n",
"const m = 0.03 ;\n",
"\n",
"y(α) = x -> (exp(-α * x) - exp(-α)) / (exp(α) - exp(-α))\n",
"dy(α) = x -> -α*exp(-α * x) / (exp(α) - exp(-α))\n",
"d2y(α) = x -> α^2*exp(-α * x) / (exp(α) - exp(-α))"
]
},
{
"cell_type": "markdown",
"id": "7f2cb2bc",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "24f7e5a9587fef6afbce965ecf99e33f",
"grade": false,
"grade_id": "cell-2b3d5a8159f80773",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"1. Tracer les courbes correspondant à plusieurs valeurs de $\\alpha \\in \\{1, \\dotsc, 10\\}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc762d5f",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "17dbcf8f1dbd5c4eac5a98c84bb483fc",
"grade": true,
"grade_id": "plot_profiles",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"throw(ErrorException(\"No code provided\"))"
]
},
{
"cell_type": "markdown",
"id": "e1908483",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "13201b6fa8613056ff332dbc6b73ec42",
"grade": false,
"grade_id": "cell-a8301fcbf0a4ebe9",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"2. Écrire la fonction $f$ sous forme d'une fonction Julia `f(u, α)`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d76a2cbc",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "46f2ca0c25cb250142a350ff14ca3b37",
"grade": false,
"grade_id": "f_u_alpha",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function f(u, α)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d077124b-d74f-47fa-a936-2da2f758a557",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "7bc26277d5e7583ded43340407ff0bbe",
"grade": true,
"grade_id": "f_u_alpha_tests",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"U, tabα = [2rand(2) .- 1 for _ in 1:10], rand(10)\n",
"F = f.(U,tabα)\n",
"@assert first.(F) == last.(U)"
]
},
{
"cell_type": "markdown",
"id": "a2c46e89",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "469714abfe65b731e86681eb6a16adfd",
"grade": false,
"grade_id": "cell-ef18dcbc116b067d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"3. Écrire une fonction `rk4(uₙ, 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 $u' = f(u)$.\n",
" Cette méthode est basée sur l'itération suivante:\n",
" $$\n",
" \\mathbf u_{n+1} = \\mathbf u_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 &= \\ \\mathbf f(\\mathbf u_n), \\\\\n",
" \\mathbf k_2 &= \\ \\mathbf f\\!\\left(\\mathbf u_n + \\frac{\\Delta}{2} \\mathbf k_1\\right), \\\\\n",
" \\mathbf k_3 &= \\ \\mathbf f\\!\\left(\\mathbf u_n + \\frac{\\Delta}{2} \\mathbf k_2\\right), \\\\\n",
" \\mathbf k_4 &= \\ \\mathbf f\\!\\left(\\mathbf u_n + \\Delta\\mathbf k_3\\right).\n",
" \\end{align*}\n",
" La fonction devra renvoyer $\\mathbf u_{n+1}$. À noter que l'argument générique `f` est considéré ici comme une fonction d'une unique variable vectorielle `u`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6ee52678",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "f1cfc1764d051cc526569c3da646c869",
"grade": false,
"grade_id": "rk4",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function rk4(uₙ, Δ, f)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "768a2cc1",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "f3aeb2f64d49653195a82e287ea54bea",
"grade": true,
"grade_id": "rk4_tests",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"@assert rk4([0.], 1., u -> [1.]) ≈ [1.0]\n",
"@assert rk4([3.,12.], 1., u -> u .* [2.,3.]) ≈ [21.,196.5]"
]
},
{
"cell_type": "markdown",
"id": "ebc05bd0",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "609ddb485c5915811c18968442790874",
"grade": false,
"grade_id": "cell-1d9afc657548a3b7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"4. Écrire une fonction `solve_ode(u₀, Δ, α)` pour une condition initiale `u₀` et un paramètre `α`,\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 `T` et un vecteur de vecteurs `U` contenant la solution à ces temps.\n",
" On demande d'interrompre l'intégration numérique dès que la valeur de la coordonnée $x$ sera devenue supérieure ou égale à 1;\n",
" il faudra donc que `U[end-1][1]` soit inférieur à 1 et `U[end][1]` soit supérieur ou égal à 1."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ced29f41",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "df38383efef212af405b1c70725a1609",
"grade": true,
"grade_id": "solve_ode",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function solve_ode(u₀, Δ, α)\n",
" fα(u) = f(u, α)\n",
" U = [u₀ .+ 0.0 * α] # pour éviter les problèmes de type\n",
" T = [0.]\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
" return T, U\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "03ebafb9",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "1e497bd38d7db32d254f7798acd2984b",
"grade": false,
"grade_id": "cell-9551dbbc0e31ebb3",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"5. Écrire une fonction `final_time(α)`,\n",
" qui retourne une approximation du temps mis par le solide pour atteindre la position $x = 1$.\n",
" Pour ce faire, résoudre l'équation différentielle avec un pas de temps $Δ = 0.01$,\n",
" et estimer le temps requis par interpolation linéaire sur le dernier pas de temps,\n",
" durant lequel la coordonnée $x$ du solide passe au delà de 1. Les conditions initiales sont $(x_0,\\dot{x}_0)=(-1,0)$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ee3c76f3",
"metadata": {
"deletable": false,
"lines_to_next_cell": 2,
"nbgrader": {
"cell_type": "code",
"checksum": "4e15c2516fdb027b4d05f02398722d4c",
"grade": true,
"grade_id": "final_time",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function final_time(α)\n",
" Δ = .01\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "85325358",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "facf41ca557255dbc273eb6e39f53ee7",
"grade": false,
"grade_id": "cell-4bf19352e5b24af2",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"6. Faire un plot du temps final en fonction de `α`.\n",
" Estimer graphiquement, par exemple à l'aide de la fonction `vline`,\n",
" la valeur de `α` permettant de minimiser le temps final."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3029100b",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "c6d239851e4e6bbb5eaed166d18f0609",
"grade": true,
"grade_id": "plot_time",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"throw(ErrorException(\"No code provided\"))"
]
},
{
"cell_type": "markdown",
"id": "0c4b1a03",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "899b814ebbb2b7b4cf69b32125f5c9be",
"grade": false,
"grade_id": "cell-1abeb30d548bbb1b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"7. On se propose ici de déterminer la ou les valeurs de $α$ permettant d'atteindre le point final en un temps donné par une méthode de Newton-Raphson s'appuyant sur les nombres duaux (afin d'être en mesure de dériver une fonction de $\\alpha$). Écrire une fonction `newton_raphson_dual(x, f, maxiter=100; ε = 1e-12)` renvoyant une racine de la fonction `f` en partant d'un point initial `x`.\n",
"\n",
" \n",
" \n",
" Indication (cliquer pour afficher)\n",
"
\n",
"\n",
" On rappelle que l'on peut obtenir simultanément la valeur et la dérivée d'une fonction `f` en `x` par la méthode suivante\n",
"\n",
" ```julia\n",
" y = f(ForwardDiff.Dual(x, 1.))\n",
" fx = y.value # renvoie f(x)\n",
" dfx = y.partials[1] # renvoie f'(x)\n",
" ```\n",
"\n",
" **Pour éviter des problèmes de confusion en manipulant plusieurs ordres de dérivation, il est préférable d'utiliser `f(ForwardDiff.Dual(x, 1.))` plutôt que `f(x+ForwardDiff.Dual(0., 1.))`.**\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "46c5fda6",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "17d60fd50544f06f17ffe128a66b1edc",
"grade": false,
"grade_id": "newton_raphson_dual",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function newton_raphson_dual(x, f, maxiter=100; ε = 1e-12)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "776e1128-6365-4225-a50f-e12e6eaa762c",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "dd6e2c2733b79b27c791b85313b429a1",
"grade": true,
"grade_id": "newton_raphson_dual_tests",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"@assert newton_raphson_dual(1, x -> x^2 - 2) ≈ √2\n",
"@assert newton_raphson_dual(-1, x -> x^2 - 2) ≈ -√2\n",
"@assert newton_raphson_dual(1, x -> x^3 - 2) ≈ cbrt(2)\n",
"@assert newton_raphson_dual(2, x -> cos(x) - .5) ≈ acos(.5)"
]
},
{
"cell_type": "markdown",
"id": "23893e1a",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "fe89d70f98459d0c87313df9deada204",
"grade": false,
"grade_id": "cell-e7c0348b3375d803",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"8. Déterminer la ou les valeurs de $\\alpha$ permettant d'atteindre le point final en $t=0.85s$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7f3a8447",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "70323f0bdf328ca5da1722f078515b07",
"grade": true,
"grade_id": "alpha_given_t",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"throw(ErrorException(\"No code provided\"))"
]
},
{
"cell_type": "markdown",
"id": "f164a63e",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "8738af76f9885f4a0bfe0a7ac2e141c3",
"grade": false,
"grade_id": "cell-12abfd3820e22d2a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"9. Ayant constaté graphiquement l'existence d'un paramètre `α` optimal,\n",
" calculer précisément ce paramètre que l'on nommera `α_opt` en utilisant la méthode de votre choix."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9397a469",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "c20c064c34723672ee7d0ccdd18d3032",
"grade": false,
"grade_id": "alpha_opt",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"throw(ErrorException(\"No code provided\"))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "07fbdc96",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "3f777b30541b1cf2da7b8346e0a58cb9",
"grade": true,
"grade_id": "alpha_opt_tests",
"locked": true,
"points": 2,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"@assert abs(α_opt-2.677)<1e-3"
]
},
{
"cell_type": "markdown",
"id": "a35c0048",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "8d89602ce8c2372f92dbbddf72eee9fe",
"grade": false,
"grade_id": "cell-04419576bd3ebddc",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"10. Écrire une fonction d'animation `animate_sol(αs, Δ)` permettant de calculer des solutions correspondant à plusieurs valeurs de $\\alpha$ dans `αs` pour un même pas de temps `Δ`,\n",
" et de superposer les trajectoires partant de conditions initiales sont $(x_0,\\dot{x}_0)=(-1,0)$.\n",
"\n",
" Bien prêter attention au fait que les vecteurs temps récupérés lors des différentes simulations ne sont pas de même longueur mais ont le même pas entre éléments consécutifs.\n",
"\n",
" En notant `α_opt` la valeur définie à la question précédente, on pourra appliquer l'animation à la liste `[α_opt/100, α_opt/2, α_opt, 2α_opt, 4α_opt]` et constater que c'est bien la trajectoire liée à `α_opt` qui arrive en premier."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "080f84f4",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "09cab64df15b1b1bab419a9705d54e99",
"grade": true,
"grade_id": "animate_sol",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function animate_sol(αs, Δ)\n",
" # YOUR CODE HERE\n",
" throw(ErrorException(\"No code provided\"))\n",
"end\n",
"\n",
"αs = [α_opt/100, α_opt/2, α_opt, 2α_opt, 4α_opt]\n",
"gif(animate_sol(αs, 0.01), fps=10, show_msg=false)"
]
}
],
"metadata": {
"citation-manager": {
"items": {}
},
"jupytext": {
"encoding": "# -*- coding: utf-8 -*-"
},
"kernelspec": {
"display_name": "Julia 1.10.4",
"language": "julia",
"name": "julia-1.10"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}