{
"cells": [
{
"cell_type": "markdown",
"id": "0136f8ef",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" v = x₀\n",
" for k=1:maxiter\n",
" v = A*v\n",
" normalize!(v)\n",
" λ = v'A*v\n",
"\n",
" if norm(A*v - λ*v) < ε\n",
" return λ,v\n",
" end\n",
"\n",
" end\n",
" error(\"Power iteration failed to converge.\")\n",
" ### END SOLUTION\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8b28140c",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" for k=1:nev\n",
" λ,v = power_iter(A,x₀ - vs*vs'x₀;ε=ε,maxiter=maxiter)\n",
" λs[k] = λ\n",
" vs[:,k] .= v\n",
" A -= λ*v*v'\n",
" end\n",
" ### END SOLUTION\n",
"\n",
" return λs,vs\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "515eedfa",
"metadata": {
"nbgrader": {
"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": {
"lines_to_next_cell": 0,
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "legendre",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function legendre(n)\n",
" ### BEGIN SOLUTION\n",
" p = Polynomial([-1, 0, 1])\n",
" return 1 / (2^n * factorial(n)) * derivative(p^n, n)\n",
" ### END SOLUTION\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "124a4c21-3a99-4b2c-a400-84f7de3730d9",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" nodes = sort(roots(legendre(n)))\n",
" weights = zero(nodes)\n",
" for i in 1:n\n",
" ℓ = fromroots(nodes[1:end .!= i])\n",
" ℓ = ℓ / ℓ(nodes[i])\n",
" weights[i] = integrate(ℓ, -1, 1)\n",
" end\n",
" ### END SOLUTION\n",
" return nodes, weights\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e7d5b3e",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" h = (b-a)/N\n",
" X = LinRange(a, b, N + 1)\n",
" z, w = get_nodes_and_weights(n)\n",
" result = 0.\n",
" for i in 1:N\n",
" nodes = X[i] + h/2 .+ z*h/2\n",
" result += h/2 * w'u.(nodes)\n",
" end\n",
" return result\n",
" ### END SOLUTION\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad1d6e71",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"### BEGIN SOLUTION\n",
"# Number of nodes\n",
"ns = [1, 2, 3]\n",
"\n",
"# Number of cells\n",
"N = 1:40\n",
"\n",
"p = plot(title=\"Convergence of Gauss Legendre quadrature\", legend=:bottomleft,\n",
" xticks=([1, 5, 10, 20, 30], [\"1\", \"5\", \"10\", \"20\", \"30\"]))\n",
"for n in ns\n",
" errors = composite_gauss_legendre.(uᵢ, a, b, n, N) .- I_exact\n",
" polyfit = fit(log.(N), log.(abs.(errors)), 1)\n",
" β = round(- polyfit[1], digits=2)\n",
" scatter!(N, abs.(errors), label=\"n=$n, β=$β\", scale=:log10)\n",
" xlabel!(L\"N\")\n",
" ylabel!(L\"|I - \\widehat I_{n,N}|\")\n",
"end\n",
"p\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "03d1fc13",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2b3d5a8159f80773",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"1. Tracer les courbes $y_{\\alpha}(x)$ correspondant à plusieurs valeurs de $\\alpha \\in \\{1, \\dotsc, 10\\}$ pour $x\\in[-1,1]$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc762d5f",
"metadata": {
"nbgrader": {
"grade": true,
"grade_id": "plot_profiles",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"pl = plot(y(0), -1, 1, label=\"\")\n",
"for α in LinRange(1., 10., 10)\n",
" plot!(pl, y(α), label=\"\")\n",
"end\n",
"pl\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "e1908483",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "f_u_alpha",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function f(u, α)\n",
" ### BEGIN SOLUTION\n",
" x, dx = u\n",
" yp, ypp = x .|> [dy(α), d2y(α)]\n",
" return [dx, -(g*yp + yp * ypp * dx^2)/(1 + yp^2)]\n",
" ### END SOLUTION\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d077124b-d74f-47fa-a936-2da2f758a557",
"metadata": {
"nbgrader": {
"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)\n",
"### BEGIN HIDDEN TESTS\n",
"_f(x, ẋ, α) = -(g*dy(α)(x)+dy(α)(x)*d2y(α)(x)*ẋ^2)/(1+dy(α)(x)^2)\n",
"@assert all(last.(F) .≈ _f.(first.(U), last.(U), tabα))\n",
"### END HIDDEN TESTS"
]
},