

\documentclass[10pt,a4paper]{report}
\usepackage[height=255mm,width=183mm,margin=1cm,bottom=2cm,top=2cm]{geometry}
\usepackage{preambule}
\usepackage[upright]{fourier}

\usepackage{bbding,xlop}
 \usepackage{wasysym,marvosym,keystroke}
\graphicspath{{/home/moi/Photos/Tehessin/}{/home/moi/Figures/FigSTI/}{/home/moi/Figures/FigMaple/}}




\entete{TD Maple n°3}{Équations différentielles}{PCSI 2 - Lycée Clemenceau}%g,d,b

\usepackage{index}
\makeindex
\renewindex{default}{idx}{ind}{Index des commandes vues aujourd'hui}

\newcommand{\indexkeywordsmaple }[1]{\index[default]{maple@\textbf{\large Commandes}!#1@\textcolor{red}{\texttt{#1}}|textbf}}

\newcommand{\indn }[1]{\index[default]{maple@{\large
      Notions }!#1@ #1}}

\newcommand{\indt }[1]{#1\index[default]{maple@{\large
      Touches }!#1@ #1}}


%%% Environnement des réponses du logiciel
\DefineVerbatimEnvironment{boxedverbatim}{Verbatim}{fontsize=\small,formatcom=\color{blue},fontshape=it,gobble=2,frame=single,framesep=4mm,label=\textit{\small \textcolor{black}{Réponse du logiciel}},xleftmargin=10pt,rulecolor=\color{black}}




\newcommand\Com[1]{%
\shorthandoff{;:}\textcolor{red}{{ \ttfamily #1}}\index[default]{maple@\textbf{\large Commandes}!#1@\textcolor{red}{\texttt{#1}}|textbf}}


%%% Com avec syntaxe
\newcommand\Comsi[2]{%
\shorthandoff{;:}\textcolor{red}{{ \ttfamily #1}}(\textcolor{red}{\texttt{#2}})\index[default]{maple@\textbf{\large Commandes}!#1@\textcolor{red}{\texttt{#1}}|textbf}}





%%%%%%  Commande personnelle ou bout de code

\renewcommand{\code}[1]{\textcolor{red}{\texttt{#1}}}





\newcommand{\MAPLE}{\lstset{numbers=none,language=Maple,xleftmargin=10pt,%
keywordstyle =\color{red}\usefont{OT1}{cmtt}{b}{n},basicstyle=\ttfamily\color{red},commentstyle=\normalfont\scriptsize\itshape,
    backgroundcolor=\color{0.9white},index=[2][keywords],indexstyle=[2]\indexkeywordsmaple,breaklines=true}
}

\begin{document}
\setlength{\parindent}{0pt}
\thispagestyle{empty}

\setcounter{NumLecon}{3}
         \TITREC{TD Maple }{Équations différentielles}{0.9white}
        

\MAPLE














\SECTION{Syntaxe Maple}

Nous allons découvrir aujourd'hui comment Maple résout les équations différentielles et comment se débrouiller quand même Maple ne parvient pas à déterminer une solution explicite : voilà qui devrait éveiller l'attention des brillant(e)s physicien(ne)s que vous êtes.

\SUBSECTION{fonction - expression}

Nous aurns besoin de faire la différence entre une fonction et une expression, car Maple la fait...

Pour définir  une fonction, il faut la rentrer sous la  forme 

\begin{lstlisting}
 fonction:=variable->expression(variable);
\end{lstlisting}


Ainsi,
\begin{lstlisting}
> a:=x^2+2;
\end{lstlisting}

 est une expression, alors que
\begin{lstlisting}
> f:=x->x^2+2;
\end{lstlisting}



est une fonction.

Première différence :
 \begin{lstlisting}
a(32); f(32);                  
                      \end{lstlisting}


On peut transformer une expression en fonction à l'aide de

 \Comsi{unapply}{expression, variable} 

Par exemple, observez 

\begin{lstlisting}
> g:=unapply(a,x); g(32);
\end{lstlisting}


On peut évidemment définir des fonctions de plusieurs variables :

\begin{lstlisting}
> k:=(r,t)->(r*cos(t),r*sin(t)); k(32,Pi/3);
\end{lstlisting}



\SUBSECTION{diff - D}

Nous n'emploierons pas la même syntaxe selon que l'objet à dériver est une fonction ou une expression. Dans le cas d'une expression, on emploie \Comsi{diff}{expression, varaible1, variable2, ...}. Par exemple,
\begin{lstlisting}
> F:=sin(4*x^3+5*x);  diff(F,x);
\end{lstlisting}

On peut dériver plusieurs fois

\begin{lstlisting}
> diff(F,x,x,x);  diff(F,x$3);
\end{lstlisting}

pour dériver une fonction, on utilise l'opérateur \Comsi{D}{fonction}. Notez bien que \Comsi{d}{fonction} est une fonction. Par exemple, 
\begin{lstlisting}
f:=unapply(F,x);  D(f); D(f)(x); D(f)(32); (D@@3)(f); (D@@3)(f)(x);
\end{lstlisting}


\SUBSECTION{dsolve}

Maple sait résoudre pas mal d'équations différentielles grâce à

 \Comsi{dsolve}{\{equa diff, conditions initiales\}, fonction(variable)}. 

Par exemple,
\begin{lstlisting}
> dsolve({diff(y(t),t,t)+y(t)=0,D(y)(0)=1,y(0)=1},y(t));
\end{lstlisting}

est sans surprise. On peut même omettre les conditions initiales

\begin{lstlisting}
> dsolve({diff(y(t),t,t)+y(t)=0},y(t));
\end{lstlisting}

Au passage, notons que Maple permet de résoudre des équations différentielles dépendant de paramètres

\begin{lstlisting}
> dsolve({R*diff(q(t),t)+(q(t))/C=E,q(0)=0},q(t));
\end{lstlisting}


Testez par exemple Maple sur les équations suivantes
\begin{nliste}{3}
\item $(x+1)y'-xy=0$ \item $(1+x^2)y'+xy=3x^3+3x$ \item $(1+x^2)y'+xy=1$ \item $xy'+2y=x/(1+x^2)$ \item $y'\sqrt{1-x^2}-y^2-1=0$ \item $y''-3y'+2y=xe^x+\sin(x)$
\end{nliste}


\SUBSECTION{plot - odeplot}

L'instruction de base pour représenter graphiquement une fonction est 

\begin{lstlisting}
plot(expression(variable), variable=mini..maxi,options graphiques);
\end{lstlisting}

Par exemple

\begin{lstlisting}
> plot(ln(sinh(x))-x,x=0.1..10);
\end{lstlisting}

C'est parfois un peu farfelu

\begin{lstlisting}
> plot(cos(x)*exp(tan(x)),x=-4..4);
\end{lstlisting}

Il faut recadrer
\begin{lstlisting}
> plot(cos(x)*exp(tan(x)),x=-4..4,-4..4);
\end{lstlisting}

C'est déjà mieux. Mais il manque la touche finale 
\begin{lstlisting}
> plot(cos(x)*exp(tan(x)),x=-4..4,-4..4,title=`C'EST BEAU`);
\end{lstlisting}


Ce qui nous intéressera plus aujourd'hui, c'est de tracer des lignes brisées

\begin{lstlisting}
> plot([[0,0],[1,3],[2,-4],[3,2]]);
\end{lstlisting}

Notez bien l'utilisation des crochets : on a construit une liste de points différente de 


\begin{lstlisting}
> plot([[0,0],[2,-4],[1,3],[3,2]]);
\end{lstlisting}

On peut comparer deux tracés sur un même graphique en utilisant la fonction \Com{display} qu'il faut aller chercher dans la bibliothèque \Com{plots}

\begin{lstlisting}
> with(plots):
> d1:=plot([[0,0],[1,3],[2,-4],[3,2]],color=blue):
> d2:=plot(sin(2*x),x=0..3,style=point):
> display(d1,d2);
\end{lstlisting}


Pour terminer ce petit tour d'horizon, parlons de la fonction \Com{odeplot} qui, elle aussi, fait partie de la bibliothèque \Com{plots}. Il faut d'abord définir une solution d'une équation différentielle, par exemple
\begin{lstlisting}
> s:=dsolve({(x+1)*diff(y(x),x)-x*y(x)=0,y(0)=1},y(x),numeric):
\end{lstlisting}

L'argument \Com{numeric} ajouté à la fin impose à Maple de résoudre l'équation numériquement. Pour comprendre ce qui se passe, demandez

\begin{lstlisting}
> s(32);
\end{lstlisting}

Maple renvoie les coordonnées d'un point de la courbe solution. Pour tracer la courbe correspondante, il suffit d'entrer

\begin{lstlisting}
> odeplot(s,[x,y(x)],-0.5..2);
\end{lstlisting}

Le tracé semble manquer de précision. On peut rajouter comme option \Com{numpoints}\code{=1000} sachant que la valeur par défaut est 50.


Vous voici armé(e) pour entamer les hostilités




\SECTION{Résolution d'équations différentielles par discrétisation}

\SUBSECTION{Principe}


On considère une équation différentielle ordinaire (EDO pour les intimes) $u'(t)=f(t,u(t))$. Lorsqu'on ne connait pas de solution exacte à cette EDO, on essaye d'en avoir une bonne approximation par des méthodes numériques.

Nous ne nous pencherons pas ici sur le problème d'existence de solutions : nous supposerons que $f$ a les propriétés qui font des solutions existent là où il faut...
 
Nous allons par exemple \og observer\fg $u$ à intervalle de temps régulier~:~notons $h$ cette période d'échantillonnage.  

Nous allons essayer d'obtenir une suite d'approximations $u_n$ de $u(t_n)$ et d'évaluer la pertinence de cette approximation, le principe étant résumé sur  la figure \vref{discret} \vpageref{discret}



\begin{center}
 \includegraphics[height=4cm]{DiscretEDO.1}\caption{Principe de discrétisation d'une EDO}\label{discret}
\end{center}


Deux grandes possibilités s'offrent à nous pour \og se débarrasser\fg de $u'$ : 

\begin{itemize}
 \item soit remplacer une courbe par sa tangente, c'est à dire utiliser une approximation affine de $u$. On remplacera alors $u'(t)$ par une fonction dépendant de $u$ et ainsi obtenir une relation de récurrence définissant la suite $(u_n)$.
\item soit intégrer l'EDO : nous n'aurons plus de $u'$ ... mais il faudra calculer une intégrale ! Nous pourrons alors utiliser des  méthodes d'approximation  d'intégrales : méthodes des rectangles, des trapèzes, etc.
\end{itemize}

\SUBSECTION{Méthode d'Euler}

\SUBSUBSECTION{ Découverte sur un exemple}

Cette méthode, que vous avez dû étudier au lycée, a pour idée directrice de remplacer localement une courbe par sa tangente. Il y a des méthodes plus efficaces\footnote{Notamment celle de Runge-Kutta qu'utilise d'ailleurs Maple}, mais pour l'instant hors de notre portée mathématique.

Soit donc $u$ une fonction dérivable sur un intervalle $I=[a\pv b]$ vérifiant $u(0)=1$ et , pour tout $t$, $u'(t)=\alpha\cdot u(t)$ (oui, oui, c'est bon, moi aussi je sais  qui c'est, on n'a pas besoin de Maple pour ça, mais mieux vaut découvrir sur un exemple simple\footnote{Par exemple l'intensité des radiations émises par un corps radioactif est estimée en mesurant la concentration $u(t)$ d'un isotope instable. Soit $\tau$ sa demi-vie. Alors $u'(t)=\alpha\cdot u(t)\quad \text{avec} \ \alpha=-\fr{\ln 2}{\tau} $}...).

%Pour se simplifier la vie, nous réaliserons un maillage régulier de $I$ : nous le découperons en $n$ intervalles de longueur $\ofr{1}{n}$.

\begin{itemize}
\item Soit $t_n$ un élément  de $]a\pv b[$ tel que, pour $h$ un réel strictement positif suffisamment \og voisin\fg  de 0, $t_n+h$ soit encore dans I. On peut alors considérer que le quotient 

\[\fr{u(t_n+h)-u(t_n)}{h}\]

réalise une approximation d'autant plus valable de $u'(t_n)$ que $h$ est petit.

Que devient alors notre EDO \og approchée\fg ? En d'autres termes, trouver une relation de récurrence définissant $(u_n)$.

 Quelle est la nature de la suite $(u_n)$ ?

Est-on sûr(e) que cette suite réalise une bonne approximation de la solution de l'EDO ? 

\item Si vous avez des doutes, regardons des exemples.  Pour se simplifier la vie, nous réaliserons un maillage régulier de $I$ : nous le découperons en $N$ intervalles de longueur $\ofr{(b-a)}{N}$. Construisez alors une procédure dépendant des paramètres $\al$ et N qui compare\footnote{C'est à dire qui calcule $d_n=u(t_n)-u_n$} les valeurs de $u(t_n)$ et $u_n$ pour différentes valeurs de $n$.

\item Pour rendre encore plus visible ce phénomène, nous allons faire apparaître sur un même graphique des courbes obtenues par la méthode d'Euler et la solution exacte.

On construira pour cela une procédure \code{eu:=}\Com{proc}\code{(alpha,N,a,b)} qui ressemblera à 
\end{itemize}

\begin{lstlisting}
> with(plots);
>eu:=proc(alpha,N,a,b,u0)    
local ...:   
S:=NULL;   // On construit une suite de points vide au départ
y:=u0;  // y=u(0) au départ
for k from ... to N  do
x:=a+k*(b-a)/N; 

...
 
P(k):=[x,y]: // on définit le point P(k)
S:=S,P(k):
od:
Cnum:=plot([S]):
snum:=dsolve(...
Cexact:=odeplot(snum,...):
display(Cnum,Cexact);
end:
\end{lstlisting}


Regardez par exemple
\begin{lstlisting}
> eu(-4,6,0,3);
> eu(-4,24,0,3);
\end{lstlisting}
et commentez...




\SUBSUBSECTION{Méthode d'Euler : généralisation}

Oui, c'est bien beau, mais l'intérêt de l'informatique c'est de régler des cas les plus généraux possibles. 

Donc réfléchissons un peu ... $u'(t)=f(u(t),t)$ et $u'(t)\simeq \fr{u(t+h)-u(t)}{h}$ ... 

\begin{lstlisting}
 for k from 0 to N do
 y:=y+f(y,x)*(b-a)/N;
 x:=a+k*(b-a)/N;  
 P(k):=[x,y]:      
 S:=S,P(k):
 od:
\end{lstlisting}

Notez l'ordre de \code{x} et \code{y} ...


\SUBSUBSECTION{Estimation de l'erreur}

C'est un point extrêmement important mais les outils théoriques nous manquent pour l'instant pour en discuter


\SUBSUBSECTION{Sans filet}

Le but de notre nouveau jeu va être de discrétiser une équation différentielle du second ordre, de tracer le graphe correspondant et de comparer le résultat obtenu avec les solutions proposées par Maple grâce à \Com{dsolve} et \Com{odeplot}. Nous utiliserons le fait que pour une fonction $y$ deux fois dérivable sur un intervalle I et $h$ proche de 0, on a 

$$y(x+h)\simeq y(x)+y'(x)\times h$$

et donc 

$$y'(x+h)\simeq y'(x)+y''(x)\times h$$

Il s'agit maintenant de mettre au point une procédure \code{EquaDiff:=}\Com{proc}\code{(f,xmin,xmax,y0,yp0,n)} où $y0=y(0)$, $yp0=y'(0)$ et $f$ est la fonction de 3 variables définie par 

$$y''=f(y',y,x)$$

Le résultat doit faire apparaître sur le même graphique la courbe obtenue par la méthode d'Euler et celle obtenue à l'aide de \Com{odeplot}. 

Il ne restera plus qu'à tester votre procédure sur les nombreuses équations différentielles du second ordre que vous avez résolues en cours avec Madame Gornet.


\SUBSECTION{Méthode du point milieu}

Les anglo-saxons l'appellent schéma \emph{Leapfrog}\footnote{Pourquoi ce nom ?...}

En fait, pour une fonction suffisamment régulière sur l'intervalle $[a\pv b]$, l'idée est d'approcher la pente entre les points de coordonnées $\big(t-h,u(t-h)\big)$
et $\big(t+h,u(t+h)\big)$ par $u'(t)$.

Quand vous aurez étudié les formules de Taylor, vous pourrez montrer que l'approximation obtenue est meilleure puisque d'ordre 2, alors que la méthode d'Euler donne une approximation d'ordre 1.

Modifiez une des  procédures précédentes pour utiliser cette méthode. Attention~! Il y a trois niveaux de récurrence...


\begin{lstlisting}
 x:=a;
 ymoins:=u0; 
 y:=ymoins+h*f(ymoins,x); 
 for k from 1 to N do
 x:=a+k*h;
 yplus:=ymoins+2*h*f(y,x); 
 P(k):=[x+h,yplus]:      
 S:=S,P(k):
 ymoins:=y; 
 y:=yplus;
 od:
\end{lstlisting}

Observez :

\begin{lstlisting}
Mil((u,t)->-2*u+t,200,0,3,1);
Mil((u,t)->-2*u+t,600,0,3,1);
Mil((u,t)->4*u+t,200,0,3,1);
\end{lstlisting}


MORALITÉ : l'approche intuitive de la notion de limite ($\d x\simeq \Delta x$) permet de modéliser facilement certaines situations, mais 
si on pousse un peu, on ne peut se passer d'une théorisation rigoureuse pour expliquer certaines anomalies. Analyse numérique ne veut donc pas dire bidouillage... L'étude de l'estimation des erreurs commises doit être menée avec le plus grand soin mais est encore hors de notre portée. 

Du point de vue informatique, un autre facteur doit être pris en compte : la vitesse de convergence, pour ne pas avoir à attendre quelques siècles le résultat de la compilation...

\SUBSECTION{À partir de l'intégration numérique}




\SUBSUBSECTION{Premier exemple}



Soit $(E)$ l'équation différentielle $u'(t)+3u(t)=2$, $u(0)=1$. Nous savons parfaitement résoudre cette équation. Nous allons encore voir une méthode qui va vous paraître beaucoup plus compliquée, mais qui s'avérera fort utile
quand il s'agira de résoudre des équations différentielles plus sophistiquées où les méthodes habituelles de résolution ne nous
permettront plus de conclure mais que vous rencontrerez au hasard de votre carrière  d'ingénieur.

Nous allons <<~observer~>> la solution de l'équation à intervalle de temps régulier : la période d'échantillonnage qu'on notera encore $h$. On
fixera par exemple $h=0,2\, s$.

Nous allons intégrer l'équation $(E)$ sur une période d'échantillonnage quelconque, c'est à dire sur un intervalle $\big[kh,(k+1)h\big]$, avec
$k$ un entier naturel. Cela donne, grâce à la linéarité de l'intégration :

$$\int_{kh}^{(k+1)h}u'(t)\d t+3\int_{kh}^{(k+1)h}u(t)\d t=2\int_{kh}^{(k+1)h}\d t$$

On définit la suite $(u_k)_{k\in\bbn}$ par $u_k=u(kh)$.

Calculez le premier terme du premier membre ainsi que le second membre de l'équation : pas de problème.

Plus compliqué va être de calculer $\int_{kh}^{(k+1)h}u(t)\d t$ car nous ne connaissons pas $u$, donc encore moins une de ses
primitives. Nous allons malgré tout nous en tirer en déterminant une approximation du résultat.


\SUBSUBSECTION{Méthode des trapèzes}



Observez la figure   \vpageref{trap}

\begin{center}
\includegraphics[height=4cm]{DiscretEDO.2}\caption{Méthode des trapèzes}\label{trap}
\end{center}

Il semble raisonnable de penser que plus  $h$ sera petit, plus l'aire sombre le sera également et donc meilleure sera l'approximation.

Maintenant, en vous souvenant de la formule $\AR(\text{trapèze})=\fr{(\text{petite base}+\text{grande base})\times\text{hauteur}}{2}$,
donnez l'aire du trapèze bleu.

\SUBSUBSECTION{\'Equation aux différences}



Déduisez-en que la suite $(u_n)$ vérifie l'équation aux différences $(F)$
$$\fr{13}{10}u_{n+1}-\fr{7}{10}u_n=\fr{2}{5}$$
Vous reconnaissez une suite arithmético-géométrique. Exprimez $u_n$ en fonction de $n$.


\SUBSUBSECTION{\'Equation différentielle vs équation aux différences}


Pour vérifiez si la solution de l'équation aux différences n'est pas très éloignée de la solution théorique, nous allons observer la
solution de $(E)$ toutes les 0,2 s : on pose $y_n=f(nh)$. Donnez l'expression de $y_n$.

\sm

On s'attend à ce que $u_n$ et $y_n$ soient très proches mais pas égaux. En effet, le passage de $(E)$ à $(F)$ s'opère à l'aide d'une
approximation (celle des trapèzes) et non d'une équivalence.

On pose $d_n=y_n-u_n$. À l'aide  de Maple, dressez un tableau où
figureront les 10 premiers termes des suites $(u_n)$, $(y_n)$ et $(d_n)$. Cela vous rassure-t-il~?

Montrez que $(d_n)$ est en fait la différence entre deux suites géométriques. Utilisez ce résultat pour confirmer ou infirmer votre
conjecture.

\SUBSUBSECTION{\'Epilogue : les maths, c'est pas Harry Potter}

\bigskip

Tout serait donc magique en maths, ou bien n'avez-vous été confronté(e) toute votre vie qu'à des problèmes truqués~?
\sm

Soit $(E')$ l'équation différentielle $2u'(t)-3u(t)=t$, $u(0)=2$. Vous pouvez demander à Maple de vérifier que la solution de cette équation est la fonction
$s$ définie par $$u(t)=-\fr{1}{3}t-\fr{2}{9}+\fr{20}{9}{\rm e}^{3/2t}$$ Vous pourrez peut-être un jour montrer que la suite solution de
l'équation aux différences associée à $(E)$ est $$u_n=\fr{20}{9}\pa{\fr{43}{37}}^n-\fr{1}{30}n-\fr{2}{9}$$

On pose encore $d_n=y_n-u_n$. Pouvez-vous conclure comme dans l'exemple précédent~? Moralité~?

\SUBSUBSECTION{Comparaison avec le méthode d'Euler}

Cette méthode est-elle plus précise que la méthode d'Euler vue précédemment ? 


\SUBSECTION{Généralisation}

En dehors des équations linéaires, la généralisation devient extrêmement difficile... Nous nous contenterons donc des équations du style 

\[u'(t)=\alpha u(t)+\beta t+\gamma\]

Écrivez donc une procédure ayant  pour arguments $\alpha$, $\beta$ et $\gamma$, les bornes $a$ et $b$ de l'intervalle, le nombre N de \og mailles du filet\fg et la valeur de $u(0)$.






\end{document}