{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "*ENSAM-Bordeaux, Mathématiques et informatique. Date : le 13/11/19. Auteur : Éric Ducasse. Version : 1.1 ;* **compatible sympy 1.0 ** " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
*On pourra faire exécuter ce notebook cellule par cellule $\\left(\\mathtt{Maj+Entrée}\\right)$ ou intégralement par $\\mathtt{\\;Kernel\\rightarrow Restart\\;\\&\\;Run\\;All}$.*
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import sympy as sb\n", "sb.init_printing()\n", "print(\"Version de sympy :\",sb.__version__) \n", "from IPython.display import display" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# **Calcul formel : TD n°3, deuxième partie

$\\hspace{5cm}$ *Équations différentielles* ** " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
**Pour chaque exercice, faire exécuter la partie *exemples* cellule par cellule $\\left(\\mathtt{Maj+Entrée}\\right)$, avant de passer à la partie *Travail à faire*.**
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## *Exercice 15* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 15.1 Objectifs " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### * Savoir analyser un problème différentiel et en particulier prédire s'il admet des solutions formelles,
trouver sa solution générale puis déterminer les constantes qui satisfont les conditions initiales ou aux limites. *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### * Problèmes de Cauchy scalaires : $\\forall t\\;\\in]0,\\infty[$, $f^{\\prime}(t)=\\Phi(f(t),t)$ avec la condition initiale $f(0)=f_0$.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 15.2 Exemples " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $a)$ Cas dont la résolution formelle ne pose pas de problème." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On cherche la fonction $f$ telle que $$\\forall t \\in\\;]0,\\infty[\\,,\\;f^{\\prime}(t)=\\frac{1+f(t)^2}{(t+a)^2}\\;,$$avec la condition initiale $f(0)=f_0$, où $a$ est un réel strictement positif." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $a.1\\;-$ Définition de l'équation différentielle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = sb.symbols(\"t\", real=True)\n", "a = sb.symbols(\"a\", positive=True)\n", "f = sb.Function(\"f\")\n", "EDO1 = sb.Eq(f(t).diff(t),(1+f(t)**2)/(a+t)**2) ; EDO1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $a.2\\;-$ Forme générale de la solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **`sympy.dsolve` renvoie une ou plusieurs équations de la forme $f(t)=\\cdots$ :**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eq_f = sb.dsolve(EDO1,f(t)) ; eq_f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **On extrait l'expression se trouvant à droite de l'égalité par `.rhs` (*right-hand side*) :**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen = eq_f.rhs # right-hand side\n", "f_gen" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Autre écriture de la solution, que l'on se contente de vérifier :**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C1 = sb.symbols(\"C1\")\n", "f_gen2 = sb.tan( C1 - 1/(t+a) )\n", "if f_gen.equals(f_gen2) :\n", " from IPython.display import display\n", " display(f_gen2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $a.3\\;-$ Recherche de la solution en fonction de la condition initiale" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f0 = sb.symbols(\"f0\", real=True)\n", "regle_C1 = sb.solve( f_gen.replace(t,0)-f0, C1 ,dict=True)[0]\n", "regle_C1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_sol = f_gen2.xreplace(regle_C1) ; f_sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $a.4\\;-$ Tracés pour différentes valeurs initiales et différentes valeurs de $a$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On pourrait discuter des propriétés de la solution selon les valeurs de $f_0$ et de $a$... mais on ne le fait pas ici." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_num = sb.lambdify((t,f0,a), f_sol, \"numpy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = np.linspace(0,10,400)\n", "for va in (2,1.5,1,0.5) :\n", " for vf0 in (0,-1,1) :\n", " X = f_num(T,vf0,va)\n", " plt.plot(T, X, \"-\", label=r\"$f_0\\approx{:.1f}$\".format(vf0))\n", " if np.abs(X).max() > 20 : plt.ylim(-20,20)\n", " plt.title(\"a = {:.2f}\".format(va),size=14) ; plt.grid()\n", " plt.legend(loc=\"best\") ; plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $b)$ Cas plus délicat, qui montre qu'un outil de calcul formel n'est pas « *miraculeux* »" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
** *Renvoyé à la fin du fichier. À ne pas regarder dans un premier temps.* **
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 15.3 Travail à faire " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $a)$ Résoudre le problème différentiel suivant $\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall t \\in\\;]0,\\infty[ & \\displaystyle f^{\\prime}(t)=\\frac{\\exp(-2\\,t)}{1-f(t)} & , \\\\\n", "\\mathrm{avec\\;la\\;condition\\;initiale\\;:} & f(0)=s\\quad(s\\neq1) & .\\end{array}\\right.$
$\\hspace{4mm}$On distinguera les cas $s>1$ et $s<1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = sb.symbols(\"t\",positive=True)\n", "f = sb.Function(\"f\")\n", "EDO = sb.Eq( f(t).diff(t), sb.exp(-2*t)/(1-f(t)) )\n", "EDO" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **On obtient 2 solutions : **" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L_eq = sb.dsolve(EDO, f(t)) ; L_eq" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen1 = L_eq[0].rhs ; f_gen2 = L_eq[1].rhs ; (f_gen1,f_gen2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **On « voit » que le premier cas correspond à $s<1$ et le second, à $s>1$.
On recherche la constante $C_1$ en fonction de $s$ : **" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s,C1 = sb.symbols(\"s,C1\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle1 = sb.solve( sb.Eq(f_gen1.replace(t,0), s), C1, dict=True )[0] ; regle1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **La méthode `refine` (non exigible) fonctionne ici pour séparer les cas : **" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen1.replace(t,0).xreplace(regle1).refine(sb.Q.negative(s-1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen1.replace(t,0).xreplace(regle1).refine(sb.Q.positive(s-1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Solution pour $s<1$ : **" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_sol1 = f_gen1.xreplace(regle1).simplify() ; f_sol1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Recherche de la constante $C_1$ pour le deuxième solution : **" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle2 = sb.solve( sb.Eq(f_gen2.replace(t,0), s), C1, dict=True )[0] ; regle2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle2 == regle1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen2.replace(t,0).xreplace(regle2).refine(sb.Q.negative(s-1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen2.replace(t,0).xreplace(regle2).refine(sb.Q.positive(s-1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Solution pour $s>1$ : **" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_sol2 = f_gen2.xreplace(regle2).simplify() ; f_sol2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $b)$ Faire tracer les courbes représentatives des solutions sur l'intervalle $[0,2]$ pour $\n", "s\\in\\lbrace -0.1,0,0.1,0.2,0.5,1.5,1.8,1.9,2,2.1\\rbrace$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_num1 = sb.lambdify( (t,s), f_sol1, \"numpy\")\n", "f_num2 = sb.lambdify( (t,s), f_sol2, \"numpy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = np.linspace(0,2,300)\n", "plt.figure(figsize=(10,5))\n", "plt.xlim(-0.05,2.45)\n", "for vs in (-0.1,0,0.05,0.1,0.2,0.5) :\n", " plt.plot(T, f_num1(T,vs), \"-\", label=\"$s={:.2f}$\".format(vs))\n", " plt.plot(T, f_num2(T,2-vs), \"-\", label=\"$s={:.2f}$\".format(2-vs))\n", "plt.grid() ; plt.legend(loc=\"best\") ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $c)$ [Question subsidiaire] Discuter le domaine de définition des solutions en déterminant, en fonction de $s$,
$\\hspace{5mm}$les valeurs de $t$ qui rende la solution égale à $1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ecart = 1 - f_sol1 ; ecart" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ecart.equals(f_sol2-1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sb.solve(ecart,sb.exp(2*t))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle = {t:-sb.log(s*(2-s))/2} ; regle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ecart.xreplace(regle).equals(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si $0 *Exercice 16* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 16.1 Objectifs " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### * Savoir analyser un problème différentiel et en particulier prédire s'il admet des solutions formelles,
trouver sa solution générale puis déterminer les constantes qui satisfont les conditions initiales ou aux limites. *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### * Problèmes linéaires du second ordre à coefficients constants avec conditions initiales :

$\n", "\\hspace{3cm}\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall t\\;\\in]0,\\infty[\\;, & f^{\\prime\\prime}(t)+a\\,f^{\\prime}(t)+b\\,f(t)=s(t) & , \\\\\n", "\\mathrm{avec\\;les\\;conditions\\;initiales\\;:} & f(0)=f_0 \\mathrm{\\;et\\;} f^{\\prime}(0)=v_0 & .\n", "\\end{array}\\right.$*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 16.2 Exemple " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $1\\;-$ Définition de l'équation différentielle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = sb.symbols(\"t\", real=True)\n", "f_e = sb.symbols(\"f_e\", positive=True)\n", "f = sb.Function(\"f\")\n", "EDO3 = sb.Eq( f(t).diff(t,2) + 2*f(t).diff(t) + (1+4*sb.pi**2)*f(t), sb.sin(2*sb.pi*f_e*t) )\n", "EDO3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $2\\;-$ Forme générale de la solution, avec deux constantes $C_1$ et $C_2$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_gen = sb.dsolve(EDO3,f(t)).rhs\n", "f_gen" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $3\\;-$ Détermination des constantes en fonction des conditions initiales" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C1,C2 = sb.symbols(\"C1,C2\")\n", "regle_C1C2 = sb.solve( [ f_gen.replace(t,0), f_gen.diff(t).replace(t,0).simplify() ], (C1,C2))\n", "regle_C1C2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_sol = f_gen.xreplace(regle_C1C2).simplify()\n", "f_sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $4\\;-$ Définition de la fonction numérique à partir de la solution exacte, pour tracés" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = sb.lambdify( (t,f_e), f_sol, \"numpy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = np.linspace(0,3.2,500)\n", "for fe in (0.8,1,1.2) :\n", " plt.plot(T, s(T,fe), \"-\", label=\"$f_e={:.2f}$ Hz\".format(fe))\n", "plt.legend(loc=\"best\") ; plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### [Non exigible] Amplitude de la réponse en régime établi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### *Réécriture sous une forme plus élégante* (*technique non exigible*)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = sb.Wild(\"X\")\n", "denom = [ 1/T for T in f_sol.args if T.func.__name__ == \"Pow\" ][0]\n", "numer = (f_sol*denom).simplify().expand()\n", "T1 = numer.coeff(sb.exp(-t)).collect([sb.sin(X),sb.cos(X)])\n", "T0 = (numer - T1*sb.exp(-t)).expand().collect([sb.sin(X),sb.cos(X)])\n", "f_sol = (T0 + T1*sb.exp(-t))/denom ; f_sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### *Solution en régime établi, en fonction de la fréquence d'excitation $f_e$*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_sol_regime_etabli = f_sol.replace(sb.exp(X),0) ; f_sol_regime_etabli" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On reconnaît une sinusoïde de fréquence $f_e$ dont on extrait l'amplitude :" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* On extrait des coefficients C, S, D et X tels que l'expression précédente s'écrivent $\\displaystyle \n", "\\frac{C\\,\\cos(X)+S\\,\\sin(X)}{D}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C,S,D = sb.Wild(\"C\"),sb.Wild(\"S\"),sb.Wild(\"D\")\n", "regle = f_sol_regime_etabli.match( (C*sb.cos(X)+S*sb.sin(X))/D )\n", "regle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* L'amplitude de la sinusoïde est $\\displaystyle \n", "\\frac{\\sqrt{C^2+S^2}}{D}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "amplitude_regime_etabli = (sb.sqrt(C**2+S**2)/D).xreplace(regle).simplify()\n", "amplitude_regime_etabli" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Recherche des extremums :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fe_max = sb.solve(amplitude_regime_etabli.diff(f_e),f_e)[0].simplify()\n", "a_max = amplitude_regime_etabli.replace(f_e,fe_max).simplify()\n", "(fe_max,a_max)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Diagramme de Bode pour la réponse en amplitude :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = sb.lambdify( f_e, amplitude_regime_etabli, \"numpy\")\n", "Vf = np.linspace(0.2,12,400)\n", "plt.loglog(Vf, a(Vf), \"-b\", linewidth=2)\n", "plt.loglog([fe_max], [a_max], \"dr\")\n", "plt.grid() ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 16.3 Travail à faire " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $a)$ Résoudre le problème différentiel suivant : $\n", "\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall t\\;\\in]0,\\infty[\\;, & f^{\\prime\\prime}(t)+20\\,f^{\\prime}(t)+500\\,f(t)=500\\,\\theta(t-1) & , \\\\\n", "\\mathrm{avec\\;les\\;conditions\\;initiales\\;:} & f(0)=1 \\mathrm{\\;et\\;} f^{\\prime}(0)=0 & ,\n", "\\end{array}\\right.$
$\\hspace{6mm}\\theta$ désignant la fonction *échelon unitaire de Heaviside* (`sb.Heaviside`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $1\\;-$ Définition de l'équation différentielle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = sb.symbols(\"t\",positive=True)\n", "f = sb.Function(\"f\")\n", "EDO = sb.Eq(f(t).diff(t,2)+20*f(t).diff(t)+500*f(t), 500*sb.Heaviside(t-1))\n", "EDO" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $2\\;-$ Forme générale de la solution, avec deux constantes $C_1$ et $C_2$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C1,C2 = sb.symbols(\"C1,C2\")\n", "if sb.__version__ >= \"1.1\" :\n", " f_gen = sb.dsolve(EDO,f(t)).rhs.expand().collect([C1,C2,sb.Heaviside(t-1)])\n", "else : # Version 1.0 : on met le second membre à zéro\n", " f_gen = sb.dsolve(EDO.replace(sb.Heaviside,lambda X:0), f(t)).rhs\n", "f_gen" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $3\\;-$ Détermination des constantes en fonction des conditions initiales" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "CI = [ f_gen.replace(t,0)-1, f_gen.diff(t).replace(t,0) ] ; CI" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle_constantes = sb.solve(CI,[C1,C2]) ; regle_constantes" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_sol = f_gen.xreplace(regle_constantes) ; f_sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $b)$ Tracer la courbe représentative de la solution sur l'intervalle $[0,2.5]$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def my_Heaviside(v) : return (v>=0)*1 # Fonction de Heaviside vectorisée" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `sb.lambdify` aura comme troisième argument : `[{\"Heaviside\" : my_Heaviside },\"numpy\"]`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_num = sb.lambdify( t, f_sol,[{\"Heaviside\" : my_Heaviside },\"numpy\"])\n", "T = np.linspace(0,2.5,500)\n", "plt.plot(T,f_num(T)) ; plt.grid() ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## *Exercice 17* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 17.1 Objectifs " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### * Savoir analyser un problème différentiel et en particulier prédire s'il admet des solutions formelles,
trouver sa solution générale puis déterminer les constantes qui satisfont les conditions initiales ou aux limites. *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### *Problème de conduction thermique (second ordre linéaire à coefficients non-constants) avec conditions aux limites* :

$\n", "\\hspace{3cm}\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall x\\;\\in]0,e[\\;, & \\displaystyle \\frac{\\mathbb{d}}{\\mathbb{d}x}\\left(k(x)\\,\\frac{\\mathbb{d}}{\\mathbb{d}x}T(x)\\right)=-s(x) & , \\\\\n", "\\mathrm{avec\\;les\\;conditions\\;aux\\;bords\\;:} & a_0\\,T^{\\prime}(0)+b_{0}\\,T(0)=c_{0} \\mathrm{\\;et\\;} \\\\\n", "& a_1\\,T^{\\prime}(e)+b_{1}\\,T(e)=c_{1} & .\n", "\\end{array}\\right.$

$T(x)$ est un écart de température par rapport à une température de référence [K], $k(x)$ est la conductivité thermique [W/K/m] et $s(x)$ est une source volumique de chaleur [W/m$\\vphantom{m}^3$]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 17.2 Exemple " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall x\\;\\in]0,e[\\;, & \\displaystyle \\frac{\\mathbb{d}}{\\mathbb{d}x}\\left((x+e)\\,\\frac{\\mathbb{d}}{\\mathbb{d}x}T(x)\\right)=-10 & , \\\\\n", "\\mathrm{avec\\;les\\;conditions\\;aux\\;bords\\;:} & T^{\\prime}(0)=T(0) \\mathrm{\\;et\\;} T^{\\prime}(e)=-T(e)\\,/\\,2 & .\n", "\\end{array}\\right.$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $1\\;-$ Définition de l'équation différentielle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ici, la conductivité thermique $k(x)=x+e$ est une fonction affine de la position $x$ et on a une source de chaleur uniforme (constante 10) sur l'intervalle $]0,e[$. Les conditions aux bords sont des conditions de convection où la température de référence est la température extérieure et le coefficient de convection est numériquement égal à la conductivité thermique sur le bord gauche." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x,e = sb.symbols(\"x,e\", positive=True)\n", "T = sb.Function(\"T\")\n", "EDO4 = sb.Eq( ((e+x)*T(x).diff(x)).diff(x) , -10 )\n", "EDO4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $2\\;-$ Forme générale de la solution, avec deux constantes $C_1$ et $C_2$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "try :\n", " T_gen = sb.dsolve(EDO4, T(x)).rhs\n", " display(T_gen)\n", " resolue = not \"Order(\" in sb.srep(T_gen) # Solution exacte trouvée\n", "except Exception as err : # Version de sympy vraisemblablement antérieure à 1.4\n", " resolue = False\n", " print(\"Erreur :\",err)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "resolue = False" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if not resolue : # Version de sympy vraisemblablement antérieure à 1.4\n", " print(\"On résout d'abord en T'(x) :\")\n", " Tprime = sb.Function(\"T'\")\n", " EDO4_ordre1 = sb.Eq( ((e+x)*Tprime(x)).diff(x) , -10 )\n", " display( EDO4_ordre1 ) \n", " Tprime_gen = sb.dsolve(EDO4_ordre1, Tprime(x)).rhs\n", " display({Tprime(x):Tprime_gen})\n", " print(\"Puis on intègre :\")\n", " T_gen = sb.symbols(\"C2\") + Tprime_gen.integrate(x)\n", " display(T_gen)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $3\\;-$ Détermination des constantes en fonction des conditions limites" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "CLG = sb.Eq(T_gen.diff(x),T_gen).replace(x,0) ; CLG" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "CLD = sb.Eq(T_gen.diff(x),-T_gen/2).replace(x,e) ; CLD" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C1,C2 = sb.symbols(\"C1,C2\")\n", "regle_constantes = sb.solve([CLG,CLD],[C1,C2],dict=\"True\")[0]\n", "regle_constantes" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_sol = T_gen.xreplace(regle_constantes).simplify() ; T_sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Outil fourni de réécriture des logarithmes :**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simplifier_log(expr) :\n", " if isinstance(expr,list) :\n", " return [ simplifier_log(elem) for elem in expr ]\n", " elif isinstance(expr,tuple) :\n", " return tuple([ simplifier_log(elem) for elem in expr ])\n", " elif isinstance(expr,dict) :\n", " return dict([ (k,simplifier_log(v)) for k,v in expr.items() ])\n", " else :\n", " X,Y,Z = [ sb.Wild(L) for L in \"XYZ\" ]\n", " expr = expr.simplify()\n", " expr = expr.replace(sb.log(X*Y*Z),sb.log(X)+sb.log(Y)+sb.log(Z))\n", " expr = expr.replace(sb.log(X*Y),sb.log(X)+sb.log(Y))\n", " return expr.replace(sb.log(X**Y),Y*sb.log(X)).factor()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simplifier_log(regle_constantes)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "simplifier_log(T_sol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $4\\;-$ Fonction numérique représentant la solution pour $e=0.2\\,\\mathrm{m}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_num = sb.lambdify( (x,e), T_sol, \"numpy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_ref = lambda x,T_num=T_num : T_num(x,0.2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "e_num = 0.2\n", "Vx = np.linspace(0,e_num,200)\n", "plt.plot( Vx, T_ref(Vx), \"-m\" )\n", "plt.xlabel( \"Position $x$ [m]\", size=14)\n", "plt.ylabel( \"$\\delta T(x)$ [K]\", color=\"m\", size=14)\n", "plt.grid() ;\n", "plt.twinx()\n", "plt.plot( Vx, e_num+Vx, \"-b\" ) \n", "plt.ylabel( \"$k(x)$ [W/K/m]\", color=\"b\", size=14);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 17.3 Travail à faire " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Pour $e=\\frac{1}{5}$, on veut résoudre le même problème que dans l'exemple :
$\\hspace{3cm}\n", "\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall x\\;\\in]0,e[\\;, & \\displaystyle \\frac{\\mathbb{d}}{\\mathbb{d}x}\\left((x+e)\\,\\frac{\\mathbb{d}}{\\mathbb{d}x}T(x)\\right)=-s(x) & , \\\\\n", "\\mathrm{avec\\;les\\;conditions\\;aux\\;bords\\;:} & T^{\\prime}(0)=T(0) \\mathrm{\\;et\\;} T^{\\prime}(e)=-T(e)\\,/\\,2 & .\n", "\\end{array}\\right.$
sauf que la source est localisée et uniforme sur l'intervalle $\\left[\\frac{e}{2},\\frac{3\\,e}{4}\\right]$ et vaut 40 W/m$\\vphantom{m}^3$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $a)$ Déterminer la forme générale de la solution pour les deux équations différentielles suivantes :
$\\hspace{1cm}\n", "\\displaystyle \\frac{\\mathbb{d}}{\\mathbb{d}x}\\left((x+e)\\,\\frac{\\mathbb{d}}{\\mathbb{d}x}T(x)\\right)=0\\quad$ et $\\quad\\displaystyle \\frac{\\mathbb{d}}{\\mathbb{d}x}\\left((x+e)\\,\\frac{\\mathbb{d}}{\\mathbb{d}x}T(x)\\right)=-40$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = sb.symbols(\"x\", positive=True)\n", "e = sb.S(1)/5\n", "T = sb.Function(\"T\")\n", "EDO0 = sb.Eq( ((e+x)*T(x).diff(x)).diff(x) , 0 )\n", "EDO0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "try :\n", " T_gen0 = sb.dsolve(EDO0,T(x)).rhs\n", " display(T_gen0)\n", " resolue = not \"Order(\" in sb.srep(T_ge0) # Solution exacte trouvée\n", "except Exception as err :\n", " resolue = False\n", " print(\"Erreur :\",err)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if not resolue : # Version de sympy vraisemblablement antérieure à 1.4\n", " print(\"On résout d'abord en T'(x) :\")\n", " Tprime = sb.Function(\"T'\")\n", " EDO0_ordre1 = sb.Eq( ((e+x)*Tprime(x)).diff(x) , 0 )\n", " display( EDO0_ordre1 ) \n", " Tprime_gen0 = sb.dsolve(EDO0_ordre1, Tprime(x)).rhs\n", " display({Tprime(x):Tprime_gen0})\n", " print(\"Puis on intègre :\")\n", " T_gen0 = sb.symbols(\"C2\") + Tprime_gen0.integrate(x)\n", " display(T_gen0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "EDOS = sb.Eq( ((e+x)*T(x).diff(x)).diff(x) , -40 )\n", "EDOS" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "try :\n", " T_genS = sb.dsolve(EDOS,T(x)).rhs\n", " display(T_genS)\n", " resolue = not \"Order(\" in sb.srep(T_genS) # Solution exacte trouvée\n", "except Exception as err :\n", " resolue = False\n", " print(\"Erreur :\",err)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if not resolue : # Version de sympy vraisemblablement antérieure à 1.4\n", " print(\"On résout d'abord en T'(x) :\")\n", " Tprime = sb.Function(\"T'\")\n", " EDOS_ordre1 = sb.Eq( ((e+x)*Tprime(x)).diff(x) , -40 )\n", " display( EDOS_ordre1 ) \n", " Tprime_genS = sb.dsolve(EDOS_ordre1, Tprime(x)).rhs\n", " display({Tprime(x):Tprime_genS})\n", " print(\"Puis on intègre :\")\n", " T_genS = sb.symbols(\"C2\") + Tprime_genS.integrate(x)\n", " display(T_genS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $b)$ Sur chacun des intervalles $\\left]0,\\frac{e}{2}\\right[$, $\\left]\\frac{e}{2},\\frac{3\\,e}{4}\\right[\n", "$ et $\\left]\\frac{3\\,e}{4},e\\right[$, définir la forme générale de la solution en remplaçant les constantes $C_1$ et $C_2$ par d'autres portant des noms différents." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C1,C2,a1,a2,a3,b1,b2,b3 = sb.symbols(\"C1,C2,a1,a2,a3,b1,b2,b3\")\n", "T_gen1 = T_gen0.xreplace({C1:a1,C2:b1})\n", "T_gen2 = T_genS.xreplace({C1:a2,C2:b2})\n", "T_gen3 = T_gen0.xreplace({C1:a3,C2:b3})\n", "(T_gen1,T_gen2,T_gen3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $c)$ En déduire un système (linéaire) de 6 équations à 6 inconnues,
$\\hspace{5mm}$construit en écrivant les conditions aux limites en $x=0$ et en $x=e$,
$\\hspace{5mm}$ainsi que les conditions de raccordement
$\\hspace{10mm}$en $x=\\frac{e}{2}$ : « $T\\!\\left(\\frac{e}{2}^-\\right)=T\\!\\left(\\frac{e}{2}^+\\right)$ » et « $T^{\\prime}\\!\\left(\\frac{e}{2}^-\\right)=T^{\\prime}\\!\\left(\\frac{e}{2}^+\\right)$ », et
$\\hspace{10mm}$en $x=\\frac{3\\,e}{4}$ : « $T\\!\\left(\\frac{3\\,e}{4}^-\\right)=T\\!\\left(\\frac{3\\,e}{4}^+\\right)$ » et « $T^{\\prime}\\!\\left(\\frac{3\\,e}{4}^-\\right)=T^{\\prime}\\!\\left(\\frac{3\\,e}{4}^+\\right)$ »." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "LEQ = [sb.Eq(T_gen1.diff(x),T_gen1).replace(x,0),\\\n", " sb.Eq(T_gen3.diff(x),-T_gen3/2).replace(x,e),\\\n", " sb.Eq(T_gen1,T_gen2).replace(x,e/2),\\\n", " sb.Eq(T_gen1.diff(x),T_gen2.diff(x)).replace(x,e/2),\\\n", " sb.Eq(T_gen2,T_gen3).replace(x,3*e/4),\\\n", " sb.Eq(T_gen2.diff(x),T_gen3.diff(x)).replace(x,3*e/4)\\\n", " ] ; LEQ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $d)$ Le résoudre, puis utiliser la fonction `simplifier_log` fournie pour récrire le résultat." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle_ctes = simplifier_log(sb.solve(LEQ,(a1,a2,a3,b1,b2,b3))) ; regle_ctes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $e)$ Définir la solution du problème différentiel sur chaque intervalle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_sol1 = simplifier_log(T_gen1.xreplace(regle_ctes)) ; T_sol1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_sol2 = simplifier_log(T_gen2.xreplace(regle_ctes)) ; T_sol2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_sol3 = simplifier_log(T_gen3.xreplace(regle_ctes)) ; T_sol3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $f)$ Définir les fonctions numériques correspondantes, puis la fonction numérique vectorisée représentant la solution globale." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T_num1 = sb.lambdify( x, T_sol1, \"numpy\")\n", "T_num2 = sb.lambdify( x, T_sol2, \"numpy\")\n", "T_num3 = sb.lambdify( x, T_sol3, \"numpy\")\n", "def T_num(x,ve=float(e)) :\n", " V = T_num1(np.clip(x,0,0.5*ve))*(x<=0.5*ve)\n", " V += T_num2(np.clip(x,0.5*ve,0.75*ve))*(x>0.5*ve)*(x<0.75*ve)\n", " V += T_num3(np.clip(x,0.75*ve,ve))*(x>=0.75*ve)\n", " return V" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $g)$ Tracer la courbe de la solution et la superposer à celle de l'exemple." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "e_num = float(e)\n", "Vx = np.linspace(0,e_num,200)\n", "plt.plot( Vx, T_num(Vx), \"-r\", label=\"solution\" )\n", "plt.xlabel( \"Position $x$ [m]\")\n", "plt.ylabel( \"$\\delta T(x)$\")\n", "#plt.plot( Vx, T_ref(Vx), \"-b\" , label=\"source uniforme\")\n", "plt.grid() ; plt.legend() ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## *Exercice 15 (suite)* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $b)$ Cas plus délicat, qui montre qu'un outil de calcul formel n'est pas « *miraculeux* »" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On cherche la fonction $f$ telle que $\\displaystyle\\left\\lbrace\\begin{array}{lll}\n", "\\forall t \\in\\;]0,\\infty[ & f^{\\prime}(t)=\\cos^2(f(t)) & , \\\\\n", "\\mathrm{avec\\;la\\;condition\\;initiale\\;:} & f(0)=f_0 & .\\end{array}\\right.$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $b.1\\;-$ Définition de l'équation différentielle" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = sb.symbols(\"t\", real=True)\n", "f = sb.Function(\"f\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* L'équation différentielle est ici du premier ordre et non-linéaire " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "EDO2 = sb.Eq(f(t).diff(t),sb.cos(f(t))**2)\n", "EDO2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $b.2\\;-$ Forme(s) générale(s) de la solution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eqs_f = sb.dsolve(EDO2,f(t)) ; eqs_f" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "verif = [ EDO2.replace(f(t),eq.rhs) for eq in eqs_f ]\n", "from IPython.display import display\n", "for eq in verif : display(eq)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[ eq.doit().simplify() for eq in verif ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Alors qu'on connaît la solution générale :
$\\hspace{6mm}$si $\\cos(f(t))\\neq 0$, $\\displaystyle\\frac{f^{\\prime}(t)}{\\cos^2(f(t))}=\n", "\\partial_t\\left(\\tan(f(t))\\right)=1\\:\\Longleftrightarrow\\;\\tan(f(t))=t+a$ ;
$\\hspace{6mm}$Sinon, si $\\cos(f_0)\\neq 0$, alors $f$ est la fonction constante $f_0$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sb.tan(f(t)).diff(t).simplify()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def solution_exacte(f0,t=t) :\n", " if (f0%sb.pi).equals(sb.pi/2) :\n", " return f0\n", " else :\n", " a = sb.tan(f0)\n", " return (f0-sb.atan(a)) + sb.atan(t+a)\n", "valeurs_initiales = [0,1,sb.pi/2,3*sb.pi/4,3*sb.pi/2,11*sb.pi/6]\n", "exemples = [ solution_exacte(v) for v in valeurs_initiales ]\n", "exemples" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[ sb.simplify(EDO2.replace(f(t),F).doit()) for F in exemples ]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[ F.replace(t,0) for F in exemples ] == valeurs_initiales" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $b.3\\;-$ Recherche de la solution en fonction de la condition initiale" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "C1,f0 = sb.symbols(\"C1,f0\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "valeurs_initiales = [ eq.rhs.replace(t,0).simplify() for eq in eqs_f ] ; valeurs_initiales" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol_C1 = [ sb.solve(v-f0, C1, dict=True) for v in valeurs_initiales ]\n", "sol_C1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F1,F2 = [ eq.rhs.xreplace(s_C1[0]) for eq,s_C1 in zip(eqs_f,sol_C1) ]\n", "F1,F2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### $b.4\\;-$ Comparaisons graphiques pour différentes valeurs initiales" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F1num,F2num = [ sb.lambdify((t,f0), F, \"numpy\") for F in (F1,F2) ]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "val_init = [0,-1,sb.pi/2,2,9] ; val_init" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = np.linspace(0,3,100)\n", "for fnum0 in val_init :\n", " if (fnum0%sb.pi).equals(sb.pi/2) :\n", " Fexacte = lambda t : np.ones_like(t)*float(fnum0)\n", " plt.plot([], [], \"r-\", label=\"$F_1(t)$ non définie\")\n", " plt.plot([], [], \"b-\", label=\"$F_2(t)$ non définie\")\n", " else :\n", " a = np.tan(fnum0)\n", " Fexacte = lambda t,a=a,f0=float(fnum0) : (f0-np.arctan(a)) + np.arctan(t+a)\n", " plt.plot(T, F1num(T,fnum0), \"r-\", label=\"$F_1(t)$\")\n", " plt.plot(T, F2num(T,fnum0), \"b-\", label=\"$F_2(t)$\")\n", " plt.plot(T, Fexacte(T), \"g:\", label=\"$F_e(t)$\", linewidth=3)\n", " plt.legend() ; plt.grid() ; \n", " plt.title(r\"$f_0\\;\\approx\\;${:.3f}\".format(float(fnum0))) ; plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }