{
"cells": [
{
"cell_type": "markdown",
"id": "657ee43e-abd8-4efe-b5e0-6a632114529a",
"metadata": {
"tags": []
},
"source": [
"# Cours ENPC - Pratique du calcul scientifique"
]
},
{
"cell_type": "markdown",
"id": "a30eaeba",
"metadata": {},
"source": [
"## Résolution numérique d'équations différentielles ordinaires"
]
},
{
"cell_type": "markdown",
"id": "d3623bef",
"metadata": {},
"source": [
"--------------"
]
},
{
"cell_type": "markdown",
"id": "a97e8668",
"metadata": {},
"source": [
"### Problème de la couche limite de Blasius\n",
"\n",
"Le problème de la vitesse d'un fluide visqueux newtonien s'écoulant parallèlement à une plaque plane, dans un régime laminaire à nombre de Reynolds élevé, se résout en raccordant une solution triviale (écoulement uniforme) en champ lointain et une solution de couche limite au voisinage de la plaque pour respecter la condition de vitesse nulle. Sans rentrer dans les détails relevant du cours de mécanique des fluides, on rappelle que ce problème de couche limite après adimensionnement, se ramène à la résolution de l'équation différentielle suivante\n",
"$$\n",
"2g'''+g\\,g''=0\n",
"\\quad\\textrm{et}\\quad\n",
"\\left\\{\n",
"\\begin{array}{rcl}\n",
"g(0)&=&0\\\\\n",
"g'(0)&=&0\\\\\n",
"\\lim_{t\\to\\infty}g'(t)&=&1\n",
"\\end{array}\n",
"\\right.\n",
"$$\n",
"\n",
"A des normalisations près, $t$ représente ici l'abscisse perpendiculaire à la plaque et $u=g'$ est la vitesse du fluide dans la direction de la plaque.\n",
"\n",
"Il apparaît de prime abord que ce problème ne relève pas *stricto sensu* des résultats vus en cours car d'une part il s'agit d'un problème impliquant des dérivées strictement supérieures à 1 et d'autre part les conditions aux limites devant être respectées par la solution impliquent à la fois les valeurs de $g$ et $g'$ à $t=0$ mais également de $g'$ à $t\\to\\infty$.\n",
"\n",
"La stratégie mise en place ici pour résoudre ce problème est le recours à la méthode dite **méthode de tir** : elle consiste à se ramener à une équation différentielle d'ordre 1 vectorielle en dimension 3 et à omettre provisoirement la condition $\\lim_{t\\to\\infty}g'(t)=1$ pour la remplacer par une nouvelle condition en $t=0$ à savoir $g''(0)=λ$ avec $λ$ un paramètre qu'il s'agira ensuite d'ajuster pour satisfaire la condition à l'infini.\n",
"\n",
"1. Montrer que l'on se ramène à une équation différentielle ordinaire d'ordre 1 de type\n",
" \n",
" $$\n",
" \\tag{1}\n",
" x'(t)=f(x(t)) \\quad ; \\quad x(0)=\\begin{pmatrix} 0 \\\\ 0 \\\\ λ \\end{pmatrix}\n",
" $$\n",
" en considérant la fonction vectorielle inconnue de $\\mathbb{R}$ dans $\\mathbb{R}^3$\n",
" $$\n",
" \\begin{array}{rccl}\n",
" x :&\n",
" [0,\\infty[&→&\\mathbb{R}^3 \\\\\n",
" &t&↦&x(t)= \\begin{pmatrix}\n",
" g(t)\\\\\n",
" g'(t)\\\\\n",
" g''(t)\n",
" \\end{pmatrix}\n",
" \\end{array}\n",
" $$\n",
" et la fonction de $\\mathbb{R}^3$ dans $\\mathbb{R}^3$\n",
" $$\n",
" \\begin{array}{rccl}\n",
" f :&\n",
" \\mathbb{R}^3&→&\\mathbb{R}^3 \\\\\n",
" &x&↦&f(x)= \\begin{pmatrix}\n",
" x_2\\\\\n",
" x_3\\\\\n",
" -\\frac{x_1\\,x_3}{2}\n",
" \\end{pmatrix}\n",
" \\end{array}\n",
" $$\n",
"\n",
" Construire les fonctions `f_Blasius` et `df_Blasius` prenant comme argument un vecteur de dimension 3 `x` et renvoyant respectivement le vecteur $f(x)$ et la matrice jacobienne $∇f(x)$.\n",
"\n",
"1. Implémenter une fonction `newton_raphson(x, f, df, maxiter=100; ε = 1e-12)` **dans le cadre général d'une fonction vectorielle** renvoyant le résultat de l'algorithme de Newton-Raphson partant d'un point initial `x` pour une fonction `f` de jacobienne `df` 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 $\\lVert f(x) \\rVert<ε$.\n",
"\n",
"1. On donne ci-dessous les codes renvoyant l'itération suivante des schémas d'Euler explicite et implicite à partir de la valeur précédente `xₙ`, l'incrément `Δ` de la variable $t$ ainsi que la fonction `f` décrivant l'équation différentielle (1) et éventuellement la jacobienne `df` de `f` si elle est nécessaire dans le schéma (on notera que si celle-ci n'est pas nécessaire on peut remplacer l'argument `df` par `_` pour éviter de le nommer mais il est important de garder ce 4ème argument pour respecter la même signature pour tous les schémas)\n",
"\n",
" ```julia\n",
" Euler_exp(xₙ, Δ, f, _) = xₙ+Δ*f(xₙ)\n",
" Euler_imp(xₙ, Δ, f, df) = newton_raphson(xₙ, x -> x-xₙ-Δ*f(x), x -> I-Δ*df(x))\n",
" ```\n",
"\n",
" Après avoir recopié ces schémas, implémenter de manière analogue les schémas de `Crank_Nicolson` et de `Heun`\n",
"\n",
"1. Implémenter un solveur générique `solve_ode(x₀, tᶠ, Δ ; alg=Euler_exp)` prenant comme arguments\n",
"\n",
" - le vecteur initial `x₀`,\n",
"\n",
" - la valeur finale `tᶠ` de l'intervalle de résolution (ne pouvant bien sûr pas prendre l'infini dans une résolution numérique on supposera dans la suite que la valeur de $10$ sera suffisante pour représenter une valeur grande en notant que l'échelle de grandeur de $t$ fournie par une analyse dimensionnelle de l'équation différentielle initiale est l'unité),\n",
"\n",
" - le pas de résolution `Δ`,\n",
"\n",
" - `alg` le choix du schéma (parmi les fonctions implémentées à la question précédente ou d'autres que vous voudriez tester...).\n",
"\n",
" Ce solveur générique devra renvoyer les tableaux des valeurs de $t$ (vecteur de scalaires) et de $x$ (vecteur de vecteurs de $\\mathbb{R}^3$).\n",
"\n",
"1. Tester la résolution en traçant pour les différents algorithmes $t$ en fonction de $g'(t)=x_2(t)$ (il suffit d'inverser les axes pour que $t$ représente l'axe des ordonnées et $g'(t)$ la vitesse du fluide parallèle à l'axe des $x$). On choisira $x_0=(0, 0, 0.2)$ puis $x_0=(0, 0, 0.5)$. Conjecturer l'existence d'un $λ$ permettant de satisfaire la condition à l'infini.\n",
"\n",
" Pour trouver la valeur adéquate de $λ$, on se propose de mettre en œuvre l'algorithme de Newton-Raphson sur une fonction scalaire prenant comme argument $λ$ et qui s'annule lorsque l'estimation de la valeur à l'infini de $g'$ est respectée (autrement dit $x_2(t^f)-1=0$). Cette fonction impose donc la résolution numérique complète de l'équation différentielle puisqu'elle fait intervenir $x_2(t^f)$. Il n'est donc pas question de pouvoir la différentier par rapport à $λ$. C'est pourquoi on se propose de s'appuyer sur la notion de différentiation automatique vue au TD précédent. On donne ci-dessous la structure `D` de nombre dual avec les fonctions nécessaires et suffisantes pour entreprendre la résolution numérique de l'équation de Blasius à l'aide de nombres duaux.\n",
"\n",
"1. Ecrire une fonction de résolution par Newton-Raphson `newton_raphson_dual(x, f, maxiter=100; ε = 1e-12)` d'une fonction scalaire `f` dans laquelle la dérivée de `f` au point courant sera obtenue par exploitation des nombres duaux.\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)))`.\n",
"\n",
"1. Implémenter la fonction `shooting_residual(λ, tᶠ, Δ; alg=Crank_Nicolson)` devant s'annuler lorsque la résolution respecte la condition estimée \"à l'infini\".\n",
"\n",
"1. Tester la résolution de la valeur de $λ$ en appliquant l'algorithme de Newton-Raphson à la fonction `shooting_residual`.\n",
"*Indication :* attention la fonction `newton_raphson_dual` est implémentée avec une fonction `f` ne dépendant que d'un seul argument donc il faut se ramener à une fonction déduite de `shooting_residual` qui ne dépend plus que du seul argument `λ`. Il suffit pour cela de fixer localement les valeurs de autres arguments `tᶠ`, `Δ` et `alg` et d'appeler `newton_raphson_dual` avec la fonction anonyme `λ -> shooting_residual(λ, tᶠ, Δ; alg=alg)`.\n",
"\n",
"1. Tracer les courbes donnant les valeurs de $λ$ en fonction de $Δ$ pour les différents schémas.\n",
"\n",
"1. En déduire une bonne estimation de $λ$ et tracer à nouveau le profil de vitesse (i.e. $t$ en fonction de $g'(t)=x_2(t)$) pour cette valeur de $λ$ pour les différents schémas.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0bd98a07",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra, Plots"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "adbfd89b",
"metadata": {},
"outputs": [],
"source": [
"f_Blasius(x) = # votre code ici\n",
"df_Blasius(x) = # votre code ici"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "afb40085",
"metadata": {},
"outputs": [],
"source": [
"function newton_raphson(x, f, df, maxiter=100; ε = 1e-12)\n",
" # votre code ici doit retourner la valeur convergée de la racine de f\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "24dc10f0",
"metadata": {},
"outputs": [],
"source": [
"Euler_exp(xₙ, Δ, f, _) = xₙ+Δ*f(xₙ)\n",
"Euler_imp(xₙ, Δ, f, df) = newton_raphson(xₙ, x -> x-xₙ-Δ*f(x), x -> I-Δ*df(x))\n",
"Crank_Nicolson(xₙ, Δ, f, df) = # votre code ici\n",
"Heun(xₙ, Δ, f, _) = # votre code ici"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7c2eaf1d",
"metadata": {},
"outputs": [],
"source": [
"function solve_ode(x₀, tᶠ, Δ ; alg=Crank_Nicolson)\n",
" X = [x₀]\n",
" T = [0.]\n",
" # votre code ici\n",
" return T, X\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ed64a293",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"tᶠ = 10.\n",
"Δ = tᶠ/101\n",
"\n",
"pl=plot(xlabel=\"u\", ylabel=\"y\")\n",
"for λ in (0.2,0.5)\n",
" x₀ = [0.,0.,λ]\n",
" for alg in (Euler_imp, Euler_exp, Crank_Nicolson, Heun)\n",
" # votre code ici pour tracer les profils de vitesse u=g' afin de tester le solveur\n",
" end\n",
"end\n",
"pl"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9006115e",
"metadata": {},
"outputs": [],
"source": [
"# Définition de la structure de nombre dual\n",
"# Ne pas modifier mais exécuter pour pouvoir l'utiliser\n",
"\n",
"import Base: +, -, *, /, inv, conj, abs, isless, AbstractFloat, convert, promote_rule\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",
"abs(x::D) = D((abs(x.f[1]), x.f[2]*sign(x.f[1])))\n",
"conj(x::D) = D(conj.(x.f))\n",
"isless(x::D, y::D) = isless(x.f[1],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{Float64}, x::D) = x.f[1]\n",
"promote_rule(::Type{D}, ::Type{<:Real}) = D\n",
"Base.show(io::IO,x::D) = print(io,x.f[1],x.f[2]<0 ? \" - \" : \" + \",abs(x.f[2]),\" ε\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b47a3759",
"metadata": {},
"outputs": [],
"source": [
"function newton_raphson_dual(x, f, maxiter=100; ε = 1e-12)\n",
" # votre code ici d'implémentation de Newton-Raphson exploitant les nombres duaux pour estimer la dérivée de f\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5763d974",
"metadata": {},
"outputs": [],
"source": [
"function shooting_residual(λ, tᶠ, Δ; alg=Crank_Nicolson)\n",
" # votre code ici devant retourner un scalaire s'annulant lorsque la condition finale est satisfaite\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6a3a9701",
"metadata": {},
"outputs": [],
"source": [
"lΔ = tᶠ ./ (2 .^(4:13))\n",
"pl=plot(xlabel=\"Δ\", ylabel=\"u′(0)\", xaxis=:log2)\n",
"for alg in (Euler_imp, Euler_exp, Crank_Nicolson, Heun)\n",
" # votre code ici pour tester le solveur\n",
"end\n",
"pl"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b86e211d",
"metadata": {},
"outputs": [],
"source": [
"λ = # votre code ici pour estimer la bonne valeur de λ"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6790a16e",
"metadata": {},
"outputs": [],
"source": [
"tᶠ = 10.\n",
"Δ = tᶠ/101\n",
"\n",
"pl=plot(xlabel=\"u\", ylabel=\"y\")\n",
"x₀ = [0.,0.,λ]\n",
"for alg in (Euler_imp, Euler_exp, Crank_Nicolson, Heun)\n",
" # votre code ici pour tracer le profil de vitesse avec la bonne valeur de λ\n",
"end\n",
"pl"
]
},