
{
"cell_type": "markdown",
"id": "a2c46e89",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "rk4",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function rk4(uₙ, Δ, f)\n",
" ### BEGIN SOLUTION\n",
" k₁ = f(uₙ)\n",
" k₂ = f(uₙ + Δ/2 * k₁)\n",
" k₃ = f(uₙ + Δ/2 * k₂)\n",
" k₄ = f(uₙ + Δ * k₃)\n",
" return uₙ + Δ/6 * (k₁ + 2k₂ + 2k₃ + k₄)\n",
" ### END SOLUTION\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "768a2cc1",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" while U[end][1] < 1.\n",
" # println(length(T))\n",
" push!(U, rk4(U[end], Δ, fα))\n",
" push!(T, T[end]+Δ)\n",
" end\n",
" ### END SOLUTION\n",
" return T, U\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "03ebafb9",
"metadata": {
"nbgrader": {
"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": {
"lines_to_next_cell": 2,
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" T, U = solve_ode([-1., 0.], Δ, α)\n",
" t₁, t₂, δx₁, δx₂ = T[end-1], T[end], U[end-1][1]-1, U[end][1]-1\n",
" return (t₁*δx₂ - t₂*δx₁)/(δx₂-δx₁)\n",
" ### END SOLUTION\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "85325358",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": true,
"grade_id": "plot_time",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"plot(LinRange(0.001, 20, 1000), final_time)\n",
"vline!([2.55])\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "0c4b1a03",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" for i in 1:maxiter\n",
" y = f(ForwardDiff.Dual(x, 1.))\n",
" x -= y.value/y.partials[1]\n",
" norm(f(x)) < ε && return x\n",
" end\n",
" error(\"Failed to converge!\")\n",
" ### END SOLUTION\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "776e1128-6365-4225-a50f-e12e6eaa762c",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": true,
"grade_id": "alpha_given_t",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"α₁ = newton_raphson_dual(1., x -> final_time(x) - 0.85)\n",
"α₂ = newton_raphson_dual(10., x -> final_time(x) - 0.85)\n",
"@show α₁, α₂ ;\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "f164a63e",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "alpha_opt",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"function bisection(f, a, b; δ = 1e-10)\n",
" @assert f(a) * f(b) ≤ 0\n",
" while abs(b - a) ≥ δ\n",
" x = (a + b) / 2\n",
" a, b = f(a) * f(x) ≤ 0 ? [a, x] : [x, b]\n",
" end\n",
" return (a + b) / 2\n",
"end\n",
"α_opt = bisection(x -> final_time(ForwardDiff.Dual(x, 1.)).partials[1], 1., 3.)\n",
"\n",
"# ou\n",
"\n",
"α_opt = newton_raphson_dual(1., x -> final_time(ForwardDiff.Dual(x, 1.)).partials[1])\n",
"### END SOLUTION"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "07fbdc96",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": true,
"grade_id": "animate_sol",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"function animate_sol(αs, Δ)\n",
" ### BEGIN SOLUTION\n",
" TUs = [solve_ode([-1, 0], Δ, α) for α in αs]\n",
" Ts, Us = first.(TUs), last.(TUs)\n",
" T = argmax(length, Ts)\n",
" anim = @animate for i in eachindex(T)\n",
" t = T[i]\n",
" plot(title=\"t = $(round(t, digits=3))\", xlims=(-1,1), ylims=(-0.5,1.2),\n",
" aspect_ratio=:equal, xaxis=false, yaxis=false, grid=false, ticks=false)\n",
" for j in eachindex(αs)\n",
" α = αs[j]\n",
" u = Us[j][min(i,length(Us[j]))]\n",
" plot!(y(α), aspect_ratio=:equal, label=nothing, color=:blue)\n",
" scatter!([u[1]], [y(α)(u[1])], label=nothing)\n",
" end\n",
" end\n",
" return anim\n",
" ### END SOLUTION\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
}