{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "*ENSAM-Bordeaux, Mathématiques et informatique. Date : le 04/11/19. Auteur : Éric Ducasse. Version : 1.1*" ] }, { "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" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sympy as sb\n", "sb.init_printing()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f\"Version de sympy : {sb.__version__}\") " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# **Calcul formel : TD n°2, seconde partie** " ] }, { "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 12* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 12.1 Objectifs " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### * Savoir définir et résoudre un système d'équations algébriques, en utilisant la fonction $\\mathtt{sympy.solve}$;
savoir utiliser le résultat d’un calcul symbolique pour créer une fonction, en maîtrisant la fonction $\\mathtt{sympy.lambdify}$. *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Documentation sur les équations et les solveurs :
https://docs.sympy.org/latest/tutorial/solvers.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Documentation sur $\\mathtt{lambdify}$ :
https://docs.sympy.org/latest/tutorial/basic_operations.html#lambdify" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 12.2 Exemples " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problème à résoudre : déterminer les paramètres $a$ et $b$ d'une fonction paramétrée pour satisfaire des conditions initiales données." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*$a)$ Définition d'une expression représentant une fonction de la variable $t$ avec 4 paramètres $a$, $b$, $\\omega$ et $\\tau$.*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a,b,t = sb.symbols(\"a,b,t\", real=True)\n", "w,tau = sb.symbols(\"omega,tau\", positive=True)\n", "U_de_t = (a*sb.cos(w*t)+b*sb.sin(w*t))*sb.exp(-t/tau)\n", "U_de_t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*$b)$ Valeurs initiales de la fonction et de sa dérivée.*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "valeur_initiale = U_de_t.replace(t,0) ; valeur_initiale" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "derivee_initiale = U_de_t.diff(t).replace(t,0) ; derivee_initiale" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*$c)$ Détermination des paramètres $a$ et $b$ qui permettent à la fonction de satisfaire des conditions initiales données.*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u0,v0 = sb.symbols(\"u_0,v_0\")\n", "equations = [sb.Eq(valeur_initiale,u0), sb.Eq(derivee_initiale,v0)]\n", "equations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle_ab = sb.solve(equations, [a,b]) ; regle_ab" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Variante, sans utiliser d'équation :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s_annulent = [valeur_initiale-u0, derivee_initiale-v0]\n", "s_annulent" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sb.solve(s_annulent, [a,b])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*$d)$ Solution du problème.*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "solu = U_de_t.xreplace(regle_ab).simplify()\n", "solu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*$e)$ Définition de fonctions utilisant le dernier résultat.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction $\\mathtt{lambdify}$ est faite pour définir une fonction à valeurs numériques :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(sb.lambdify.__doc__[:487])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Fonction à valeurs symboliques : on utilise $\\mathtt{xreplace}$**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "alpha,f = sb.symbols(\"alpha,f\", positive=True)\n", "solu.xreplace({tau:1/alpha, w:2*sb.pi*f}).simplify()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Us(time,initial_value=1,initial_derivative=0,frequency=1,time_constant=1,expression=solu) : \n", " regle = {u0:initial_value, v0:initial_derivative, w:sb.pi*2*frequency, tau:time_constant}\n", " return expression.xreplace(regle).simplify()\n", "# Remarque : on peut aussi définir une fonction lambda, avec des performances identiques\n", "Us(t)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Us(t,frequency=f,time_constant=1/alpha,initial_value=0,initial_derivative=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Fonction vectorisée à valeurs numériques**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Avec $\\mathtt{lambdify}$, on ne peut pas définir de valeurs par défaut.\n", "Si on veut le faire, il faut procéder en deux temps, comme ci-dessous." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u_numerique = sb.lambdify( (t,u0,v0,w,tau), solu, \"numpy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = np.arange(0,1.01,0.2)\n", "print(np.array( [T,u_numerique(T,1,0,2*np.pi,1)] ))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Un = lambda t,u0=1,v0=0,w=2*np.pi,tau=1 : u_numerique(t,u0,v0,w,tau)\n", "print(np.array( [T,Un(T)] ))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Vt = np.linspace(0,3.2,500)\n", "plt.plot(Vt,Un(Vt),\"-r\",linewidth=1.6,label=r\"$f$ = 1 Hz ; $\\tau$ = 1 s\")\n", "for fq,ct in ((1,2),(2,1)) :\n", " plt.plot(Vt,Un(Vt,w=2*np.pi*fq,tau=ct),\"-\",linewidth=1.6,\\\n", " label=r\"$f$ = {} Hz ; $\\tau$ = {} s\".format(fq,ct))\n", "plt.grid() ; plt.legend() ; plt.xlabel(\"Temps $t$ [s]\") ; plt.ylabel(\"$u(t)$\") ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 12.3 Travail à faire " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'objet de cet exercice est de définir 4 fonctions $p$ polynomiales de degré 3, que l'on considère uniquement sur l'intervalle $[0,1]$ telles que l'une des valeurs parmi $[p(0),p^{\\prime}(0),p(1),p^{\\prime}(1)]$ est 1 et les trois autres sont 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$a)$ Soit l'expression $P=p(X)=a+b\\,X+c\\,X^2+d\\,X^3$. Déterminer la liste $[p(0),p^{\\prime}(0),p(1),p^{\\prime}(1)]$ en fonction des paramètres $a$, $b$, $c$ et $d$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a,b,c,d,X = sb.symbols(\"a,b,c,d,X\")\n", "P = ( a + X*(b + X*(c + X*d))).expand(); P" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "valeurs_limites = [P.replace(X,0), P.diff(X).replace(X,0), P.replace(X,1), P.diff(X).replace(X,1)]\n", "valeurs_limites" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$b)$ En déduire les expressions correspondants aux 4 fonctions cherchées." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.display import display\n", "expr = []\n", "for p0,q0,p1,q1 in np.eye(4, dtype=np.int) :\n", " print(f\"+++ p(0)={p0}, p'(0)={q0}, p(1)={p1}, p'(1)={q1} +++\") \n", " equations = [ sb.Eq(p,v) for p,v in zip((p0,q0,p1,q1),valeurs_limites) ]\n", " display(equations)\n", " regle_abcd = sb.solve(equations,(a,b,c,d))\n", " display(regle_abcd)\n", " expr.append( P.xreplace(regle_abcd).simplify() )\n", " display(expr[-1])\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$c)$ Convertir ces expressions en 4 fonctions numériques et faire tracer leurs courbes représentatives sur l'intervalle $[0,1]$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "LF = [ sb.lambdify(X,E,\"numpy\") for E in expr ]\n", "LF" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Vx = np.linspace(0,1,200)\n", "for n,f in enumerate(LF,1) :\n", " plt.plot(Vx,f(Vx),\"-\",linewidth=2.0, label=f\"$p_{n}(x)$\")\n", "plt.grid() ; plt.legend(loc=\"best\") ; plt.xlabel(\"$x$\") ; plt.ylabel(\"$p(x)$\") ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## *Exercice 13* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 13.1 Travail à faire " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On s'intéresse à la solution $f$ du système différentiel $(1)$ suivant : $\\displaystyle \\left\\lbrace\\begin{array}{lll} \\forall r\\in]0,2], & \\displaystyle f^{\\prime\\prime}(r)+\\frac{1}{r}\\,f^{\\prime}(r)+25\\,f(r)=0&; \\\\ \\forall r\\in]2,5[ & \n", "\\displaystyle f^{\\prime\\prime}(r)+\\frac{1}{r}\\,f^{\\prime}(r)+4\\,f(r)=0&.\\end{array}\\right.$
la solution $f$ étant de classe $\\mathcal{C}^{1}$ sur l'intervalle $[0,5]$ et vérifiant les conditions aux limites : $f^{\\prime}(0)=0$ et $f(5)=1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On admet que la solution générale de l'équation différentielle $\\displaystyle f^{\\prime\\prime}(r)+\\frac{1}{r}\\,f^{\\prime}(r)+k^2\\,f(r)=0$ est une combinaison linéaire des deux fonctions $r\\mapsto J_0(k\\,r)$ et\n", "$r\\mapsto Y_0(k\\,r)$, où $J_0$ et $Y_0$ désignent les fonctions de Bessel de première et de seconde espèces :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "r,k = sb.symbols(\"r,k\", real=True, nonnegative=True)\n", "S1,S2 = sb.besselj(0,k*r),sb.bessely(0,k*r)\n", "(S1,S2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[ (S.diff(r,2) + (1/r)*S.diff(r) + k**2*S ).equals(0) for S in (S1,S2) ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'étude de ces deux fonctions en $r=0$ donne les résultats suivants :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "S1.replace(r,0), S1.diff(r).replace(r,0), S2.limit(r,0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On en déduit la solution générale du problème différentiel :" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Sur l'intervalle $[0,2]$ :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a,b,c = sb.symbols(\"a,b,c\", real=True)\n", "F1 = a*sb.besselj(0,5*r) ; F1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Sur l'intervalle $[2,5]$ :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F2 = b*sb.besselj(0,2*r)+c*sb.bessely(0,2*r) ; F2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$a)$ Vérifier que la solution proposée satisfait bien la condition limite à gauche. Définir ensuite un système de 3 équations à 3 inconnues $a$, $b$ et $c$ à partir des conditions de raccordement en $r=2$ de la fonction et de sa dérivée, et de la condition limite à droite." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sb.Eq(F1.diff(r).replace(r,0),0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "systeme = [ sb.Eq(F1,F2).replace(r,2), \\\n", " sb.Eq(F1.diff(r),F2.diff(r)).replace(r,2), \\\n", " sb.Eq(F2.replace(r,5),1) ]\n", "systeme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$b)$ En déduire les coefficients $a$, $b$ et $c$ qui satisfont le problème différentiel $(1)$, puis la solution sur chacun des intervalles $[0,2]$ et $[2,5]$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regle_abc = sb.solve(systeme, (a,b,c)) ; regle_abc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F_sol1 = F1.xreplace(regle_abc).simplify() ; F_sol1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F_sol2 = F2.xreplace(regle_abc).simplify() ; F_sol2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$c)$ Définir la solution sous forme d'une fonction numérique vectorisée $\\mathtt{f\\underline{\\phantom{m}}num}$ puis faire tracer sa courbe représentative. Les fonctions de Bessel numériques se trouvent dans le module $\\mathtt{scipy.special}$ et se nomment $\\mathtt{jn}$ et $\\mathtt{yn}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.special as sf" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mes_modules = [\"numpy\",{\"besselj\":sf.jn, \"bessely\":sf.yn}]\n", "F_num1 = sb.lambdify( r, F_sol1, modules=mes_modules )\n", "F_num2 = sb.lambdify( r, F_sol2, modules=mes_modules )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F_num_brouillon = lambda r : ((0<=r)&(r<=2))*F_num1(r) + ((2 **Compléments non exigibles** " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## *Exercice 12 : compléments sur $\\mathtt{lambdify}$* " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 12.2 Exemples complémentaires " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*$f)$ Autres exemples avec des fonctions qui ne sont pas dans $\\mathtt{numpy}$ (voir aussi le mémento).*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Fonction $\\mathtt{sympy.erfc}$**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x,t = sb.symbols(\"x,t\", real=True)\n", "Ix = sb.integrate( sb.exp(-t**2), (t,x,sb.oo) ).rewrite(sb.erfc).simplify()\n", "Ix" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "try :\n", " F = sb.lambdify( x, Ix, modules=[\"numpy\"])\n", " F( np.array([0.2,0.4]) )\n", "except Exception as e :\n", " print(\"Erreur :\",e)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(sb.erfc.__doc__[:161])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.special as sf # Special Functions\n", "print(sf.erfc.__doc__[:168])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Vérification :\n", "T = np.linspace(0,2.5,100)\n", "V1 = sf.erfc(T)\n", "V2 = np.array( [float(sb.erfc(x)) for x in T] )\n", "np.allclose(V1,V2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(T,V1-V2,\".m\") ; plt.grid() ;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F = sb.lambdify( x, Ix, modules=[{\"erfc\":sf.erfc}, \"numpy\"])\n", "F( np.array([0.2,0.4]) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "plt.plot(T,F(T),\"-g\") ; plt.grid() ; plt.xlabel(\"$t$\") ; plt.ylabel(\"$F(t)$\") ; " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* **Fonction $\\mathtt{sympy.lowergamma}$**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x,t = sb.symbols(\"x,t\", positive=True)\n", "Jx = sb.integrate( sb.sqrt(t)*sb.exp(-t**2), (t,0,x) ).simplify()\n", "Jx" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "try :\n", " G = sb.lambdify( x, Jx, modules=[\"numpy\"])\n", " G( np.array([0.2,0.4]) )\n", "except Exception as e :\n", " print(\"Erreur :\",e)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(sb.lowergamma.__doc__[:202])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(sf.gammainc.__doc__[:259])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def my_lowergamma(s,x) :\n", " return sf.gamma(s)*sf.gammainc(s,x)\n", "T = np.linspace(0,2.5,100)\n", "W1 = my_lowergamma(0.75,T) ; W1[:5]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "W2 = np.array([ float(sb.lowergamma(0.75,x)) for x in T ]) ; W2[:5]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.allclose(W1,W2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "G = sb.lambdify( x, Jx, modules=[{\"lowergamma\":my_lowergamma},\"numpy\"])\n", "G( np.array([0.2,0.4]) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(T,G(T),\"-r\") ; plt.grid() ; plt.xlabel(\"$t$\") ; plt.ylabel(\"$G(t)$\") ; " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Variante :\n", "G_bis = sb.lambdify( x, Jx, modules=[{\"lowergamma\": lambda s,x : sf.gamma(s)*sf.gammainc(s,x)},\"numpy\"])\n", "G_bis( np.array([0.2,0.4]) )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.allclose(G(T),G_bis(T))" ] } ], "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 }