
{
"cell_type": "markdown",
"id": "42d828bf",
"metadata": {},
"source": [
"### Solution d'un problème aux valeurs initiales d'ordre trois.\n",
"\n",
"\n",
"$$\n",
"\\tag{1}\n",
"\\begin{align*}\n",
" &v''' + v'' + 4v' + 4v = 4t^2 + 8t - 10, \\\\\n",
" &v(0) = -3, \\quad v'(0) = -2, \\quad v''(0) = 2.\n",
"\\end{align*}\n",
"$$\n",
"\n",
"#### Question 1 ####\n",
"\n",
"1. Vérifier que $v_\\mathrm{e}(t) = -\\sin(2t)+t^2-3$ est la solution exacte du problème (1).\n",
"\n",
"1. Réécrire le problème comme un système d'équations du premier ordre de la forme\n",
"\n",
" $$ \\mathbf{u}' = \\mathbf{F}(t,\\mathbf{u}), \\qquad \\mathbf{u}(t) \\colon \\mathbb{R}^+ \\cup \\{0\\} \\mapsto \\mathbb{R}^3, \\qquad \\mathbf{F}(t,\\mathbf{u}) \\colon \\mathbb{R}^+ \\cup \\{0\\} \\times \\mathbb{R}^3 \\mapsto \\mathbb{R}^3. $$\n",
"\n",
" Spécifier aussi la condition initiale comme un 3-vecteur.\n",
"\n",
"1. Implémenter la méthode d'Euler explicite et la méthode de Runge-Kutta d'ordre 4 (RK4) en `Julia`.\n",
" Cette dernière méthode est basée sur l'itération suivante:\n",
" $$\n",
" \\mathbf u_{n+1} = \\mathbf u_n + \\frac{h}{6}\\left(\\mathbf k_1 + 2\\mathbf k_2 + 2\\mathbf k_3 + \\mathbf k_4 \\right),\n",
" $$\n",
" où\n",
" \\begin{align*}\n",
" \\mathbf k_1 &= \\ f(t_n, \\mathbf u_n), \\\\\n",
" \\mathbf k_2 &= \\ f\\!\\left(t_n + \\frac{h}{2}, \\mathbf u_n + \\frac{h}{2} \\mathbf k_1\\right), \\\\\n",
" \\mathbf k_3 &= \\ f\\!\\left(t_n + \\frac{h}{2}, \\mathbf u_n + \\frac{h}{2} \\mathbf k_2\\right), \\\\\n",
" \\mathbf k_4 &= \\ f\\!\\left(t_n + h, \\mathbf u_n + h\\mathbf k_3\\right).\n",
" \\end{align*}\n",
" Les méthodes doivent prendre en arguments les entrées suivantes : la fonction `f`, le temps initial `ti`, le temps final `tf`, le pas de temps `h`, et la valeur initiale `u0`.\n",
" Elles doivent renvoyer la solution `uout` et le vecteur de temps correspondant `tout`.\n",
"\n",
"1. Appliquez les deux méthodes au système obtenu dans la partie 2.\n",
" Tracez sur le même graphique la solution exacte et deux solutions approximatives (une pour chaque méthode).\n",
"\n",
"1. Effectuer un test de convergence pour conclure que la méthode d'Euler est précise au premier ordre, alors que la méthode RK4 est précise au quatrième ordre.\n",
" Pour ce faire, diviser $h$ par deux, c'est-à-dire $h \\to h/2$, 8 fois.\n",
" Sur une échelle logarithmique, tracer en fonction de la taille du pas de temps l'erreur $\\|\\mathrm{err}\\|_\\infty$ (la norme $\\infty$ de l'erreur), définie comme suit\n",
" $$ \\|\\mathrm{err}\\|_\\infty = \\max_{0\\leq n\\leq t_f/h} |v_\\mathrm{e}(t_n) - V_n|, $$\n",
" où $V_n$ est la solution approchée au temps $t_n$ et où on prendra $t_f = 2$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c76adf1",
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"using LaTeXStrings\n",
"using LinearAlgebra"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5747b7f6",
"metadata": {},
"outputs": [],
"source": [
"# Solution to question 3\n",
"function euler_fwd(f, ti, tf, h, u0)\n",
" # Fill me\n",
"end\n",
"\n",
"function rk4(f, ti, tf, h, u0)\n",
" # Fill me\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "380e2487",
"metadata": {},
"outputs": [],
"source": [
"# Solution to question 4"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d84f341e",
"metadata": {},
"outputs": [],
"source": [
"# Solution to question 5"
]
}
],
"metadata": {
"citation-manager": {
"items": {}
},
"jupytext": {
"encoding": "# -*- coding: utf-8 -*-"
},
"kernelspec": {
"display_name": "Julia 1.9.0",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}