Skip to content
Snippets Groups Projects
Forked from orestis.malaspin / math_tech_info
339 commits behind the upstream repository.
cours.md 180.02 KiB
author:
- Orestis Malaspinas
title: Résumé du cours de Mathématiques
autoSectionLabels: true
autoEqnLabels: true
eqnPrefix: 
  - "éq."
  - "éqs."
chapters: true
numberSections: true
chaptersDepth: 1
sectionsDepth: 3
lang: fr
documentclass: book
papersize: A4
cref: false
urlcolor: blue

Rappel

Fonctions

Une fonction f de façon générale est un objet qui prend un (ou plusieurs) paramètres et qui lui (leur) associent un (ou plusieurs) résultats. \mbox{résultat}=f(\mbox{paramètres}).

  1. La tension U est une fonction de la résistance R et du courant I \begin{aligned} U=f(R,I)=R\cdot I.\end{aligned}

  2. Une fonction peut être quelque chose de beaucoup plus général (qu’on ne peut pas forcément représenter simplement avec des opérateurs mathématiques). Prenons le cas de la fonction qui pour un nombre entier x rend le prochain entier qui commence par la même lettre que x. f(2)=10,\ f(3)=13,\ ...

Dans ce cours nous allons nous intéresser à des fonctions à un seul paramètre (aussi appelé variable). Si on note la variable x et le résultat y, de façon générale on peut écrire y = f(x). Si par ailleurs on a une fonction g et une fonction f, on peut effectuer des compositions de fonction, qu’on note g\circ f, ou encore y=g(f(x)).

  1. Soit f(x)=2\cdot x et g(x)=\sqrt{x}, alors la composition des deux fonctions f(g(x))=f(\sqrt{x})=2\sqrt{x}.

  2. On peut composer un nombre arbitraire de fonctions. Voyons le cas avec trois fonctions f(x)=2x^2+3, g(x)=\cos(2\cdot x), et h(x)=1/x f(g(h(x)))=f(g(1/x))=f(\cos(2/x))=2\cos^2(2/x)+3.

Pour certaines fonctions, notons les f(x), on peut également définir une fonction inverse que l’on note f^{-1}(x) dont la composition donne la variable de départ f(f^{-1}(x))=x.

  1. Soient f(x)=2\cdot x et f^{-1}(x)=x/2, alors la composition des deux fonctions f(f^{-1}(x))=f(x/2)=2x/2=x.

  2. Soient f(x)=x^2 et f^{-1}(x)=\sqrt{x}, alors la composition des deux fonctions f(f^{-1}(x))=f(\sqrt{x})=|x|. On a donc que \sqrt{x} est l’inverse de x^2 uniquement pour les réels positifs. f(x)=x^2 n’a pas d’inverse pour les x négatifs.

Domaine de définition

Le domaine de définition, noté D\subseteq{\mathbb{R}}, d’une fonction f, est l’ensemble de valeurs où f admet une image.

  1. Le domaine de définition de f(x)=x est D={\mathbb{R}}.

  2. Le domaine de définition de f(x)=1/x est D={\mathbb{R}}^\ast.

  3. Le domaine de définition de f(x)=\sqrt{x+1}/(x-10) est D=[-1;10[\cup]10;\infty[.

Limites

Soit f une fonction et D\subseteq{\mathbb{R}} non-vide et non réduit à un point et soient a et b deux réels.

Limite

Pour f définie en D, sauf peut-être en a, on dit que b est la limite de x en a si \lim\limits_{x\rightarrow a}f(x)=b. C’est-à-dire si on a un voisinage de b qui contient toutes la valeurs de f(x) pour x proche de a.

Si f est définie en a alors on a \lim\limits_{x\rightarrow a}=f(a).

  1. Si f(x)=x, alors \lim\limits_{x\rightarrow 0}f(x)=0.

Pour f définie en D, sauf peut-être en a, et c un réel positif. On dit que la limite de f(x) en a tend vers l’infini si l’intervalle [c;\infty[ contient toutes les valeurs de f(x) pour x proche de a.

  1. Si f(x)=1/x^2, alors \lim\limits_{x\rightarrow 0}f(x)=\infty.

Limite à gauche, limite à droite

Pour certaines fonctions, il est possible que le comportement de celles-ci soit différent selon qu’on approche par la gauche ou par la droite (i.e. f(x)=1/x).

On note la limite à droite \lim\limits{x\rightarrow a^+} f(x) ou \lim\limits_{x\rightarrow a,x>a} f(x) et \lim\limits{x\rightarrow a^-} f(x) ou \lim\limits_{x\rightarrow a,x<a} f(x) la limite à gauche de la fonction f en a.

Si la fonction f admet une limite en a, alors les deux limites doivent être égales.

  1. Si f(x)=1/x, alors \lim\limits_{x\rightarrow 0^+} f(x)=\infty et \lim\limits_{x\rightarrow 0^-} f(x)=-\infty.

Asymptotes

Dans certains cas il peut être intéressant d’étudier le comportement des fonctions quand x\rightarrow\pm\infty. Dans ces cas-là on dit qu’on s’intéresse au comportement asymptotique d’une fonction. Ce concept est particulièrement relevant quand on étudie une fonction que a la forme d’une fraction h(x)=\frac{f(x)}{g(x)}. Si on s’intéresse au comportement à l’infini de cette fonction on va prendre sa “limite” lorsque x\rightarrow\infty \lim_{x\rightarrow\infty} h(x)=\lim_{x\rightarrow\infty}\left(\frac{f(x)}{g(x)}\right). Un exemple peut être f(x)=x-1, g(x)=x+1 et donc h(x)=(x-1)/(x+1) \lim_{x\rightarrow\infty} \frac{x-1}{x+1}=\lim_{x\rightarrow\infty} \frac{x}{x}=1. De même quand on a f(x)=3x^4-5x^3+1, g(x)=1 et donc h(x)=3x^4-5x^3+1. Il vient donc \lim_{x\rightarrow\infty} 3x^4-5x^3+1=\lim_{x\rightarrow\infty}3x^4=\infty.

Si nous compliquons un peu l’exemple, et que nous avons f(x)=x^3+3x^2+1, g(x)=x^2 et donc h(x)=(x^3+3x^2+1)/x^2 \lim_{x\rightarrow\infty} (x^3+3x^2+1)/x^2=\lim_{x\rightarrow\infty} x=\infty. Ce genre d’estimations est imporant en informatique lors de l’analyse de performance des algorithmes. On peut prendre l’exemple des algorithmes de tri “bubble sort” et “quick sort”. Leur complexité respective moyenne est de n^2 et de n\log(n), quand n est le nombre d’éléments de la chaîne à trier. Si on fait le rapport pour de ces deux complexités on a \lim_{n\rightarrow\infty} \frac{n^2}{n\log(n)}=\lim_{n\rightarrow\infty} \frac{n}{\log(n)}. On peut simplement voir que ce rapport va tendre vers l’infini en dessinant la courbe n/\log(n). Il existe un moyen “analytique” d’évaluerce rapport. Tout nombre n peut s’écrire avec une précision p comme n=A\cdot 10^{p-1},p est le nombre de chiffres significatifs qu’on veut représenter, et 1\leq A< 10. On a également que1 \log(A)=\log\left(\frac{1+y}{1-y}\right)=2\sum_{k=0}^\infty \frac{y^{2k+1}}{2k+1}, avec y=(A-1)/(A+1). On a finalement que \log(n)=\log(A\cdot 10^{p-1})=(p-1)\log(10)+2\sum_{k=0}^\infty \frac{y^{2k+1}}{2k+1}. La valeur de y étant quelque chose de proche de 0, la somme converge vite vers une valeur finie et on peut faire l’approximation \log(n)\cong(p-1)\log(10), pour n grand (ce qui est équivalent à p grand). On a donc que finalement le rapport n/\log(n) va comme \lim_{n\rightarrow\infty}\frac{n}{\log(n)}=\frac{A}{\log(10)}\cdot\lim_{n\rightarrow\infty}\frac{10^{p-1}}{(p-1)}=\frac{A}{\log(10)}\cdot\lim_{n\rightarrow\infty}\frac{10^{p-1}}{p}=\infty.

Continuité

Soit f une fonction définie sur un intervalle ouvert D contenant a. On dit que f est continue en a si et seulement si \lim\limits_{x\rightarrow a}f(x)=f(a).

Soient f et g deux fonctions continues en a alors et b un réel:

  1. f+g est continue en a.

  2. b f est continue en a.

  3. si g(a)\neq 0, f/g est continue en a.

  4. h=g\circ f est continue en a.

Une fonction f est dite continue dans un intervalle D=]a;b[ si et seulement si elle est continue en tout point de D. De plus, elle est continue sur D=[a,b] si elle est continue sur ]a;b[ et continue à droite en a et à gauche en b.

[Théorème des valeurs intermédiaires] Soit f une fonction continue sur D, et a,b deux points contenus dans D tels que a<b et f(a)<f(b), alors \forall y\in [f(a);f(b)],\ \exists\ c|f(c)=y.

Dérivées

Soit f une fonction définie sur D et a\in D. On dit que f est dérivable en a s’il existe un b (appelé la dérivée de f en a) tel que \begin{aligned} &\lim\limits_{h\rightarrow 0}\frac{f(a+h)-f(a)}{h}=b,\hbox{ ou}\\ &\lim\limits_{x\rightarrow a}\frac{f(x)-f(a)}{x-a}=b.\end{aligned}

Si f est dérivable en tout point de D=]a;b[, alors on définit f' la fonction dérivée de f dans l’intervalle D qui associe en tout point x de D la valeur dérivée de f.

Si f est dérivable en a alors f est continue en a.

Soient f et g deux fonctions dérivables (dont les dérivées sont f' et g'), et a\in{\mathbb{R}}, alors

  1. (f+g)'=f'+g'.

  2. (af)'=a f'.

  3. (f\cdot g)'=f'g+fg'.

  4. Si g ne s’annule pas (f/g)'=(f'g-fg')/g^2.

  5. (g\circ f)'=(g'\circ f)\cdot f'.

Il existe quelques dérivées importantes que nous allons utiliser régulièrement dans la suite de ce cours. En supposons que C\in {\mathbb{R}}, nous avons

  1. f(x)=x^n, f'(x)=nx^{n-1} .

  2. f(x)=e^{C x}, f'(x)=Ce^{Cx}.

  3. f(x)=\ln(x), f'(x)=\frac{1}{x}.

  4. $f(x)=C, $f’(x)=0.

  5. f(x)=\sin(x), f'(x)=\cos(x).

  6. f(x)=\cos(x), f'(x)=-\sin(x).

Si f' est dérivable sur D, alors sa dérivée, notée f'', est appelée la dérivée seconde de f.

Variation des fonctions

Soit f' la fonction dérivée de f sur D

  1. Si f'>0 sur D, alors f est croissante sur D.

  2. Si f'<0 sur D, alors f est décroissante sur D.

  3. Si f'=0 sur D, alors f est constante sur D.

Une fonction admet un maximum local (respectivement minimum local) sur un intervalle D=]a;b[ s’il existe un x_0 tel que f(x_0)\geq f(x) (respectivement f(x_0)\leq f(x)) pour tout x\in D.

Soit f une fonction dérivable sur D=]a;b[ et x_0\in D. Si f admet un maximum ou un minimum en x_0 alors f'(x_0)=0. De plus si f'(x_0)=0 et f' change de signe en x_0 alors f(x_0) est un maximum ou un minimum de f.

Etude de fonction

Effectuer l’étude de fonction de la fonction suivante f(x)=\frac{x^3}{x^2-4}.

  1. Déterminer le domaine de définition.

  2. Déterminer la parité de la fonction. Rappel: \begin{aligned} f(-x)&=f(x),\ \mbox{paire},\\ f(-x)&=-f(x),\ \mbox{impaire}. \end{aligned}

  3. Trouver les zéros de la fonction (Indication: trouver les x tels que f(x)=0).

  4. Trouver les éventuelles asymptotes verticales ou disconsinuités, ainsi que les asymptotes affines.

  5. Caluler f'(x) et déterminer sa croissance et points critiques (déterminer où la fonction est croissante, décroissante, atteint un extremum, etc).

  6. Faire un croquis de f(x).

Intégrales

Interprétation géométrique

Dans ce chapitre nous nous intéressons au calcul d’aires sous une fonction f. La fonction f satisfait les hypothèses suivantes.

  1. f(x) est bornée dans l’intervalle [a,b]\in{\mathbb{R}}.

  2. f(x) est continue presque partout.

Nous définissions également l’infimum de f sur un intervalle [x_0,x_1], noté \inf\limits_{[x_0,x_1]} f(x) comme étant la valeur bornant par dessous toutes les valeurs prises par f(x) dans l’intervalle [x_0,x_1]. Le suprémum sur un intervalle [x_0,x_1], noté \sup\limits_{[x_0,x_1]} f(x) comme étant la valeur bornant par dessus toutes les valeurs prises par f(x) dans l’intervalle [x_0,x_1].

Finalement nous définissons une subdivision \Delta_n=\{a=x_0<x_1<...<x_{n-1}<x_{n}=b\} est une suite finie contenant n termes dans [a,b].

On peut à présent approximer l’aire sous la fonction f(x) dans l’intervalle [a,b] de deux façon:

  1. A^i(n)=\sum_{i=0}^{n-1} \inf\limits_{[x_i,x_{i+1}]} f(x) (x_{i+1}-x_i) comme étant l’aire inférieure.

  2. A^s(n)=\sum_{i=0}^{n-1} \sup\limits_{[x_i,x_{i+1}]} f(x) (x_{i+1}-x_i) comme étant l’aire supérieure.

L’aire de sous la fonction f(x) est donnée par la limite pour n\rightarrow\infty de A^i ou A^s (si elle existe).

  1. Ces sommes peuvent être positives ou négatives en fonction du signe de f.

  2. Une implantation informatique est immédiate.

Une fonction est dite intégrable au sens de Riemann si \lim\limits_{n\rightarrow\infty}A^i(n)=\lim\limits_{n\rightarrow\infty}A^s(n)=\int_a^b f(x){\mathrm{d}}x.

Dans la formule \int_a^b f(x){\mathrm{d}}x Ici x est appelée variable d’intégration, a et b sont les bornes d’intégration. Pour des raisons de consistance dans les notations la variable d’intégration ne peut être désignée avec le même symbole qu’une des bornes d’intégration.

Intégration de f(x)=x dans intervalle [0,1].

Il est trivial de calculer que cette aire vaut 1/2 (c’est l’aire d’un triangle rectangle de côté 1). Néanmoins, évaluons également cette aire à l’aide de A^i et A^s. Commençons par subdiviser [0,1] en n intervalles égaux de longueur \delta=1/n. Comme f(x) est strictement croissante, on a que \inf\limits_{[x_i,x_{i+1}]}f(x)=f(x_i) et que \sup\limits_{[x_i,x_{i+1}]}f(x)=f(x_{i+1}). On a donc que

  1. A^i(n)=\delta\sum_{i=0}^{n-1} x_i=\delta\sum_{i=0}^{n-1}\frac{i}{n}=\frac{n(n-1)}{2n^2}=\frac{n-1}{2n}2. Et donc en prenant la limite pour n\rightarrow\infty il vient A^i=\lim\limits_{n\rightarrow\infty}\frac{n-1}{2n}=\frac{1}{2}.

  2. A^s(n)=\delta\sum_{i=0}^{n-1} x_{i+1}=\delta\sum_{i=0}^{n-1}\frac{i+1}{n}=\delta\sum_{i=0}^{n}\frac{i}{n}=\frac{n(n+1)}{2n^2}=\frac{n+1}{2n}. Et donc en prenant la limite pour n\rightarrow\infty il vient A^s=\lim\limits_{n\rightarrow\infty}\frac{n+1}{2n}=\frac{1}{2}.

Calculer l’aire sous la courbe de f(x)=x^2 dans intervalle [0,1].

Indication: \sum_{i=0}^n i^2=\frac{1}{6}n(n+1)(2n+1).

Interprétation physique

Supposons que nous ayons une fonction, x(t), qui donne la position d’un objet pour un intervalle de temps t\in[a,b]. Nous pouvons aisément en déduire la vitesse v(t) de l’objet, comme étant la variation de x(t) pour tout t. Autrement dit v(t)=x'(t).

Supposons à présent que nous ne connaissions que la vitesse v(t) de notre objet. Afin de déduire sa position nous prendrions un certain nombre d’intervalles de temps \delta t_i=t_{i+1}-t_i que nous multiplierions par v(t_i) afin de retrouver la distance parcourue pendant l’intervalle \delta t_i et ainsi de suite. Afin d’améliorer l’approximation de la distance parcourue nous diminuerions la valeur de \delta t_i jusqu’à ce que \delta t_i\rightarrow 0.

Nous voyons donc que cette méthode, n’est autre qu’une façon “intuitive” d’intégrer la vitesse afin de trouver la position. Et que donc l’intégrale et la dérivée sont étroitement liée. La vitesse étant la dérivée de la position et la position étant l’intégrale de la vitesse.

Primitive

Si maintenant nous essayons de généraliser le calcul de l’intégrale d’une fonction, il s’avère que le calcul d’une intégrale est l’inverse du calcul d’une dérivée.

Soit f une fonction. On dit que F est la primitive de f sur l’intervalle D\subseteq{\mathbb{R}} si F'(x)=f(x) \forall x\in D.

Si F est une primitive de f, alors on peut définir la fonction G telle que G(x)=F(x)+C, \forall C\in{\mathbb{R}} qui est aussi une primitive de f. On dit donc que la primitive de f est définie à une constante additive près. En effet, si F'=f on a G'=F'+\underbrace{C'}_{=0}=F'=f.

S’il existe a\in D et b\in{\mathbb{R}} alors il existe une unique primitive F telle que F(a)=b.

Soit f(x)=x, alors l’ensemble de primitives correspondantes est G=x^2/2+C. Si nous cherchons la primitive telle que G(0)=0, il vient que C=0 et donc la primitive est unique et vaut F(x)=x^2/2.

Calculez les primitives suivantes (Indication: Il s’agit de trouver les fonctions F(x) telles que F'(x)=f(x)):

  1. F(x)=\int x^2{\mathrm{d}}x.

  2. F(x)=\int x^n{\mathrm{d}}x, n\in {\mathbb{R}}\backslash\{-1\}.

  3. F(x)=\int \sqrt{x}{\mathrm{d}}x.

  4. F(x)=\int \frac{1}{x}{\mathrm{d}}x.

  5. F(x)=\int \exp(x){\mathrm{d}}x.

  6. F(x)=\int \sin(x){\mathrm{d}}x.

Maintenant que vous avez calculé toutes ces primitives de base, nous pouvons récapituler des formules qui seront importantes pour la suite:

  1. \int x^n{\mathrm{d}}x=\frac{1}{n+1}x^{n+1}+C, n\in {\mathbb{R}}\backslash\{-1\}.

  2. \int \frac{1}{x}{\mathrm{d}}x=\ln(x)+C.

  3. \int \exp(x){\mathrm{d}}x=\exp(x)+C.

  4. \int \sin(x){\mathrm{d}}x=-\cos(x)+C.

  5. \int \cos(x){\mathrm{d}}x=\sin(x)+C.

[def_prim] En définissant à présent l’intégrale à l’aide de la notion de primitive, nous avons que pour a,b\in{\mathbb{R}} et a<b \int_a^b f(x){\mathrm{d}}x=\left.F\right|_a^b=F(b)-F(a).

On dit que x est la variable d’intégration. Elle est dite “muette” car elle disparaît après que l’intégrale ait été effectuée. On peut donc écrire l’équation ci-dessus de façon équivalente en remplaçant le symbole x par n’importe quelle autre lettre (sauf a,b,f,F).

On notera que la constante additive C a disparu de cette formule. En effet, remplaçons F par G=F+C, il vient \int_a^b f(x){\mathrm{d}}x=G(b)-G(a)=F(b)+C-F(a)-C=F(b)-F(a).

De la définition [def_prim], il vient immédiatement que \int_a^af(x){\mathrm{d}}x=F(a)-F(a)=0 et que \int_a^bf(x){\mathrm{d}}x= -\int_b^af(x){\mathrm{d}}x

Nous pouvons à présent définir la fonction G(x) telle que G(x)=\int_a^xf(y){\mathrm{d}}y=F(x)-F(a). Nous avons donc que G(x) est la primitive de f telle que G(a)=0.

Soient f et g deux fonctions intégrables sur un intervalle D=[a,b]\subseteq{\mathbb{R}}, c\in[a,b], et \alpha\in{\mathbb{R}}. On a

  1. La dérivée et l’intégrale “s’annulent” \left(\int_a^x f(x){\mathrm{d}}x\right)'=\left(F(x)-F(a)\right)'=F'(x)-\left(F(a)\right)'=F'(x)=f(x).

  2. La fonction h=f+g admet aussi une primitive sur D, et on a \int_a^b(f(x)+g(x)){\mathrm{d}}x=\int_a^b f(x){\mathrm{d}}x+\int_a^b g(x){\mathrm{d}}x.

  3. La fonction h=\alpha f admet aussi une primitive sur D, et on a \int_a^b\alpha f(x){\mathrm{d}}x=\alpha\int_a^b f(x){\mathrm{d}}x.

  4. Relation de Chasles (faire la démonstration en exercice) \int_a^c f(x){\mathrm{d}}x=\int_a^b f(x){\mathrm{d}}x+\int_b^c f(x){\mathrm{d}}x. De cette relation on déduit qu’on peut calculer l’intégrale d’une fonction continue par morceaux sur [a,b].

  5. Si f est paire alors \int_{-a}^a f(x){\mathrm{d}}x = 2\int_0^a f(x){\mathrm{d}}x.

  6. Si f est impaire alors \int_{-a}^a f(x){\mathrm{d}}x = 0.

Intégrales impropres

Si une des bornes d’intégration ou si la fonction à intégrer admet une discontinuité à des points bien définis, nous parlons intégrales impropres.

Lorsqu’une borne d’intégration est infinie, alors nous pouvons avoir les cas de figures suivants \begin{aligned} &\int_a^\infty f(x){\mathrm{d}}x=\lim\limits_{b\rightarrow\infty}\int_a^b f(x){\mathrm{d}}x,\\ &\int_{-\infty}^b f(x){\mathrm{d}}x=\lim\limits_{a\rightarrow\infty}\int_{-a}^b f(x){\mathrm{d}}x,\\ &\int_{-\infty}^\infty f(x){\mathrm{d}}x=\lim\limits_{a\rightarrow\infty}\int_{-a}^a f(x){\mathrm{d}}x.\end{aligned}

Calculer l’intégrale suivante \int_0^\infty e^{-ax}{\mathrm{d}}x,\quad a>0. Nous pouvons réécrire l’intégrale ci-dessus comme \int_0^\infty e^{-ax}{\mathrm{d}}x=\lim\limits_{b\rightarrow \infty}\int_0^b e^{-ax}{\mathrm{d}}x=-\frac{1}{a}\lim\limits_{b\rightarrow\infty}\left[e^{-ax}\right]_0^b=-\frac{1}{a}\left[\lim\limits_{b\rightarrow \infty}e^{-ab}-1\right]=\frac{1}{a}.

Calculer l’intégrale suivante \int_1^\infty \frac{1}{x^2}{\mathrm{d}}x.

Lorsque nous avons une discontinuité dans la fonction f au point c\in[a,b] nous avons \int_a^b f(x){\mathrm{d}}x = \lim\limits_{\varepsilon\rightarrow 0}\int_a^{c-\varepsilon} f(x){\mathrm{d}}x +\int_{c+\varepsilon}^b f(x){\mathrm{d}}x.

Montrer que \int_{-1}^2\frac{1}{x}=\ln{2}.

Soit une fonction f admettant une primitive sur [a,b] et a<b, alors on appelle la valeur moyenne de cette fonction sur [a,b], \bar{f}, \bar{f}=\frac{1}{b-a}\int_a^bf(x){\mathrm{d}}x.

Méthodes d’intégration

Dans cette section, nous allons étudier différentes méthodes pour intégrer des fonctions.

Intégration de fonctions usuelles et cas particuliers

Le calcul d’une primitive ou d’une intégrale n’est en général pas une chose aisée. Nous connaissons les formules d’intégration pour certaines fonctions particulières.

Polynômes

Les polynômes s’intègrent terme à terme. Pour \{a_i\}_{i=0}^{n}\in{\mathbb{R}} \begin{aligned} &\int a_0 + a_1 x + a_2 x^2+\cdots+a_{n-1} x^{n-1}+a_{n} x^{n}{\mathrm{d}}x\\ =&\int a_0{\mathrm{d}}x + \int a_1 x{\mathrm{d}}x + \int a_2 x^2{\mathrm{d}}x+\cdots+\int a_{n-1} x^{n-1}{\mathrm{d}}x+\int a_{n} x^{n}{\mathrm{d}}x\\ =&a_0 x + \frac{a_1}{2}x^2+\frac{a_2}{3}x^3+\cdots+\frac{a_n}{n+1}x^{n+1}+c.\end{aligned}

Intégrer la fonction suivante \int (x+2)(x^3+3x^2+4x-3){\mathrm{d}}x.

Application de la règle de chaîne pour l’intégration

Une primitive de la forme \int f(x)f'(x){\mathrm{d}}x=\frac{1}{2}f(x)^2+c.

Le calcul de la primitive suivante \int \sin(x)\cos(x){\mathrm{d}}x=\frac{1}{2}\sin^2(x)+c.

Inverse de la dérivation logarithmique

Une primitive de la forme \int \frac{f'(x)}{f(x)}{\mathrm{d}}x=\ln(f(x))+c.

Le calcul de la primitive de suivante \int \frac{1}{x}{\mathrm{d}}x=\int \frac{(x)'}{x}{\mathrm{d}}x=\ln(x)+c.

Règle de chaîne

De façon une des façons les plus simples de calculer une primitive est de reconnaître la règle de chaîne dans le terme à intégrer \int g'(f(x))f'(x){\mathrm{d}}x=\int [g(f(x))]' {\mathrm{d}}x=g(f(x))+c.

Si g est définie comme g(x)=x^{-1} et f(x)=3x^2+2, alors la primitive \int \frac{f'(x)}{g'(f(x))}{\mathrm{d}}x=\int -\frac{6 x}{(3x^2+2)^2}{\mathrm{d}}x=\frac{1}{3x^2+2}.

Intégration par parties

La dérivation d’un produit de fonction f\cdot g s’écrit (f(x)g(x))'=f'(x) g(x)+f(x) g'(x). En intégrant cette équation on obtient f(x)g(x)=\int f'(x) g(x){\mathrm{d}}x+\int f(x) g'(x){\mathrm{d}}x. Une primitive de la forme \int f'(x) g(x){\mathrm{d}}x peut se calculer de la façon suivante \int f'(x) g(x){\mathrm{d}}x=f(x)g(x)-\int f(x) g'(x){\mathrm{d}}x. De façon similaire si nous nous intéressons à une intégrale définie \int_a^b f'(x) g(x){\mathrm{d}}x=\left.(f(x)g(x))\right|_a^b-\int_a^b f(x) g'(x){\mathrm{d}}x. Le choix des fonctions est complètement arbitraire. Néanmoins, le but de cette transformation est de remplacer une intégrale par une autre dont on connaîtrait la solution.

Des “règles” pour utiliser cette technique seraient que

  1. g' soit facile à calculer et aurait une forme plus simple que g.

  2. f=\int f'{\mathrm{d}}x soit facile à calculer et aurait une forme plus simple que f'.

Calculer les primitives suivantes

  1. \int x e^x{\mathrm{d}}x. g(x)=x, f'(x)=e^x et donc g'(x)=1, f(x)=e^x. Il vient \int x e^x=x e^x-\int e^x{\mathrm{d}}x=x e^x-e^x+c.

  2. \int \cos(x)\sin(x){\mathrm{d}}x. g= \cos(x), f'(x)=\sin(x) et donc g'(x)=-\sin(x), f(x)=-\cos(x). Il vient \begin{aligned} &\int \cos(x)\sin(x){\mathrm{d}}x=\sin^2(x)-\int \cos(x)\sin(x){\mathrm{d}}x\nonumber\\ \Rightarrow &\int \cos(x)\sin(x){\mathrm{d}}x=\frac{1}{2}\sin^2(x). \end{aligned} On voit que le résultat de l’intégration par partie nous redonne l’intégrale de départ. Ceci nous permet d’évaluer directement la dite intégrale.

Il est également possible d’enchaîner plusieurs intégrations par parties.

L’intégrale de \int x^2 e^x{\mathrm{d}}x. En posant g(x)=x^2, f'(x)=e^x et donc g'(x)=2x, f(x)=e^x. Il vient \int x^2 e^x{\mathrm{d}}x=x^2e^x-2\int x e^x{\mathrm{d}}x. On pose de façon similaire g(x)=x, f'(x)=e^x et donc g'(x)=1, f(x)=e^x et il vient \int x^2 e^x{\mathrm{d}}x=x^2e^x-2\left(x e^x -\int e^x{\mathrm{d}}x\right)=x^2e^x-2x e^x +2e^x+c.

Calculer les primitives suivantes

  1. \int \ln(x){\mathrm{d}}x

  2. \int x^2 \sin(x){\mathrm{d}}x

  3. \int e^x\sin(x){\mathrm{d}}x

Intégration par changement de variables

On observe que la dérivation de la composition de deux fonctions F et g est donnée par (F\circ g)'=(f\circ g)\cdot g',\mbox{ ou } [F(g(y))]'=f(g(y))\cdot g'(y),f=F'. Si nous intégrons cette relation on obtient \begin{aligned} \int_a^b f(g(y))g'(y){\mathrm{d}}y = \int_a^b [F(g(y))]'{\mathrm{d}}y=\left.F(g(y))\right|_a^b=F(g(b))-F(g(a))=\int_{g(a)}^{g(b)}f(x){\mathrm{d}}x.\end{aligned} Cette relation nous mène au théorème suivant.

Soit f une fonction continue presque partout, et g une fonction dont la dérivée est continue presque partout sur un intervalle [a,b]. Soit également l’image de g contenue dans le domaine de définition de f. Alors \int_a^b f(g(x))g'(x){\mathrm{d}}x = \int_{g(a)}^{g(b)}f(z){\mathrm{d}}z.

Nous utilisons ce théorème de la façon suivante. L’idée est de remplacer la fonction g(x) par z. Puis il faut également remplacer {\mathrm{d}}x par {\mathrm{d}}z où nous avons que {\mathrm{d}}x={\mathrm{d}}z/g'(x). Finalement, il faut changer les bornes d’intégration par a\rightarrow g(a) et b\rightarrow g(b). Si on ne calcule pas l’intégrale mais la primitive, on ne modifie (évidemment) pas les bornes d’intégration, mais en revanche pour trouver la primitive il faut également appliquer la transformation x=g^{-1}(z) sur la solution.

Intégrer par changement de variables \int_1^3 6x\ln(x^2){\mathrm{d}}x.

En définissant z=x^2, nous avons {\mathrm{d}}x={\mathrm{d}}z/(2x). Les bornes d’intégration deviennent z(1)=1^2=1 et z(3)=3^2=9. On obtient donc \begin{aligned} \int_1^3 6x\ln(x^2){\mathrm{d}}x&=\int_1^9 6x\ln(z)\frac{1}{2x}{\mathrm{d}}z=\int_1^9\ln(z){\mathrm{d}}z\nonumber\\ &=3\left[z\ln(z)-z\right]_1^9=3(9\ln(9)-9-\ln(1)+1)=27\ln(9)-24. \end{aligned}

Calculer les primitives suivantes par changement de variable

  1. \int \frac{1}{5x-7}{\mathrm{d}}x

  2. \int \sin(3-7x){\mathrm{d}}x

  3. \int x e^{x^2}{\mathrm{d}}x

Intégration numérique

Dans certains cas, il est impossible d’évaluer analytiquement une intégrale ou alors elle est très compliquée à calculer. Dans ce cas, on va approximer l’intégrale et donc commettre une erreur.

Pour ce faire on va subdiviser l’espace d’intégration [a,b] en N pas équidistants (pour simplifier) \delta x=(b-a)/N, et approximer l’intégrale par une somme finie \int_a^bf(x){\mathrm{d}}x=\sum_{i=0}^{N-1} \delta x f(a+i\delta x) g_i+E(a,b,\delta x)\cong\sum_{i=0}^{N-1} \delta x f(a+i\delta x) g_i,g_i est un coefficient qui va dépendre de la méthode d’intégration que nous allons utiliser, E est l’erreur commise par l’intégration numérique et va dépendre des bornes d’intégration, de \delta x (du nombre de pas d’intégration), de la forme de f(x) (combien est “gentille”) et finalement de la méthode d’intégration.

Erreur d’une méthode d’intégration

D’une façon générale plus \delta x est petit (N est grand) plus l’erreur sera petite et donc l’intégration sera précise (et plus le calcul sera long). Néanmoins, comme la précision des machines sur lesquelles nous évaluons les intégrale est finie, si \delta x devient proche de la précision de la machine des erreurs d’arrondi vont dégrader dramatiquement la précision de l’intégration.

De façon générale il est difficile de connaître à l’avance la valeur exacte de E. En revanche on est en capable de déterminer l’ordre de l’erreur.

On dit qu’une méthode d’intégration est d’ordre k, si l’erreur commise par la méthode varie proportionnellement à \delta x^k. On note qu’une erreur est d’ordre k par le symbole \mathcal{O}(\delta x^k). Exemple: si une méthode est d’ordre deux, alors en diminuant \delta x d’un facteur 2, l’erreur sera elle divisée par 2^2=4. Si une méthode est d’ordre 3, alors en diminuant \delta x d’un facteur 2, nous aurons que l’erreur est divisée par un facteur 2^3=8. Etc.

Comme le calcul d’une intégrale de façon numérique ne donne en général pas un résultat exact, mais un résultat qui va dépendre d’un certain nombre de paramètres utilisés pour l’intégration, il faut définir un critère qui va nous dire si notre intégrale est calculée avec une précision suffisante.

Si nous notons I(N,a,b,f,g) l’approximation du calcul de l’intégrale entre a et b de la fonction f avec une résolution N pour la méthode d’intégration g I(N,a,b,f,g)=\sum_{i=0}^{N-1} \delta x f(a+i\delta x) g_i,g_i est encore à préciser. Afin de déterminer si le nombre de points que nous avons choisi est suffisant, après avoir évalué I(N,a,b,f,g), nous évaluons I(2\cdot N,a,b,f,g). En d’autres termes nous évaluons l’intégrales de la même fonction avec la même méthode mais avec un nombre de points deux fois plus élevé. Puis, nous pouvons définir \varepsilon(N) comme étant l’erreur relative de notre intégration avec une résolution N et 2\cdot N \varepsilon(N)\equiv\left|\frac{I(2N)-I(N)}{I(2N)}\right|. Si à présent nous choisissons un \varepsilon_0>0 (mais plus grand que la précision machine), nous pouvons dire que le calcul numérique de notre intégrale a convergé (on parle de convergence du calcul également) pour une résolution N quand \varepsilon(N)<\varepsilon_0.

Méthode des rectangles

Pour la méthode des rectangles, nous allons calculer l’intégrale en approximant l’aire sous la fonction par une somme de rectangles, comme nous l’avons fait pour la définition de l’intégration au sens de Riemann. La différence principale est que nous ne regarderons pas les valeurs minimales ou maximales de f sur les subdivisions de l’espace, mais uniquement les valeurs sur les bornes. Cette approximation donne donc la formule suivante \begin{aligned} \int_a^bf(x){\mathrm{d}}x&\cong\sum_{i=0}^{N-1} \delta x f(a+i\cdot\delta x)+\mathcal{O}(\delta x),\\ &\cong\sum_{i=1}^{N} \delta x f(a+i\cdot\delta x)+\mathcal{O}(\delta x)\end{aligned}{#eq:rect_gauche} Cette méthode est d’ordre 1. Une exception s’applique cependant concernant l’ordre de l’intégration. Si la fonction à intégrer est une constante f(x)=c, alors l’intégration est exacte.

Dans les deux cas ci-dessus on a évalué la fonction sur une des bornes. On peut améliorer la précision en utilisant le “point du milieu” pour évaluer l’aire du rectangle. L’approximation devient alors \begin{aligned} \int_a^bf(x){\mathrm{d}}x&\cong\sum_{i=0}^{N-1} \delta x f(a+(i+1/2)\cdot\delta x)+\mathcal{O}(\delta x^2).\end{aligned} Cette astuce permet d’améliorer la précision de la méthode à très faible coût. En effet, la précision de la méthode des rectangles est améliorée et devient d’ordre 2.

Méthode des trapèzes

Pour la méthode des trapèzes, nous allons calculer l’intégrale en approximant l’aire sous la fonction par une somme de trapèzes. Pour rappel l’aire d’un trapèze, dont les côtés parallèles sont de longueurs c et d et la hauteur h, est donnée pas A=(c+d)h/2. Cette approximation donne donc la formule suivante \int_a^bf(x){\mathrm{d}}x\cong\sum_{i=0}^{N-1} \delta x \frac{f(a+i\cdot\delta x)+f(a+(i+1)\cdot\delta x)}{2}+\mathcal{O}(\delta x^2). Cette méthode est d’ordre 2. Cette méthode d’intégration est exacte pour les fonctions linéaires f(x)=c\cdot x + d.

Méthode de Simpson

Pour cette approximation, on approxime la fonction à intégrer dans un intervalle par une parabole.

Commençons par évaluer l’intégrale à l’aide d’une subdivision dans l’ensemble [a,b].

L’idée est la suivante. On pose f(x)=c\cdot x^2+d\cdot x+e. Donc, il nous faut déterminer c, d, et e. Il nous faut donc choisir 3 points dans l’intervalle [a,b] pour déterminer ces constantes. On choisit comme précédemment f(a), f(b), et le troisième point est pris comme étant le point du milieu (f(a+b)/2). On se retrouve donc avec trois équations à trois inconnues \begin{aligned} f(a)&=c\cdot a^2+d\cdot a+e,\\ f(b)&=c\cdot b^2+d\cdot b+e,\\ f((a+b)/2)&=\frac{c}{4}\cdot (a+b)^2+\frac{d}{2}\cdot (a+b)+e.\end{aligned} En résolvant ce système (nous n’écrivons pas la solution ici) nous pouvons à présent évaluer l’intégrale \begin{aligned} I&=\int_a^b f(x){\mathrm{d}}x\cong\int_a^b (cx^2+dx+e){\mathrm{d}}x,\nonumber\\ &=\frac{b-a}{6}(f(a)+f(b)+4f((a+b)/2))+\mathcal{O}(\delta x^4).\end{aligned}

On peut donc généraliser affiner cette formule en rajoutant des intervalles comme précédemment et en répétant cette opération pour chaque intervalle.

Il vient donc que \begin{aligned} I&=\frac{\delta x}{6}\sum_{i=0}^{N-1}\left[f(a+i\cdot \delta x)+f(a+(i+1)\cdot\delta x)\right.\nonumber\\ &\left.+4f(a+(i+1/2)\cdot\delta x)\right]+\mathcal{O}(\delta x^4).\end{aligned}

Cette méthode permet d’évaluer exactement des polynômes d’ordre 4, f(x)=ax^3+bx^2+cx+d.

Équations différentielles ordinaires

Introduction

Pour illustrer le concept d’équations différentielles, nous allons considérer pour commencer des systèmes qui évoluent dans le temps (évolution d’une population, taux d’intérêts, circuits électriques, ...).

Mouvement rectiligne uniforme

Imaginons que nous connaissons la fonction décrivant le vitesse d’une particle au cours du temps et notons la v(t). Nous savons également que la vitesse d’une particule est relié à l’évolution au cours du temps de sa position. Cette dernière peut être notée, x(t). En particulier, nous avons que la vitesse n’est rien d’autre que la dérivée de la position. On peut onc écrire une équation reliant la vitesse à la position x'(t)=v(t). Cette équation est appelée équation différentielle, car elle fait intervernir non seulement les fonctions x(t) et v(t), mais également la dérivée de la fonction x(t). Si maintenant nous précisons ce que vaut la fonction v(t) nous pourrons résoudre cette équation. Comme le nom de la sous-section le laisse entendre, nous nous intéressons à un mouvement rectiligne uniforme, qui a la particularité de décrire le mouvement d’un objet qui se déplace à vitesse constante. On a donc v(t)=v. Nous cherchons donc à résoudre l’équation différentielle x'(t)=v. Ou en d’autres termes, nous cherchons la fonction dont la dérivée donne C3. Vous savez sans doute que l’ensemble de fonctions satisfaisant la contrainte précédente est x(t)=v\cdot t+B,B est une constante arbitraire. On a donc la solution générale de cette équation différentielle qui n’est pas unique, mais qui donne une infinité de solution (comme quand nous avons calculé la primitive d’une fonction au chapitre précédent). Afin de trouver une solution unique, nous devons imposer une “condition intiale” à notre équation différentielle. En effet, si nous imposons la condition intiale x(t_0)=x_0, il vient x(t_0)=x_0=v\cdot t_0+B \Leftrightarrow B=x_0-v\cdot t_0. Finalement, la solution de l’équation différentielle est donnée par x(t)=v\cdot (t-t_0)+x_0.

La solution de l’équation différentielle x'(t)=v,\ x(t_0)=x_0, revient à calculer \begin{aligned} \int x'(t){\mathrm{d}}t=\int v {\mathrm{d}}t,\\ x(t)=v\cdot t + B.\end{aligned}

Mouvement rectiligne uniformément accéléré

Dans le cas du mouvement rectiligne d’un objet dont on le connaît que l’accélération, a(t), on peut également écrire une équation différentielle qui décrirait l’évolution de la position de l’objet en fonction du temps. En effet, l’accélération d’un objet est la deuxième dérivée de la position, soit x''(t)=a(t), ou encore la première dérivée de la vitesse. \begin{aligned} v'(t)&=a(t),\\ x'(t)&=v(t).\end{aligned}

Par simplicité supposons que l’accélération est constante, a(t)=a. On doit donc résoudre4 x''(t)=a, ou \begin{aligned} v'(t)&=a,\\ x'(t)&=v(t).\end{aligned}{#eq:xpv} Commençons pas le système d’équations ci-dessus. On commence par résoudre la première équation pour v(t) et on a v(t)=a\cdot t+C. En substituant ce résultat dans l’@eq:xpv, on a x'(t)=a\cdot t+C. On peut donc directement intégrer des deux côtés comme vu dans la sous-section précédente \begin{aligned} \int x'(t){\mathrm{d}}t&=\int a\cdot t+C{\mathrm{d}}t,\nonumber\\ x(t)&=\frac{a}{2}\cdot t^2+C\cdot t + D.\end{aligned} On a donc que la position d’un objet en mouvement rectiligne uniformément accéléré est donné par une parabole. Cette équation ci-dessus néanmoins a encore deux constante indéterminées. Pour les déterminer, on doit donc imposer deux conditions intiales. Une possibilité est d’imposer une condition initiale par équation v(t_0)=v_0,\mbox{ et } x(t_0)=x_0. On obtient donc v(t_0)=v_0=a\cdot t_0+C \Leftrightarrow C=v_0-a\cdot t_0, et x(t_0)=x_0=\frac{a}{2}\cdot t_0^2+D \Leftrightarrow D=x_0-\frac{a}{2}\cdot t_0^2. Finalement la solution est donc x(t)=\frac{a}{2}\cdot (t^2-t_0^2)+v_0\cdot (t-t_0)+x_0.

La solution de l’équation différentielle peut également se calculer de la façon suivante x''(t)=av,\ x(t_0)=x_0,\ v(t_0)=v_0. revient à calculer \begin{aligned} \int \int x''(t){\mathrm{d}}t{\mathrm{d}}t=\int \int a {\mathrm{d}}t{\mathrm{d}}t,\\ x(t)=\frac{a}{2}t^2+C\cdot t + D.\end{aligned}

Évolution d’une population

Imaginons une colonie de bactéries dont nous connaissons le taux de reproduction r. Nous connaissons le nombre de ces bactéries au temps t, qui est donné par n(t). Nous souhaitons connaître la population au temps t+\delta t. On a donc n(t+\delta t)=n(t)+(r\delta t)\cdot n(t)=n(t)(1+r\delta t).{#eq:evolpop} Imaginons que le taux de reproduction r=1/3600 s^{-1}, que la population à un temps donné t_0 est de n(t_0)=1000, et qu’on veuille connaître la population après \delta t=1h=3600s. Il vient alors n(t_0+3600)=(1+1/3600 \cdot 3600)\cdot n(t_0)=2\cdot1000=2000. Imaginons maintenant que nous voulions calculer la population après \delta t=2h=7200s. Nous avons deux façons de faire. Soit nous utilisons le résultat précédent n(t_1)=2000 avec t_1=t_0+3600 et évaluons la population après une heure supplémentaire (\delta t_1=3600s) n(t_1+3600)=(1+1/3600 \cdot 3600)\cdot n(t_1)=2\cdot 2000=4000.{#eq:comp} Soit nous reprenons l’équation de départ (voir l'@eq:evolpop) et nous obtenons n(t_0+7200)=(1+1/3600 \cdot 7200)\cdot n(t_0)=3\cdot 1000=3000. On voit que ces deux résultats ne sont pas égaux. Effectuer deux itérations de notre algorithme discret avec un pas d’itération de \delta t, ne correspond pas à effectuer une seule itération avec un pas deux fois plus grand (2\delta t). Néanmoins cela devrait être le cas plus \delta t\rightarrow 0.

Pour nous en convaincre faisons l’exercice suivant. Reprenons l’@eq:comp que vous pouvons réécrire comme n(t_0+2\delta t)=n(t_1+\delta t)=(1+r\delta t) n(t_1)=(1+r \delta t)(1+r \delta t) n(t_0)=(1+r\delta t)^2 n(t_0). Si à présent nous comparons les résultats obtenus pour \delta t_1=2\delta t dans l’@eq:evolpop on a \begin{aligned} n_1&=(1+r\delta t)^2 n(t_0)=(1+2r\delta t+(r\delta t)^2) n(t_0),\\ n_2&=(1+2r\delta t) n(t_0).\end{aligned} On trouve donc finalement que n_2-n_1=(r\delta t)^2n(t_0). On a donc que la différence tend bien vers 0 quand \delta t tend vers 0.

Afin de voir plus en détail ce qu’il se passe lorsque \delta t\rightarrow 0, reprenons l’équation de départ (l'@eq:evolpop), divisons la par \delta t et arrangeons les termes. Il vient \frac{n(t+\delta t)-n(t)}{\delta t}=r\cdot n(t). En prenant la limite \delta t\rightarrow 0 on voit apparaître la dérivée dans le membre de gauche de l’équation ci-dessus \lim\limits_{\delta t\rightarrow 0} \frac{n(t+\delta t)-n(t)}{\delta t}=n'(t)=r\cdot n(t).{#eq:cont} On voit qu’on a construit ici une équation différentielle à partir d’un système discret.

Nous pouvons à présent résoudre l’équation différentielle ci-dessus en se souvenant que la fonction dont la dérivée est proportionnelle à la fonction de départ est l’exponentielle. Il vient n(t)=C\exp(r t),{#eq:sol_pop} où C est une constante. Il est en effet trivial de montrer que cette solution satisfait l’@eq:cont. On voit également qu’il nous manque une condition pour avoir l’unicité de la solution ci-dessus (on ne connaît toujours pas C). La constante peut-être obtenue à l’aide d’une condition initiale (correspondant au n(t_0) de tout à l’heure). Si n(t_0)=n_0, on a pour C n(t_0)=C\exp(r t_0)=n_0 \Leftrightarrow C=n_0\exp(-r t_0). En substituant cette relation dans l'@eq:sol_pop, on obtient n(t)=n_0\exp(r (t-t_0)).

Autres illustrations de l’utilisation des équations différentielles

La plupart des systèmes naturels (ou moins naturels) peuvent être décrits à l’aide d’équations différentielles. Nous allons en écrire quelques exemples ci-dessous.

Systèmes proies-prédateurs

Considérons un système où nous avons des prédateurs (des guépards) et des proies (des antilopes)5. Supposons que les antilopes se reproduisent exponentiellement vite et que leur seul moyen de mourir est de se faire manger par les guépards et que la chance de se faire manger est proportionnelle au nombre de guépards. Les guépards meurent exponentiellement vite de faim et se reproduisent proportionnellement au nombre d’antilopes se trouvant dans le système.

Avec ces hypothèses, on peut écrire le système d’équations suivant (a est le nombre d’antilopes, et g le nombre de guépards) \begin{aligned} \frac{{\mathrm{d}}a}{{\mathrm{d}}t}&= \underbrace{k_a a(t)}_{(1)}-\underbrace{k_{g,a}g(t) a(t)}_{(2)},\\ \frac{{\mathrm{d}}g}{{\mathrm{d}}t}&= -\underbrace{k_g g(t)}_{(3)} +\underbrace{k_{a,g} a(t)g(t)}_{(4)}\end{aligned} Le terme (1) représente la reproduction des antilopes avec taux k_a. Le terme (2) représente la mort des antilopes qui se font manger par les guépards avec un taux k_{g,a} (la chance qu’un guépard rencontre une antilope). Le terme (3) est la mort des guépards avec un taux k_g. Finalement le terme (4) est la reproduction des guépards proportionnelle au nombre d’antilopes avec un taux k_{a,g}.

Nous avons à faire ici à un système d’équations différentielles. Nous n’allons pas nous intéresser à la résolution de ce système mais simplement étudier la solution à ce problème (voir la @fig:lkA et @fig:lkB).

L’évolution au cours du temps de la population d’antilopes et de guépards. Représentation paramétrique de l’évolution population d’antilopes et de guépards.

Deux représentation du système de Lotka--Volterra.

Circuits électriques: le circuit RC

Supposons que nous ayons le circuit RC de la Fig. @fig:rc, où nous avons une résistance (de résistance R) branchée en série avec une capacité (de capacité électrique C). Sur ce circuit nous avons une source qui délivre une tension U. Nous avons également un interrupteur qui quand il est en position (a) relie le circuit RC à la source, ce qui a pour effet de chargé la capacité. En position (b) la capacité se décharge et son énergie est dissipée dans la résistance.

Le circuit RC.

Nous souhaitons étudier la variation de la chute de tension dans la capacité U_c lorsque:

  1. nous mettons l’interrupteur en position (a).

  2. puis lorsque la capacité est chargée, nous mettons l’interrupteur en position (b).

Les chutes de tension dans la capacité et la résistance sont respectivement données par U_C=Q/C,\quad U_R=R I,Q est la charge de la capacité et I le courant traversant la résistance. Nous avons par la loi de Kirchoff que U=U_C+U_R.{#eq:tot_tension} De plus le courant traversant la résistance est donné par I(t)=Q'(t). En combinant ces équations, nous obtenons U_C'(t)+\frac{U_C(t)}{RC}=\frac{U}{RC}. Nous avons également la condition initiale U_C(0)=0 (la tension au moment de la mise de l’interrupteur en position (a) est nulle).

Lors de la mise de l’interrupteur en position (b) nous avons simplement que l'@eq:tot_tension devient 0=U_C+U_R.{#eq:tot_tension_0} On a donc que l’équation différentielle pour l’évolution de la chute de tension dans la capacité devient U_C'(t)+\frac{U_C(t)}{RC}=0. Et la condition initiale devient U_C(0)=U.

Pour cette dernière équation nous avons déjà calculé une solution très similaire et on a U_C(t)=U\exp(-t/(RC)). La tension dans la capacité va décroître exponentiellement vite. Pour le cas de l’interrupteur en position (a) la solution est U_C(t)=U(1-\exp(-t/(RC))). La tension augmente exponentiellement au début, puis au fur et à mesure que la capacité se charge il devient de plus en plus difficile de la charger. L’augmentation de la tension se fait donc de plus en plus lentement jusqu’à ce qu’on tende vers une asymptote horizontale en U.

Taux d’intérêts composés

Nous voulons étudier l’augmentation d’un capital c(t) au cours du temps qui est soumis à un taux d’intérêt annuel r qui est composé après chaque intervalle \delta t. On peut également inclure des dépôts/retraits d sur l’intervalle \delta t. La valeur du capital après un intervalle \delta t est de c(t+\delta t)=c(t)+(r\delta t )c(t)+d\delta t.{#eq:cap_discr} Supposons qu’on a un capital de départ 1000 \mathrm{CHF}, un taux d’intérêts annuel de 1\% et un dépôt annuel de 100\mathrm{CHF}. Après deux mois (\delta t=2/12=1/6) on a donc que le capital devient c(1/6)=1000+0.01/6\cdot 1000 +100/6=1018.3\mathrm{CHF}. Si maintenant, nous voulons avoir la valeur du capital à n’importe quel moment dans le temps, nous allons prendre \delta t\rightarrow 0. En divisant l'@eq:cap_discr par \delta t, et en réarrangeant les termes, on obtient c'(t)=rc(t)+d. En supposant que c(t=0)=c_0 (le capital initial), cette équation différentielle a pour solution c(t)=\frac{d}{r}(e^{rt}-1)+c_0e^{r t}. Cette solution a pour les paramètres précédent la forme suivante sur une période de 100 ans.

L’évolution du capital c en fonction du temps su 100 ans.

Définitions et théorèmes principaux

Soit y une fonction dérivable n fois et dépendant d’une seule variable. Une équation différentielle ordinaire est un équation de la forme F(x,y,y',y'',...,y^{(n)})=0,F est une fonction, et y', y'', ..., y^{(n)} sont les dérivées premières, deuxièmes, ..., n-ème de y.

L’équation suivante est une équation différentielle ordinaire y''+4y'+8y+3x^2+9=0.

Afin de résoudre cette équation, nous cherchons une solution de la forme y=f(x). On dit également que nous cherchons à intégrer l’équation différentielle.

Afin de classifier les équation différentielles, considérons les deux définitions suivantes

L’ordre d’une équation différentielle est l’ordre le plus haut des dérivées de y qui y apparaissent. L’ordre de l’équation différentielle F(x,y,y',y'',...,y^{(n)})=0 est de n, si n\neq 0.

L’équation différentielle suivante est d’ordre 3 4y'''+x\cdot y'+4y+6x=0.

Une condition initiale pour une équation différentielle d’ordre n, est un ensemble de valeurs, y_0, y_1, ..., y_{n-1} donnée telles que pour une valeur x_0 donnée on a y(x_0)=y_0,\ y'(x_0)=y_1,\ ...,\ y^{(n-1)}=y_{n-1}.

Nous souhaitons maintenant savoir sous quelles conditions une équation différentielle admet une solution et si elle est unique. Nous n’allons pas vraiment écrire ni démontrer le théorème d’existence et d’unicité des équations différentielles ordinaires, mais simplement en donner une version approximative et la discuter

Soit D\subseteq{\mathbb{R}} le domaine de définition de la fonction y. Soit y:D\rightarrow E\subseteq {\mathbb{R}} une fonction à valeur réelle continue et dérivable sur D, et f:D\times E\rightarrow F\subseteq{\mathbb{R}} une fonction continue sur D\times E. Alors, le système suivant (également appelé problème de Cauchy) \begin{aligned} &y'=f(y,x),\\ &y(x=x_0)=y_0, \end{aligned} admet une unique unique solution y(x).

Ce théorème peut être étendu à une équation d’un ordre arbitraire n possédant n-1 conditions initiales. En effet, n’importe quel équation différentielle d’un ordre n peut être réécrite sous la forme de n équations différentielles d’ordre 1. Pour illustrer cette propriété considérons l’équation différentielle suivante y''+3y'+y+3x=0. Si nous définissons z=y', nous avons le système suivant à résoudre \begin{aligned} y'=z,\\ z'+3y'+y+3x=0.\end{aligned} Nous voyons que ce système est d’ordre 1, mais que nous avons augmenté le nombre d’équations à résoudre.

Cette propriété peut se généraliser de la façon suivante. Soit une équation différentielle d’ordre n F(x,y,y',...,y^{(n)})=0. Nous pouvons définir z_i=y^{(i-1)} et on aura donc que z_{i+1}=z_i'. On peut donc réécrire l’équation différentielle d’ordre n comme étant \begin{aligned} &z_{i+1}=z_i',\ i=1,...,n-1\\ F(x,y,y',..,y^{(n)})=0 \Rightarrow &G(x,z_1,z_2,...,z_n)=0.\end{aligned}

Jusqu’ici F peut être totalement arbitraire. Essayons de classifier un peu les équations différentielles en fonction des propriétés du F.

Une équation différentielle ordinaire d’ordre n est dite linéaire si on peut l’écrire sous la forme a_0(x)\cdot y(x)+a_1(x)\cdot y'(x)+...+a_n(x)\cdot y^{(n)}(x)=b(x). Si les coefficients a_i ne dépendent pas de x, alors l’équation est dite à coefficients constants.

L’équation ci-dessus a les deux propriétés suivantes

  1. Les a_i ne dépendent que de x (ils ne peuvent pas dépendre de y).

  2. Les y et toutes leur dérivées ont un degré 1.

L’équation suivante est linéaire y''+4x\cdot y'=e^x. L’équation suivante n’est pas linéaire y\cdot y''+4x\cdot y'=e^x.

Une équation différentielle ordinaire est dite homogène si le terme dépendant uniquement de x est nul. Dans le cas où nous avons à faire à une équation différentielle linéaire, cela revient à dire que b(x)=0.

Les équations suivantes sont homogènes \begin{aligned} &y''+4x\cdot y\cdot y'+3x^2\cdot y^3=0,\\ &2y'''+5x^2\cdot y'=0. \end{aligned} Les équations suivantes ne le sont pas \begin{aligned} &y''+4x\cdot y\cdot y'+3x^2\cdot y^3=4x+2,\\ &2y'''+5x^2\cdot y'=1. \end{aligned}

Pour chacune de ces équations différentielles ordinaires suivantes donner tous les qualificatifs possibles. Si l’équation est inhomogène donner l’équation homogène associée. \begin{aligned} &y^{(4)}+4x^2 y=0,\\ &y'+4x^2 y^2=3x+2,\\ &\frac{1}{y+1}y''+4x^2 y^2=0,\\ &y'=y,\\ &4y''+4x y=1. \end{aligned}

Lors de la résolution d’équation différence inhomogène la solution se trouve de la façon suivante.

  1. Trouver la solution de l’équation différentielle homogène associée, notons-la y_h(x).

  2. Trouver une solution particulière à l’équation inhomogène, notons-la y_0(x).

La solution sera donnée par la somme de ces deux solutions y=y_h+y_0.

Techniques de résolution d’équations différentielles ordinaires d’ordre 1

Ici nous considérerons uniquement les équations différentielles ordinaires d’ordre 1. Pour certains types d’équations différentielles, il existe des techniques standard pour les résoudre. Nous allons en voir un certain nombre.

Équations à variables séparables

On dit qu’une équation différentielle d’ordre 1 est à variables séparables, si elle peut s’écrire sous la forme suivante y' a(y)=b(x).

L’équation suivante est à variables séparables e^{x^2+y^2(x)}y'(x)=1.

Pour ce genre d’équations, la solution se trouve de la façon suivante. Nous commençons par écrire la dérivée, y'={\mathrm{d}}y/{\mathrm{d}}x et on obtient \begin{aligned} \frac{{\mathrm{d}}y}{{\mathrm{d}}x} a(y)=b(x),\\ a(y){\mathrm{d}}y=b(x){\mathrm{d}}x.\end{aligned} On peut donc simplement intégrer des deux côtés et on obtient \int a(y){\mathrm{d}}y=\int b(x){\mathrm{d}}x. Si nous parvenons à résoudre les intégrales nous obtenons une solution pour y(x) (cette solution n’est peut-être pas explicite). Il existe le cas simple où a(y)=1 et il vient y=\int b(x){\mathrm{d}}x.

Résoudre l’équation différentielle suivante n'(t)=r\cdot n(t). En écrivant n'={\mathrm{d}}n /{\mathrm{d}}t, on réécrit l’équation différentielle sous la forme \frac{1}{n} {\mathrm{d}}n=r{\mathrm{d}}t, qu’on intègre \begin{aligned} \int \frac{1}{n} {\mathrm{d}}n&=\int r{\mathrm{d}}t,\nonumber\\ \ln(n)&=r\cdot t+C,\nonumber\\ n(t)&=e^{r\cdot t+C}=A\cdot e^{r\cdot t},\end{aligned}A=e^C.

  1. Résoudre l’équation différentielle suivante c'(t)=rc(t)+d.

  2. Résoudre l’équation différentielle suivante x\cdot y(x) \cdot y'(x)=1.

Équations linéaires {#sec:eq_lin}

Pour une équation du type y'(x)=a(x)\cdot y(x)+b(x),{#eq:lin} on doit résoudre le problème en deux parties.

Pour résoudre ce genre d’équation supposons que nous connaissons une solution “particulière” à cette équation. Notons la y_p. Si nous faisons maintenant le changement de variables y=y_h+y_p et remplaçons ce changement de variables dans l’équation ci-dessus y_p'(x)+y_h'(x)=a(x)\cdot y_p(x)+a(x)\cdot y_h(x)+b(x).{#eq:lin_hp} Comme y_p est solution de l'@eq:lin on a y_p'(x)=a(x)\cdot y_p(x)+b(x). En remplaçant cette relation dans l'@eq:lin_hp il vient y_h'(x)=a(x)\cdot y_h(x). Cette équation différentielle n’est rien d’autre que l’équation homogène correspondant à @eq:lin.

Nous voyons donc qu’une équation inhomogène se résout en trouvant la solution générale à l’équation homogène correspondante et en y ajoutant une solution particulière.

Revenons donc à la résolution de l’équation différentielle linéaire d’ordre un. La première partie de la résolution consiste à résoudre l’équation homogène associée à l'@eq:lin y'(x)=a(x)\cdot y(x). Cette équation se résout par séparation des variables. La solution est donc y_h(x)=Ce^{\int a(x){\mathrm{d}}x}. Puis nous devons chercher une solution dite particulière de l’équation inhomogène. Pour ce faire nous utilisons la méthode de la variation de la constante. Il s’agit de trouver une solution particulière qui aura la même forme que la solution de l’équation homogène, où C dépendra de x (méthode de variation de la constante) y_p(x)=C(x)e^{\int a(x){\mathrm{d}}x}. En remplaçant cette équation dans l'@eq:lin, on obtient \begin{aligned} C'(x)e^{\int a(x){\mathrm{d}}x}+C(x)\cdot a(x)e^{\int a(x){\mathrm{d}}x}&=a(x)\cdot C(x) e^{\int a(x){\mathrm{d}}x}+b(x),\nonumber\\ C'(x)&=\frac{b(x)}{e^{\int a(x){\mathrm{d}}x}}. \end{aligned} Il nous reste donc à résoudre cette équation différentielle pour C(x) qui est une équation à variable séparable où on aurait un a(c)=1. On intègre donc directement cette équation et on obtient C(x)=\int \frac{b(x)}{e^{\int a(x){\mathrm{d}}x}}{\mathrm{d}}x. Finalement, on a que la solution de l’équation générale de l’équation inhomogène est y=y_p+y_h=\left(\int \frac{b(x)}{e^{\int a(x){\mathrm{d}}x}}{\mathrm{d}}x+C\right)e^{\int a(x){\mathrm{d}}x}.

Résoudre l’équation suivante U_C'(t)+\frac{U_C(t)}{RC}=\frac{U}{RC}.{#eq:rc_inhom} On commence par résoudre l’équation homogène {U_C}_h'(t)+\frac{{U_C}_h(t)}{RC}=0. D’où on obtient {U_C}_h=A\cdot e^{-\frac{1}{RC} t}. Puis par variations des constantes, on essaie de déterminer la solution particulière de la forme {U_C}_p=B(t)\cdot e^{-\frac{1}{RC} t}. En remplaçant cette forme de solution dans l'@eq:rc_inhom, on obtient B'(t)=\frac{U}{RC}\cdot e^{\frac{1}{RC} t}. Qui donne par intégration B(t)=U e^{\frac{1}{RC} t}+D. Finalement, il vient que U_c(t)=\left(U e^{\frac{1}{RC} t}+D+A\right)e^{-\frac{1}{RC}t}=U+(D+A)e^{-\frac{1}{RC}t}=U+Ce^{-\frac{1}{RC}t},C=D+A. Pour le cas de la charge du condensateur, on a de plus U_c(0)=0. On peut donc fixer la constante C=-U.

Résoudre les équations différentielles suivantes

  1. y'+2y=t^2

  2. y'+y=\frac{1}{1+e^t}.

Équations de Bernouilli

Il existe des équations particulières qui peuvent se ramener à des équations linéaires via des changements de variables.

Une classe particulière sont les équations de Bernouilli, qui s’écrit y'(x)+a(x)\cdot y(x)+b(x)\cdot y^n(x)=0,{#eq:bernouilli} où r\in{\mathbb{R}}.

Cette équation peut également être réécrite sous la forme \frac{y'(x)}{y^n(x)}+\frac{a(x)}{y^{n-1}(x)}+b(x)=0.{#eq:bernouilli_2}

Dans ce cas là, en effectuant le changement de variable suivant z=y^{1-n}, on obtient (exercice) z'(x)+(1-n)a(x)\cdot z(x)+(1-n)b(x)=0. On a donc ramené l’équation de Bernouilli à une équation linéaire que nous savons résoudre à l’aide de la méthode de la section @sec:eq_lin.

Résoudre l’équation de Bernouilli suivante y'-y-x\cdot y^6=0. Avec la substitution z=y^5, on obtient z'-5z+5x=0. Cette équation se résout avec se résout en trouvant d’abord la solution de l’équation homogène z_h'-5z_h=0, qui est donnée par z_h=Ae^{5x}. En remarquant qu’une solution particulière à z_p'-5z_p+5x=0, peut être de la forme z_p=x+B (avec B une constante) on obtient \begin{aligned} 1-5(x+B)+5x=0,\nonumber\\ 1-5B=0\Rightarrow B=\frac{1}{5}.\end{aligned} Et finalement z=z_h+z_p=Ae^{5x}+x+\frac{1}{5}. Il nous reste à présent à calculer y=z^{1/5} et on a y=\left(Ae^{5x}+x+\frac{1}{5}\right)^{1/5}.

Équation de Riccati

L’équation de Riccati qui est de la forme y'(x)+a(x)+b(x)\cdot y(x)+c(x)\cdot y^2(x)=0,{#eq:riccati} et est donc quadratique en y. On notera que c’est une équation de Bernouilli (avec n=2 et qui est inhomogène).

Cette équation a une propriété intéressante. Si nous connaissons une solution particulière à l’équation inhomogène, notons la y_p, alors la solution générale peut être trouvée de la façon suivante.

Faisons le changement de variable suivant y=y_h+y_p. L’équation ce-dessus devient donc y_p'+y_h'+a(x)+b(x)\cdot y_p+b(x)\cdot y_h+c(x)\cdot (y_p^2+2y_p(x)y_h(x)+y_h^2)=0. En utilisant que y_p est solution de l’équation de Riccati, on a y_h'+a(x)+(b(x)+2y_p(x)c(x))\cdot y_h+c(x)\cdot y_h^2=0. Cette équation est une équation de Bernouilli avec n=2. On sait donc comment la résoudre.

Résoudre l’équation de Riccati suivante y'+y^2-\frac{2}{x^2}=0. Indication: la solution particulière a la forme y=\frac{a}{x}, avec a une constante.

De plus, ce genre d’équation peut-être transformée via un changement de variables en une équation linéaire d’ordre deux. Si c(x) est dérivable, alors on peut faire le changement de variables suivant v=y\cdot c(x), et on a donc que v'=y' c+y c'. En insérant ces relations dans l'@eq:riccati, il vient v'(x)+d(x)+e(x)\cdot v(x)+v^2(x)=0,{#eq:riccati_2} où nous avons nommé d(x)=a(x)\cdot c(x) et e(x)=\frac{c'(x)}{c(x)}+b(x). Si à présent nous faisons un autre changement de variables v(x)=-\frac{z'(x)}{z(x)}, on obtient que l’équation ci-dessus peut se réécrire comme z''(x)+e(x)\cdot z'(x)+d(x)\cdot z(x)=0.{#eq:riccati_3} L’équation de Riccati (une équation d’ordre un non-linéaire et inhomogène) est donc transformée en une équation linéaire d’ordre deux.

Equations différentielles ordinaires d’ordre deux

Dans cette section, nous allons étudier des cas particuliers d’équations différentielles que nous savons résoudre. Cela sera toujours des équations linéaires.

De façon générale ces équations s’écrivent a(x)y''(x)+b(x)y'(x)+c(x)y(x)=d(x),a,b,c,d:{\mathbb{R}}\rightarrow{\mathbb{R}} sont des fonctions réelles. Avant de résoudre l’équation générale, nous allons considérer des plus simples.

EDO d’ordre deux à coefficients constants homogènes

Ce genre d’équations s’écrit sous la forme a y''(x)+by'(x)+cy(x)=0.{#eq:edo2_cch} Voyons maintenant comment résoudre cette équation.

Ces équations ont des propriétés intéressantes dûes à la linéarité de l’équation différentielle.

Ces propriétés sont à démontrer en exercice.

  1. Soit f(x) une solution de l'@eq:edo2_cch, alors on a aussi que pour C\in{\mathbb{R}} Cf(x) est également solution de @eq:edo2_cch.

  2. Soient f(x) et g(x) deux solutions de l’équation @eq:edo2_cch, alors on a aussi que pour h(x)=f(x)+g(x) est également solution de @eq:edo2_cch.

  3. De ces deux propriétés, on déduit la propriété suivante. Soient f(x) et g(x) deux solutions de l'@eq:edo2_cch, et C_1,C_2\in{\mathbb{R}}, alors on a que h(x)=C_1f(x)+C_2g(x) sont solution de l'@eq:edo2_cch.

Afin de simplifier la discussion prenons une EDO d’ordre deux à coefficients constants particulière y''+3y'+2y=0.{#eq:edo2_ex} On va supposer que cette équation a pour solution une fonction de la forme y(x)=e^{\lambda x}. Substituons cette forme de solution dans l’équation de départ, on obtient \begin{aligned} \lambda^2 e^{\lambda x}+3\lambda e^{\lambda x}+2\lambda^2 e^{\lambda x}=0,\nonumber\\ \lambda^2+3\lambda +2=0,\end{aligned}s où on a utilisé que e^{\lambda x} ne peut jamais s’annuler pour le simplifier entre les deux lignes. La seconde ligne ci-dessus, s’appelle le polynôme caractéristique de notre EDO d’ordre 2.

Il nous reste à présent à déterminer \lambda ce qui est un simple problème d’algèbre. Le polynome ci-dessus se factorise simplement en (\lambda+1)(\lambda+2)=0, on a donc pour solution \lambda=-1, et \lambda=-2.

On a donc immédiatement deux solutions à notre équation différentielle y_1(x)=e^{-x},\quad y_2(x)=e^{-2x}. On vérifie aisément que ces deux équations vérifient l'@eq:edo2_ex. Précédemment, nous avons vu que la linéarité de ces équations différentielles, faisait qu’on pouvait contrsuire des solutions plus générales. En effet, on peut montrer que la solution la plus générale à cette EDO est y(x)=C_1 y_1(x)+C_2y_2(x)=C_1e^{-x}+C_2e^{-2x}. On constate qu’il y a deux constantes à déterminer pour avoir une solution unique. Pour ce faire il faudra donner deux conditions initiales. Une sur y(x) et une sur y'(x). Par exemple on pourrait avoir y(0)=1 et y'(0)=0 et on obtient \begin{aligned} C_1+C_2&=1,\\ -C_1-2C_2&=0.\end{aligned} Ce système d’équations ordinaires a pour solution C_1=2,\quad C_2=-1. On a donc finalement y(x)=2e^{-x}-e^{-2x}.

A présent, nous pouvons généraliser cette méthode pour l’équation @eq:edo2_cch a y''(x)+by'(x)+cy(x)=0. En faisans la même subsitution que précédemment, y=e^{\lambda x}, on a \begin{aligned} &a y''(x)+by'(x)+cy(x)=0,\\ &a \lambda^2+\lambda b+c=0.\end{aligned} L’équation ci-dessus doit être résolue pour \lambda. Nous savons comment résoudre ce genre d’équation du second degré. La solution est donnée par \lambda=\frac{-b\pm\sqrt{\Delta}}{2a},\Delta = b^2-4ac. On a donc deux solutions \begin{aligned} \lambda_1=\frac{-b-\sqrt{\Delta}}{2a},\\ \lambda_2=\frac{-b+\sqrt{\Delta}}{2a}.\end{aligned}

Il y a donc trois cas différents possibles: \Delta > 0, \Delta = 0, \Delta < 0.

Le cas \Delta>0

Dans ce cas, on a que \lambda_1,\lambda_2\in{\mathbb{R}} sont réels. La solution est donc donnée par (comme on l’a vu au paravent) y(x)=C_1e^{\lambda_1 x}+C_2e^{\lambda_2 x}.

Le cas \Delta=0

Dans ce cas, on a que \lambda_1=\lambda_2=\lambda=-b/(2a) et est réel. Dans ce cas-là les choses se compliquent un peu. Si on utilisait directement la formule ci-dessus, on aurait y(x)=Ce^{\lambda x}, avec C\in{\mathbb{R}}. Par contre, cette solution ne peut pas satisfaire deux conditions initiales comme nous avons vu précédemment. Il fau donc travailler un peu plus. Supposons que y(x) est donné par la fonction suivante y(x)=z(x)e^{\lambda x}, avec z(x) une fonction réelle. En substituant cela dans l’équation générale, on a az''+(2\lambda a+b)z'+(a\lambda^2+b\lambda+c)z=0. En utilant que \lambda=-b/(2a) et \Delta =0 il vient z''=0. La solution de cette équation est z=C_1+xC_2. On obtient donc comme solution générale de l’équation différentielle y(x)=(C_1+C_2 x)e^{\lambda x}.

Le cas \Delta<0

Dans ce cas-là, on a deux solutions complexes (la racine d’une nombre négative n’est pas réelle). Les racines sont de la forme \begin{aligned} \lambda_1=\frac{-b+i\sqrt{|b^2-4ac|}}{2a}, \lambda_2=\frac{-b-i\sqrt{|b^2-4ac|}}{2a},\end{aligned}i est le nombre imaginaire. En écrivant u=-b/(2a) et v=\sqrt{|b^2-4ac|}/(2a), on peut écrire \lambda_1=u+iv et \lambda_2=u-iv. On a donc que \lambda_2 est le complexe conjugué de \lambda_1, ou \lambda_1=\bar{\lambda}_2. En utilisant ces notations dans notre exponentielle, on a \begin{aligned} y_1&=e^{(u+iv)x}=e^{ux}e^{ivx},\\ y_2&=e^{(u-iv)x}=e^{ux}e^{-ivx}.\end{aligned} En se rappelant de la linéarité des solutions des EDO linéaires, on peut écrire la forme générale de la solution comme (C_1,C_2\in {\mathbb{R}}) y=C_1y_1+C_2y_2=C_1e^{ux}e^{ivx}+C_2e^{ux}e^{-ivx}=e^{ux}(C_1e^{ivx}+C_2e^{-ivx}).{#eq:sol2}

En utilisant la formule d’Euler \begin{aligned} e^{ivx}&=(\cos(vx)+i\sin(vx)),\\ e^{-ivx}&=e^{ux}(\cos(vx)-i\sin(vx)),\end{aligned} on peut réécrire l'@eq:sol2 comme \begin{aligned} y&=e^{ux}\left(C_1(\cos(vx)+i\sin(vx))+C_2(\cos(vx)-i\sin(vx))\right),\nonumber\\ &=e^{ux}\left((C_1+C_2)\cos(vx)+i(C_1-C_2)\sin(vx))\right),\nonumber\\ &=e^{ux}\left(C_3\cos(vx)+C_4\sin(vx))\right),\end{aligned} où on a définit C_3\equiv C_1+C_2 et C_4\equiv i(C_1-C_2).

Résoudre les EDO d’ordre 2 à coefficiens constants suivantes:

  1. y''+y'+y=0,

  2. y''+4y'+5y=0, y(0)=1, y'(0)=0.

  3. y''+5y'+6y=0, y(0)=2, y'(0)=3.

  4. 2y''-5y'+2y=0, y(0)=0, y'(0)=1.

Résolution numérique d’équations différentielles ordinaires

Pour la plupart des problèmes d’ingénierie classique, les solutions des équations différentielles sont trop compliquées à calculer analytiquement (si elles sont calculables). Il est donc nécessaire d’en obtenir des solutions approximées numériquement.

Problématique

Le problème à résoudre est une EDO avec condition initiale qui peut s’écrire de la façon suivante y'=F(t,y),\quad y(t_0)=y_0,F est une fonction de y et de t, et où y_0 est la condition initiale. Nous cherchons donc à connaître l’évolution de y(t) pour t>t_0.

Méthode de résolution: la méthode d’Euler

Afin de résoudre ce genre de problème numériquement il existe une grande quantité de techniques. Ici nous allons en considérer une relativement simple, afin d’illustrer la méthodologie (vous en verrez une autre dans le TP).

Nous cherchons donc à évaluer y(t=t_0+\delta t), étant donné y_0, \delta t et F(t,y). Intégrons donc simplement notre EDO entre t_0 et t_0+\delta t dans un premier temps et on obtient \int_{t_0}^{t_0+\delta t} y' {\mathrm{d}}t=\int_{t_0}^{t_0+\delta t} F(t,y){\mathrm{d}}t. Le théorème fondamental du calcul intégral nous dit que cette équation peut s’écrire y(t_0+\delta t)-y(t_0)=\int_{t_0}^{t_0+\delta t} F(t,y){\mathrm{d}}t, y(t_0+\delta t)-y_0=\int_{t_0}^{t_0+\delta t} F(t,y){\mathrm{d}}t.{#eq:edo_app_gen} Ont doit donc intégrer le membre de droite de cette équation. Pour ce faire nous pouvons utiliser une des techniques vues dans le chapitre précédent. Par exemple, on peut choisir la méthode des rectangle à gauche. Cette équation devient \begin{aligned} &y(t_0+\delta t)-y_0=\delta t F(t_0,y(t_0)),\nonumber\\ &y(t_0+\delta t)=y_0+\delta t F(t_0,y(t_0)).\end{aligned} Cette dernière équation nous permet donc d’évaluer y(t_0+\delta t) connaissant y_0. Cette méthode s’appelle “méthode d’Euler” et est une dite explicite, car y(t_0+\delta t) ne dépend que de la valeur de y évaluée au temps t_0.

Si plutôt que d’utiliser la méthode des rectangle à gauche pour approximer l’intégrale de l'@eq:edo_app_gen, nous utilisons la méthodes des rectangles à droite on a y(t_0+\delta t)=y_0+\delta t F(t_0+\delta t,y(t_0+\delta t)). Dans ce cas, on voit que la valeur y(t_0+\delta t) est calculée par rapport à la valeur d’elle même. Dépendant de la forme de F on ne peut pas résoudre cette équation explicitement. On a donc à faire à une équation sous forme implicite. Cette façon d’approximer une EDO est dite méthode d’Euler implicite.

Sans entrer dans les détails, la différence entre une méthode explicite et une méthode implicite est une question de stabilité numérique. En effet, les méthodes explicites peuvent devenir numériquement instables (la solution numérique s’éloigne exponentiellement vite de la solution de l’EDO) si \delta t devient “trop grand” (la contrainte du la taille de \delta t s’appelle CFL, pour Courant-Friedrich-Lévy). Les méthodes implicites ne souffrent pas de ce problème de stabilité, en revanche elles sont plus coûteuses en temps de calcul et en complexité algorithmique, étant donné qu’elles requièrent la résolution d’une équation implicite.

Notre but initial était de connaître l’évolution de y(t) pour t>t_0. Pour déterminer la valeur de y(t_1) avec t_N=t_0+N\delta t, il suffit donc d’effectuer N pas de la méthode d’intégration choisie (ici la méthode d’Euler explicite). On a donc que y(t_0+N\delta t)=y_0+\delta t\sum_{i=1}^{N}F(t_i,y_i),t_i=t_0+i\cdot\delta t et y_i=y(t_i). Le deuxième terme du membre de droite de cette équation est la même que la formule d’intégration en plusieurs pas pour la méthode du rectangle (voir l’équation @eq:rect_gauche). On a vu que cette méthode a une erreur d’ordre \delta t. On peut en conclure que l’erreur que la précision de la méthode d’Euler est également d’ordre \mathcal{O}(\delta t).

Méthode de résolution: la méthode de Verlet

Cette méthode d’intégration est utilisée pour l’intégration numérique d’EDO d’ordre deux avec une forme particulière qui est donnée par x''(t)=a(x(t)),{#eq:x2} où F est une fonction de x(t). On a également les conditions initiales x(t_0)=x_0 et x'(t_0)=v_0. Cette forme d’équation différentielle est bien connue en physique sous la forme \vec F=m\vec a, qui peut s’écrire \begin{aligned} &\vec{F}=m \vec a(t)=m \vec x''(t),\nonumber\\ &\frac{\vec{F}}{m}= \vec x''(t),\end{aligned} qui est de la forme de l’EDO de départ de l'@eq:x2. La force peut avoir différentes forme. Cela peut être la forme de gravité \vec F=m \vec g, de frottement \vec F=-\zeta \vec v=-\zeta x'(t), etc ou une combinaison de toutes ces forces.

Dans la section précédente, nous avons vu l’algorithme d’Euler pour résoudre des EDO. Cette méthode a pour avantage sa simplicité de codage, son faible coût de calcul, mais a pour désavantage son manque de précision. Dans un certain nombres d’applications, telles que les moteurs physiques pour les graphismes dans les jeux vidéos, ce manque de précision est inacceptable et une meilleure méthode doit être utilisée. Dans le TP vous avez vu les méthodes de Runge-Kutta. Ces méthodes améliorent la précision de façon spectaculaire, mais ont en général un coû de calcul trop élevé.

La méthode de Verlet qu’on va voir ci-dessous est augmente combine un faible coût de calcul et une amélioration notable de la précision. Elle est en effet très répandue dans l’industrie du jeu vidéo pour intégrer les équations différentielles omniprésentes dans les moteurs physiques.

La méthode de Verlet s’écrit (en utilisant les notations de la section précédente) x(t_{n+1})=x(t_n)+\delta t v(t_n)+\frac{1}{2}\delta t^2 a(x(t_n)).{#eq:verlet_gen} Considérons d’abord le terme v(t_n). Ce terme est approximé ici comme v(t_n) = \frac{x(t_{n+1})-x(t_{n-1})}{2\delta t}. En remplaçant cette approximation dans l’équation ci-dessus, il vient x(t_{n+1})=x(t_n)+\frac{x(t_{n+1})-x(t_{n-1})}{2}+\frac{1}{2}\delta t^2 a(x(t_n)), 2x(t_{n+1})=2x(t_n)+x(t_{n+1})-x(t_{n-1})+\delta t^2 a(x(t_n)), x(t_{n+1})=2x(t_n)-x(t_{n-1})+\delta t^2 a(x(t_n)).{#eq:verlet_novel}

On voit ici que cette formule est inutilisable pour évaluer x(t_1) (ce qui veut dire que n=0 dans le cas ce-dessus), car elle fait intervenir x(t_{-1}) dans le membre de droite. Pour résoudre ce problème il suffit d’évaluer x(t_1) grâce à l'@eq:verlet_gen où n=0 \begin{aligned} x(t_{1})&=x(t_0)+\delta t v(t_0)+\frac{1}{2}\delta t^2 a(x(t_0)),\nonumber\\ x(t_{1})&=x_0+\delta t v_0+\frac{1}{2}\delta t^2 a(x_0),\end{aligned}x_0 et v_0 sont les conditions initiales de notre problème. Esuite les itérations suivantes (n>0) sont calculables directement avec l'@eq:verlet_novel. Un autre avantage considérable de ce modèle est qu’il est très simple d’y inclure une force de frottement proportionnelle à la vitesse. Sans entrer dans les détails de la dérivation du schéma on a x(t_{n+1})=(2-\delta t\zeta)x(t_n)-(1-\delta t\zeta)x(t_{n-1})+\delta t^2 a(x(t_n)).

Transformées de Fourier

Rappel sur les nombres complexes

Dans cette section, on fait un rappel sur les nombres complexes qui seront beaucoup utilisés dans cette section.

Les nombres réels

L’ensemble des nombres réels, noté {\mathbb{R}}, possède un certain nombre de fonctions (opérateurs) tels que l’addition, la soustraction, la multiplication, la division, etc qui prennent un couple de nombres réels et rendent un autre nombre réel \begin{aligned} & +:{\mathbb{R}}\times{\mathbb{R}}\rightarrow{\mathbb{R}},\\ & \cdot:{\mathbb{R}}\times{\mathbb{R}}\rightarrow{\mathbb{R}},\\\end{aligned} On peut donc noter l’addition de deux nombres réels 7 et 2 et de la définition de l’addition on a +(7,2)=9. On lui préfère la notation +(7,2)=7+2=9. Intéressons nous plus particulièrement à la multiplication et à l’addition. Ces opérations ont les propriétés d’associativité et de commutativité. Cela veut dire que \begin{aligned} &(a+b)+c=a+(b+c), &(a\cdot b)\cdot c=a\cdot(b\cdot c),\\ &\quad\quad\quad\quad\quad\quad\mbox{ et }&\nonumber\\ &a+b=b+a,&a\cdot b=b\cdot a.\end{aligned}

Les couples de nombres réels

Intéressons-nous à présent à un ensemble plus grand que {\mathbb{R}}, soit {\mathbb{R}}^2\equiv{\mathbb{R}}\times{\mathbb{R}}. Cet ensemble est l’ensemble des des couples de nombres réels. Notons les nombres z\in{\mathbb{R}}^2 comme z=(a,b)\mbox{ tel que } a\in{\mathbb{R}}, \mbox{ et } b\in{\mathbb{R}}. Sur ces nombres on peut définir à nouveau l’addition, la soustraction, la multiplication, ... \begin{aligned} & +:{\mathbb{R}}^2\times{\mathbb{R}}^2\rightarrow{\mathbb{R}}^2,\\ & \cdot:{\mathbb{R}}^2\times{\mathbb{R}}^2\rightarrow{\mathbb{R}}^2.\end{aligned} On peut les écrire sous la forme de leurs équivalents des nombres réels sous la forme (a,b)+(c,d)=(a+c,b+d),{#eq:add} (a,b)\cdot(c,d)=(a\cdot c-b\cdot d,a\cdot d+b\cdot c).{#eq:mult} On voit assez facilement que l’addition sur {\mathbb{R}}^2 a une forme très similaire à celle sur {\mathbb{R}} du point de vue de ses propriétés telles que la commutativité ou l’associativité. Cela est moins clair pour la multiplication. Il est néanmoins assez simple de vérifier la commutativité \begin{aligned} (a,b)\cdot(c,d)&=(a\cdot c-b\cdot d,a\cdot d+b\cdot c)\nonumber\\ &=(c\cdot a-d\cdot b,d\cdot a+c\cdot b)=(c,d)\cdot (a,b).\end{aligned}

Vérifier l’associativité du produit sur notre ensemble {\mathbb{R}}^2.

Regardons à présent ce qu’il se passe si on étudie les ensemble de nombres où le deuxième nombre est nul tels que (a,0). Si on additionne deux tels nombres ont obtient (a,0)+(b,0)=(a+b,0).On constate donc que ce genre de nombre se comporte exactement comme un nombre réel normal du point de vue de l’addition. Que se passe-t-il quand on multiplie deux tels nombres (a,0)\cdot(b,0)=(a\cdot b-0\cdot 0,a\cdot 0+0\cdot b)=(a\cdot b,0). On voit que pour la multiplication également les ensembles de nombres dont le deuxième est nul, se comporte comme un nombre réel standard.

En fait on peut montrer que ce sous-ensemble de {\mathbb{R}}^2 se comporte exactement comme {\mathbb{R}}. Il se trouve que {\mathbb{R}}^2 est un ensemble de nombre plus grand que {\mathbb{R}} et qui le contient entièrement.

Les nombres complexes

Afin de simplifier les notations et les calculs, on peut introduire une notation différente. Introduisons donc le nombre imaginaire i tel que (a,b)=a+i\cdot b. On va maintenant définir l’ensemble des nombres complexes z\in{\mathbb{C}} comme tout nombre qui peut s’écrire sous la forme z=a+i\cdot b. Avec l’addition que nous avons définie à l'@eq:add, nous avons avec la nouvelle notation (a,b)+(c,d)=(a+c,b+d)\Leftrightarrow(a+i\cdot b)+(c+i\cdot d)=(a+c)+i(b+d). On constate que les nombres multipliés par i sépare nos couples de nombres (les empêche “de se mélanger”),

Pour la multiplication nous avons de même par la définition (équation @eq:mult) (a,b)\cdot(c,d)=(ac-bd,ad+bc)\Leftrightarrow(a+i\cdot b)\cdot(c+i\cdot d)=(ac-bd)+i(ad+bc).{#eq:res_mult} Si maintenant nous utilisons la multiplication de manière classique avec notre nouvelle notation (on distribue le produit comme pour les réels) (a+i\cdot b)\cdot(c+i\cdot d)=ac+i^2\cdot bd+i(ad+bc). On constate donc que pour que cette équation soit égale à l’équation @eq:res_mult on doit avoir que i^2=-1. Il se trouve que c’est la définition formelle du nombre imaginaire. Dans les réels i ne peut pas exister. En revanche dans l’espace plus grand des complexes i a une existence tout à fait naturelle et raisonnable. En fait le nombre i est associé au couple (0,1).

On appelle partie réelle d’un nombre complexe z, la partie pas multipliée par i (on la note {\mathrm{Re}}(z)) et la partie imaginaire celle multipliée par i (on la note {\mathrm{Im}}(z)). Pour z=a+ib, on a donc {\mathrm{Re}}(z)=a et {\mathrm{Im}}(z)=b.

Interprétation géométrique

Comme on l’a vu précédemment, les nombres complexes peuvent se voir comme une “notation” de {\mathbb{R}}^2. On peut donc les représenter sur un plan bidimensionnel (voir la @fig:complexPlane).

Représentation du nombre complexe z=a+ib.

La somme de deux nombres complexes s’interprête également facilement de façon graphique. On peut le voir sur la @fig:complexPlaneSum. Il s’agit en fait de simplement faire la somme des vecteurs représentant chacun des nombres complexes à sommer.

Représentation de la somme de deux nombres complexes z_1=a+ib et z_2=c+id. Le résultat est donné par z_3=a+c+i(b+d).

Pour la multiplication cela s’avère un peu plus difficile à interpréter. Pour cela il est plus simple de passer par une représentation via des sinus et des cosinus (en coordonnées cylindriques) des nombres complexes (voir la @fig:complexPlaneCyl.

Représentation du nombre complexe z=a+ib.

En utilisant la représentation en termes de \vartheta et r, on a que z=r(\cos\vartheta+i\sin\vartheta)=a+ib. On a immédiatement les relations suivantes entre ces deux représentations \begin{aligned} r=\sqrt{a^2+b^2},\\ \cos\vartheta=\frac{a}{r},\\ \sin\vartheta=\frac{b}{r}.\end{aligned} On dit que r est le module de z (aussi noté |z|) et que \vartheta est son argument (aussi noté \arg(z)).

Si à présent on définit z_1=r_1(\cos\vartheta_1+i\sin\vartheta_1) et z_2=r_2(\cos\vartheta_2+i\sin\vartheta_2), on a que z_3=z_1\cdot z_2 devient \begin{aligned} z_3=r_1r_2\left(\cos\vartheta_1\cos\vartheta_2-\sin\vartheta_1\sin\vartheta_2+i\left(\cos\vartheta_1\sin\vartheta_2+\cos\vartheta_2\sin\vartheta_1\right)\right).\end{aligned} En utilisant les relations trigonométriques suivantes \begin{aligned} \cos\vartheta_1\cos\vartheta_2-\sin\vartheta_1\sin\vartheta_2&=\cos(\theta_1+\theta_2),\\ \cos\vartheta_1\sin\vartheta_2+\cos\vartheta_2\sin\vartheta_1&=\sin(\theta_1+\theta_2),\end{aligned} il vient \begin{aligned} z_3=r_1r_2\left(\cos(\vartheta_1+\vartheta_2)+i(\sin(\vartheta_1+\vartheta_2)\right).\end{aligned} On a donc comme interprétation géométrique que le produit de deux nombres complexe donne un nombre complexe dont la longueur est le produit des longueurs des nombres complexes originaux et dont l’orientation est la somme des angles des nombres complexes originaux.

Cette propriété du produit nous amène à la notation sous forme d’exponentielle des nombres complexes. L’exponentielle, possède la propriété intéressante suivante e^a e^b=e^{a+b}. Ou encore quand on multiplie deux nombres représentés par une exponentielle, on peut représenter le résultat par l’exponentielle de la somme de leurs arguments. Comme pour les nombre complexes en somme. On a donc que z=re^{i\vartheta}=r(\cos\vartheta+i\sin\vartheta).

On peut démontrer de façon plus rigoureuse cette relation grâce aux équations différentielles. On a vu dans le chapitre précédent que l’équation différentielle f'(x)=\alpha f(x),\quad f(0)=r. a pour solution f(x)=e^{\alpha x} (\alpha\in{\mathbb{C}}). Si on remplace \alpha par i, on a f=e^{ix}. Par ailleurs, avec \alpha=i, on peut également vérifier que f(x)=r(\cos x+i\sin x) satisfait l’équation différentielle ci-dessus. On a donc bien que les deux formes sont égales.

Quelques notations et définitions

Pour la suite de ce cours, nous allons avoir besoin d’un certain nombre de notations et de définition. En particulier, nous allons noter \bar{z} le nombre complexe conjugué de z. Soit z=a+ib, son complexe conjugué {\bar{z}} est donné par {\bar{z}}=a-ib. On voit que le complexe conjugué a la même partie réelle que le nombre de départ, mais une partie imaginaire opposée.

Lors de l’utilisation de la notation polaire d’un nombre complexe, nous avons que le nombre complexe conjugué est de module égal, mais d’argument opposé. En d’autres termes, si z=re^{i\vartheta}, alors {\bar{z}}=re^{-i\vartheta}.

On peut également écrire le module d’un nombre complexe à l’aide de la notation du complexe conjugué. Il est donné par |z|=\sqrt{z{\bar{z}}}. Finalement, on peut également exprimer les parties réelle et imaginaires d’un nombre complexe à l’aide de la notation du complexe conjugué {\mathrm{Re}}(z)=\frac{1}{2}(z+{\bar{z}}),\quad {\mathrm{Im}}(z)=\frac{1}{2i}(z-{\bar{z}}).

Démontrer les trois relations précédentes.

Rajoutons encore la relation entre e^{i\theta} et les \cos,\sin. \begin{aligned} \cos(\theta)=\frac{e^{i\theta}+e^{-i\theta}}{2},\\ \sin(\theta)=\frac{e^{i\theta}-e^{-i\theta}}{2i}.\end{aligned}

Démontrer ces deux relations.

Espaces vectoriels

Ici nous introduisons de façon très simplifiée les concepts d’espaces vectoriels et certaines notions d’algèbre linéaire. Pour ce faire nous allons considérer un ensemble V, sous ensemble d’un espace plus grand E (muni d’une addition et d’une multiplication). Dans notre cas E sera {\mathbb{R}} ou {\mathbb{C}} principalement.

On appelle espace vectoriel sur E, un ensemble V, dont les éléments appelés vecteur et notés v, sont sont munis des opérations opérations + (l’addition) et \cdot (la multiplication par un scalaire) ont les propriétés suivantes

  •  

    1. L’addition est associative et commutative. Soient u,v,w\in V, alors u+v=v+u,\quad \mbox{ et }\quad (u+v)+w=u+(v+w).

    2. L’addition admet un élément neutre additif, noté 0_V, tel que 0_V+v=v.

    3. Tout v admet un opposé, noté -v tel que v+(-v)=0_V.

  •  

    1. La multiplication par un scalaire est distributive à gauche sur l’addition (et à droite sur E). Pour u,v\in V et \alpha\in E, on a \alpha\cdot(u+v)=\alpha\cdot u+\alpha\cdot v.

    2. La multiplication est associative par rapport à la multiplication de E. Soient \alpha,\beta\in E (\alpha\cdot\beta)\cdot v=\alpha\cdot(\beta\cdot v).

    3. La multiplication par un scalaire admet un élément neutre, noté 1, pour la multiplication à gauche 1 \cdot v=v.

  1. L’espace nul, v=0.

  2. L’ensemble V=E lui-même. En particulier V={\mathbb{R}} ou V={\mathbb{C}} avec l’addition et la multiplication usuelle.

  3. L’espace des n-uplets. Pour tout n>0, l’ensemble des n-uplets d’éléments de E, v=(v_1,v_2,...,v_n),\ \{v_i\}_{i=1}^n\in E, noté E^n. Sur cet espace l’addition se définit (u,v\in E^n) u+v=(u_1+v_1,u_2+v_2,...,u_v+v_n), et la mutliplication par un scalaire comme (\alpha\in E) \alpha v=(\alpha v_1,\alpha v_2,...,\alpha v_n). On a donc que l’élément neutre de l’addition est le vecteur 0_{E^n}=\underbrace{(0,0,...,0)}_{n}. L’élément opposé de v est -v=(-v_1,-v_2,...,-v_n).

    Si E={\mathbb{R}}, alors on a l’espace Euclidien. Vous avez l’habitude de l’utiliser en 2D ou 3D quand vous considérez des vecteurs. Dans ce cas {\mathbb{R}}^2 ou {\mathbb{R}}^3 avec l’addition classique et la multiplication par un scalaire standard forme un espace vectoriel.

  4. Dans ce qui suit dans ce cours, nous allons utiliser encore un autre espace vectoriel un peu moins intuitif que ceux que nous avons vus jusqu’ici. Il s’agit de l’espace des fonctions, ou espace fonctionnel. Nous définissons les applications de W dans V comme un espace vectoriel dans E avec l’addition et la multiplication par un scalaire définis commme suit. Soient f:W\rightarrow V et g:W\rightarrow V, avec \alpha\in E, alors \begin{aligned} &(f+g)(x)=f(x)+g(x), \quad \forall x\in W,\\ &f(\alpha\cdot x)=\alpha\cdot f(x), \quad \forall x\in W. \end{aligned}

  5. Espace des applications linéaires. Soit f une fonction de f:W\rightarrow V, avec W,V\in E des espaces vectoriels, alors une application est dite linéaire si \begin{aligned} &f(x+y)=f(x)+f(y),\quad \forall x,y\in W,\\ &f(\alpha \cdot x)=\alpha \cdot f(x),\quad \forall \alpha\in E,\ \mbox{et}\ x\in W. \end{aligned}

Base

Nous avons introduit la notion très générale d’espace vectorielle et nous avons présenté quelques exemples. Reprenons l’exemple de l’espace Euclidien, soit l’espace des vecteurs comme vous en avez l’habitude. Limitons nous au cas où les vecteur sont bidimensionnels, soit v=(v_1,v_2) avec v_1,v_2\in{\mathbb{R}}. D’habitude ces vecteurs sont représentés dans le système de coordonnées cartésien où on a deux vecteurs (de base) définis comme e_1=(1,0) et e_2=(0,1) qui sont implicites. Par exemple, si u=(4,5) cela signifie implicitement qu’on a u=4\cdot e_1+5\cdot e_2.

Le vecteur v dans la représentation cartésienne.

De façon générale le vecteur v=(v_1,v_2) est représenté implicitement par (voir la @fig:baseCart) v=v_1\cdot e_1+v_2\cdot e_2. On dit que e_1 et e_2 forme une base de l’espace {\mathbb{R}}^2. En d’autre terme n’importe quel vecteur v\in{\mathbb{R}}^2 peut être exprimé comme une combinaison linéaire de e_1 et e_2.

Néanmoins, le choix de la base e_1 et e_2 est totalement arbitraire. N’importe quelle autre paire de vecteurs (qui n’on pas la même direction) peut être utilisée pour représenter un vecteur quelconque dans le plan (voir la @fig:baseNonCart).

Le vecteur v dans une représentation non cartésienne.

Cette écriture en fonction de vecteurs de base, permet de faire facilement les additions de vecteurs w=u+v=u_1\cdot e_1+u_2\cdot e_2+v_1\cdot e_1+v_2\cdot e_2=(u_1+v_1)\cdot e_1+(u_2+v_2)\cdot e_2.

  1. Pour l’espace des fonctions polynomiales f(x)=\sum_{i=0}^Na_ix^i les fonction e_i=x^i forment une base.

  2. Pour l’espace vectoriel des fonctions périodiques les fonctions \sin et \cos forment une base (voir plus de détails dans ce qui suit).

Plus formellement nous allons introduire un certain nombre de concepts mathémqtiques pour définir une base. Considérons toujours V un espace vectoriel sur E.

Soient \{\alpha_i\}_{i=1}^n\in E. On dit qu’un ensemble de vecteurs \{v_i\}_{i=1}^n\in V est une famille libre si \sum_{i=1}^n \alpha_iv_i=0 \Rightarrow \alpha_i=0,\ \forall i.

  1. \{e_1\} est une famille libre de {\mathbb{R}}^2.

  2. \{e_1,e_2\} est une famille libre de {\mathbb{R}}^2.

  3. \{e_1,e_2,v\}, avec v=(1,1) n’est pas une famille libre de {\mathbb{R}}^2. En effet, 1\cdot e_1+1\cdot e_2-1\cdot v=(0,0).

  4. \{\sin(x),\cos(x)\} est une famille libre. On em peut pas écrire \sin(x)=\alpha\cos(x)+\beta. Il n’y a pas de relation linéaire qui relie les deux. La relation est non-linéaire \sin(x)=\sqrt{1-\cos^2(x)}.

On dit qu’un ensemble de vecteurs \{e_i\}_{i=1}^n\in V est une famille génératrice si \forall\ v\in V,\quad \exists \{\alpha_i\}_{i=1}^n\in E,\quad \mbox{t.q.}\quad v=\sum_{i=1}^n\alpha_i\cdot e_i. En d’autres termes, tout v\in V peut s’exprimer comme une combinaison linéaire des vecteur e_i.

  1. \{e_1\} n’est une famille génératrice de {\mathbb{R}}^2. On ne peut pas représenter tous les vecteurs de la forme v=(0,v_2), v_2\neq 0.

  2. \{e_1,e_2\} est une famille génératrice de {\mathbb{R}}^2.

  3. \{e_1,e_2,v\}, avec v=(1,1) est une famille génératrice de {\mathbb{R}}^2.

Un ensemble de vecteurs B=\{e_i\}_{i=1}^n forme une base si c’est une famille génératrice et une famille libre. En d’autres termes cela signifie qu’un vecteur v\in V peut se représenter comme une combinaison linéaire de \{e_i\}_{i=1}^n et que cette représentation est unique \forall v\in V, \quad !\exists \{\alpha_i\}_{i=1}^n\in E,\quad t.q.\quad v=\sum_{i=1}^n\alpha_i v_i. Les \alpha_i sont appelé les coordonnées de v dans la base B.

  1. \{e_1,e_2\} est une base de {\mathbb{R}}^2.

  2. \{e_1,e_2,e_3\}, avec e_3=(1,1), n’est pas une base de {\mathbb{R}}^2, car ce n’est pas une famille libre. On a par exemple que l’élément v=(0,0) peut se représenter avec les coordonnées \alpha=(0,0,0) et également les coordonnées \beta=(1,1,-1).

Introduction générale sur les séries de Fourier

Dans cette sous section, nous allons voir de façon très générale les concepts de la représentation de série de Fourier de fonctions.

Considérations historiques

Historiquement, les séries de Fourier sont apparues lorsque les mathématiciens/physiciens du 18-19ème siècles ont essayé de résoudre des équations différentielles particulières. En particulier, il y avait l’équation de la propagation d’ondes {\frac{\partial^2 \rho}{\partial t^2}}=\alpha^2\left({\frac{\partial^2 \rho}{\partial x^2}}+{\frac{\partial^2 \rho}{\partial y^2}}+{\frac{\partial^2 \rho}{\partial z^2}}\right),{#eq:ondes} où \rho est l’amplitude de l’onde et \alpha la vitesse de propagation. On a également l’équation de la chaleur {\frac{\partial T}{\partial t}}=\kappa\left({\frac{\partial^2 T}{\partial x^2}}+{\frac{\partial^2 T}{\partial y^2}}+{\frac{\partial^2 T}{\partial z^2}}\right),T est la température et \kappa la diffusivité thermique.

Ces équations ont une structure particulière. En effet, d’une part elles sont linéaires. Soient \rho_1 et \rho_2 deux solutions de l’équation @eq:ondes, on a que la somme \rho_1+\rho_2 est également solution de @eq:ondes. Cette structure d’équation différentielle impose des contraintes assez fortes sur la forme des solutions.

Par ailleurs, le fait que les dérivées à différents ordres apparaissent dans la même équation, cela impose que les fonctions et leurs dérivées à différents ordres soient reliées entre elles. Les fonctions qu’on connaît qui ont ces propriétés sont l’exponentielle et les fonctions sinus ou cosinus. Dans le cas de propagation d’ondes, on voit qu’on a uniquement des deuxièmes dérivées, et on en déduit que les fonctions importantes seront des sinus et des cosinus.

On constate que le choix du sinus ou du cosinus pour représenter ces solutions ne tombe pas du ciel. Il est dicté par les propriétés des équations que nous tentons de résoudre. En fait, nous mettons à notre disposition des outils mathématiques appropriés pour résoudre des problèmes physiques existant et qui ont des contraintes particulières.

Décomposition de signaux périodiques

Nous allons considérer une fonction f(t) qui est une fonction périodique, de période T, de pulsation \omega=2\pi/T et de fréquence \nu=1/T. Ce genre de fonction a la propriété suivante f(t+T)=f(t),\quad \forall t. Nous cherchons à décomposer f en un ensemble potentiellement infini de fonctions périodiques. Notons cet ensemble de fonctions \{g_j\}_{j=0}^\infty, où g_j est une fonction périodique. En fait on cherche une décomposition où pour un ensemble unique de \{\alpha_j\}_{j=0}^\infty f(t)=\sum_{j=0}^\infty \alpha_j g_j(t). Cette décomposition nous fait penser furieusement à une décomposition dans une base particulière, où les g_j sont les vecteurs de la base et les \alpha_j sont les coordonnées de f dans la base des g_j.

La fonction de départ f ayant une période T, on a obligatoirement que les fonctions g_j ont une période qui doit être une fraction entière de la période, T/j. Ces fonctions g_j(t) peuvent en général avoir une forme quelconque, avec l’unique contrainte qu’elles sont périodiques avec période T/j. Ça pourrait être un signal carré, triangulaire, etc. Dans les cas qui nous intéresse, on a un choix naturel qui s’impose comme fonctions périodiques: les sinus et cosinus.

Pour commencer, imaginons que nous voulions décomposer (approximer) f en une somme de g_j\sim A_j\sin(j\omega t+\phi_j). On peut jouer sur deux degrés de libertés des sinus dont la période est imposée, soit l’amplitude A_j et la phase \phi_j. On va donc écrire f(t) comme f(t)=\sum_{j=0}^\infty A_j\sin(j\omega t+\phi_j).{#eq:sin_phase_ampl} Cette forme n’est pas pratique du tout comme décomposition, en particulier à cause de la phase \phi_j. On utilise donc la relation trigonométrique (déjà utilisée pour interpréter le produit de deux nombres complexes) \sin(\theta+\phi)=\sin(\theta)\cos(\phi)+\cos(\theta)\sin(\phi). Il vient donc \begin{aligned} f(t)=\sum_{j=0}^\infty A_j\left(\sin(j\omega t)\cos(\phi_j)+\cos(j\omega t)\sin(\phi_j)\right).\end{aligned} En renommant \begin{aligned} a_j&\equiv A_j\cos(\phi_j),\\ b_j&\equiv A_j\sin(\phi_j),\end{aligned} on obtient f(t)=\sum_{j=0}^\infty \left(a_j\sin(j\omega t)+b_j\cos(j\omega t)\right). {#eq:decomp_sincos} On a donc transformé une équation où on devait déterminer une amplitude et une phase, ce qui est très compliqué, en une autre équation où on doit déterminer uniquement deux amplitude. Par ailleurs, comme \cos et \sin sont indépendants, on peut calculer les a_j et b_j de façon également indépendantes.

Nous voulons donc à présent calculer a_n et b_n pour avoir les coordonnées de f dans la base des \sin et des \cos. Pour ce faire, on va tenter de trouver les amplitudes a_j,b_j tels que les a_j\cos(j\omega t) et b_j\sin(j\omega t) approximent au mieux la fonction f.

On va donc considérer les fonctions d’erreur suivantes E^s_j=\int_0^T(f(t)-a_j\sin(j\omega t))^2{\mathrm{d}}t,\quad E^c_j=\int_0^T(f(t)-b_j\cos(j\omega t))^2{\mathrm{d}}t. Puis on va déterminer a_j,b_j tels que E_j^s et E_j^c sont minimales. Pour ce faire on va utiliser les dérivées et déterminer nos coefficients en résolvant les équations {\frac{{\mathrm{d}}E^s_j}{{\mathrm{d}}b_j}}=0,{#eq:deriv_bj} {\frac{{\mathrm{d}}E^c_j}{{\mathrm{d}}a_j}}=0.{#eq:deriv_aj} Pour l'@eq:deriv_aj, on a \begin{aligned} {\frac{{\mathrm{d}}E^c_j}{{\mathrm{d}}b_j}}&={\frac{{\mathrm{d}}\int_0^T(f(t)-b_j\cos(j\omega t))^2{\mathrm{d}}t}{{\mathrm{d}}b_j}},\nonumber\\ &=\underbrace{{\frac{{\mathrm{d}}(\int_0^Tf^2(t){\mathrm{d}}t)}{{\mathrm{d}}b_j}}}_{=0}+{\frac{{\mathrm{d}}(b_j^2\int_0^T(\cos^2(j\omega t){\mathrm{d}}t))}{{\mathrm{d}}b_j}}-{\frac{{\mathrm{d}}(2b_j\int_0^T(f(t)\cos(j\omega t){\mathrm{d}}t))}{{\mathrm{d}}b_j}},\nonumber\\ &=2b_j\int_0^T\cos^2(j\omega t){\mathrm{d}}t-2\int_0^Tf(t)\cos(j\omega t){\mathrm{d}}t,\nonumber\\ &=2b_j\frac{T}{2}-2\int_0^T\cos(j\omega t)f(t){\mathrm{d}}t.\end{aligned} Finalement on obtient b_j=\frac{2}{T}\int_0^T\cos(j\omega t)f(t){\mathrm{d}}t. Pour a_j on a de façon similaire a_j=\frac{2}{T}\int_0^T\sin(j\omega t)f(t){\mathrm{d}}t. En particulier si j=0, on a a_0=0,\quad b_0=\frac{2}{T}\int_0^T f(t){\mathrm{d}}t. On constate que b_0/2 correspond à la valeur moyenne de f(t) dans [0,T]. Cela permet d’approximer des fonctions dont la valeur moyenne n’est pas nulle (les sinus et cosinus ont toujours des moyennes nulles).

Les coefficients a_j,b_j peuvent être calculés directement à partir de f(t), comme nous venons de le voir. Nous pouvons obtenir le même résultat, en utilisant les relations suivantes (exercice) \begin{aligned} \int_0^T \sin(k \omega t)\sin(j \omega t){\mathrm{d}}t&=\delta_{jk} \frac{T}{2},\\ \int_0^T \cos(k \omega t)\cos(j \omega t){\mathrm{d}}t&=\delta_{jk} \frac{T}{2},\\ \int_0^T \sin(k \omega t)\cos(j \omega t){\mathrm{d}}t&=0,\end{aligned} qui s’obtiennent en utilisant les relations trigonométriques suivantes \begin{aligned} \sin\theta\sin\phi&= \frac{1}{2}\left(\cos(\theta-\phi)-\cos(\theta+\phi)\right),\\ \cos\theta\cos\phi&= \frac{1}{2}\left(\cos(\theta-\phi)+\cos(\theta+\phi)\right),\\ \sin\theta\cos\phi&= \frac{1}{2}\left(\sin(\theta+\phi)+\sin(\theta-\phi)\right).\end{aligned}

Cela est dû à la propriété d’othorgonalité des fonctions sinus/cosinus. En multipliant l'@eq:decomp_sincos par \frac{2}{T}\sin(k \omega t) et en intégrant entre 0 et T, on obtient \begin{aligned} \frac{2}{T}\int_0^T f(t)\sin(k\omega t){\mathrm{d}}t&=\frac{2}{T}\sum_{j=0}^\infty \left(b_j\underbrace{\int_0^T\cos(j\omega t)\sin(k\omega t){\mathrm{d}}t}_{=0}+a_j\underbrace{\int_0^T\sin(j\omega t)\sin(k \omega t){\mathrm{d}}t}_{=\frac{T}{2}\delta_{jk}}\right),\nonumber\\ \frac{2}{T}\int_0^T f(t)\sin(k\omega t){\mathrm{d}}t&=\sum_{j=0}^\infty a_j \delta_{jk}=a_k,\end{aligned}\delta_{jk} est le “delta de Kronecker”, dont la définition est \delta_{jk}=\left\{\begin{array}{ll} $1,$&$\mbox{ si }j=k$\\ $0,$&$\mbox{ sinon.}$ \end{array}\right.

En multipliant l'@eq:decomp_sincos par \frac{2}{T}\cos(k \omega t) et en intégrant entre 0 et T, on obtient \begin{aligned} \frac{2}{T}\int_0^T f(t)\cos(k\omega t){\mathrm{d}}t&=\frac{2}{T}\sum_{j=0}^\infty \left(a_j\underbrace{\int_0^T\cos(j\omega t)\sin(k\omega t){\mathrm{d}}t}_{=0}+b_j\underbrace{\int_0^T\cos(j\omega t)\cos(k \omega t){\mathrm{d}}t}_{=\frac{T}{2}\delta_{jk}}\right),\nonumber\\ \frac{2}{T}\int_0^T f(t)\cos(k\omega t){\mathrm{d}}t&=\sum_{j=0}^\infty b_j \delta_{jk}=b_k.\end{aligned}

Les séries de Fourier en notations complexes

Comme on le voit dans l'@eq:decomp_sincos, on décompose f(t) en une somme contenant des sinus et des cosinus. Cette écriture nous fait penser qu’il pourrait être possible de réécrire cette somme de façon plus concise à l’aide des nombres complexes (e^{i\theta}=\cos\theta+i\cdot\sin\theta). Effectivement cette réécriture est possible. Pour ce faire il faut définir de nouveaux coefficients c_n, c_n=\left\{\begin{array}{ll} \frac{b_n+ia_n}{2}, & $\mbox{ si }$n<0\\ \frac{b_0}{2}, & $\mbox{ si }$n=0\\ \frac{b_n-ia_n}{2}, & $\mbox{ si }$n>0 \end{array}\right. Avec cette notation, on peut réécrire l'@eq:decomp_sincos (exercice) comme f(t)=\sum_{j=-\infty}^\infty c_je^{ij\omega t}. En multipliant cette relation par \frac{1}{T}e^{-ik\omega t} et en intégrant entre -\frac{T}{2} et \frac{T}{2}, on obtient \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-ik\omega t}{\mathrm{d}}t=\frac{1}{T}\sum_{j=-\infty}^\infty c_j\int_{-\frac{T}{2}}^{\frac{T}{2}}e^{ij\omega t}e^{-ik\omega t}{\mathrm{d}}t. Pour évaluer le membre de droite de cette équation nous transformer les exponentielles en sinus/cosinus. L’intégrale du membre de droite devient \begin{aligned} \int_{-\frac{T}{2}}^{\frac{T}{2}}e^{ij\omega t}e^{-ik\omega t}{\mathrm{d}}t&=\int_{-\frac{T}{2}}^{\frac{T}{2}}\left(\cos(j\omega t)+i\sin(j\omega t)\right)\left(\cos(-k\omega t)+i\sin(-k\omega t)\right){\mathrm{d}}t,\nonumber\\ &=\int_{-\frac{T}{2}}^{\frac{T}{2}}\left(\cos(j\omega t)\cos(k\omega t)+\sin(j\omega t)\sin(k\omega t)\right.\nonumber\\ &\quad\quad\left.-i(\cos(j\omega t)\sin(k\omega t)+\cos(k\omega t)\sin(j\omega t))\right){\mathrm{d}}t,\nonumber\\ &=T\delta_{jk}.\end{aligned} En remplaçant cette relation dans l’équation ci-dessus6, on a \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-ik\omega t}{\mathrm{d}}t=\sum_{j=-\infty}^\infty c_j\delta_{jk}=c_k.{#eq:ck} Cette relation nous dit comment évaluer les coefficients c_k de la série de Fourier de f(t).

On notera que pour une fonction périodique, on obtient des coefficients de la série de Fourier qui sont discrets.

La série de Fourier pour une fonction quelconque: la transformée de Fourier

Il est possible d’écrire de telles séries pour des fonctions non-périodiques. Pour ce faire, il faut prendre la limite T\rightarrow\infty. Pour ce faire on va écrire f(t)=\sum_{j=-\infty}^\infty c_je^{ij\omega t}, où on remplace le coefficient c_j par l'@eq:ck. On obtient f(t)=\sum_{j=-\infty}^\infty \left(\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-ij\omega t}{\mathrm{d}}t\right) e^{ij\omega t}. En utilisant la relation \frac{1}{T}=\frac{\omega}{2\pi}=\frac{\omega(j-j+1)}{2\pi}=\frac{\omega(j+1)}{2\pi}-\frac{\omega j}{2\pi}, ainsi que la notation \omega_j=j\omega, on peut réécrire cette équation \begin{aligned} f(t)&=\sum_{j=-\infty}^\infty \frac{1}{2\pi}(\omega_{j+1}-\omega_j)\underbrace{\left(\int_{-\frac{\pi}{\Delta \omega_j}}^{\frac{\pi}{\Delta \omega_j}}f(t)e^{-i\omega_j t}{\mathrm{d}}t\right)}_{\equiv {\hat{f}}(\omega_j)} e^{i\omega_j t},\nonumber\\ &=\frac{1}{2\pi}\sum_{j=-\infty}^\infty (\Delta \omega_j){\hat{f}}(\omega_j) e^{i\omega_j t}.\end{aligned} Maintenant pour passer dans le cas où la fonction n’est pas périodique (la période est infinie), nous devons prendre la limite \Delta \omega_j\rightarrow 0 dans l’équation précédente, et on voit apparaître une somme de Riemann \begin{aligned} f(t)&=\frac{1}{2\pi}\sum_{j=-\infty}^\infty \lim\limits_{\Delta \omega_j\rightarrow 0}\Delta \omega_j{\hat{f}}(\omega_j) e^{i\omega_j t},\nonumber\\ &=\frac{1}{2\pi}\int_{-\infty}^\infty {\hat{f}}(\omega) e^{i\omega t}{\mathrm{d}}\omega.\end{aligned} A présent, nous avons deux opérateurs que nous allons nommer. Nous avons la transformée de Fourier {\hat{f}}(\omega)=\int_{-\infty}^{\infty}f(t)e^{-i\omega t}{\mathrm{d}}t,{#eq:fourier_transform} et la transformée de Fourier inverse f(t)=\frac{1}{2\pi}\int_{-\infty}^\infty {\hat{f}}(\omega) e^{i\omega t}{\mathrm{d}}\omega.{#eq:inverse_fourier_transform} On a immédiatement qu’appliquer la transformée de Fourier et la transformée de Fourier inverse sur une fonction f(t), nous donne la fonction originale f(t).

La fonction f(t) doit satisfaire un certain nombre de contraintes pour pouvoir calculer sa transformée de Fourier:

  1. Elle doit être de carré intégrable \int_{-\infty}^\infty |f(t)|^2{\mathrm{d}}t < \infty

  2. Elle doit avoir un nombre fini d’extrema (ne doit pas varier trop vite).

  3. Elle doit avoir un nombre fini de discontinuités.

Calculer les transformées de Fourier des fonctions suivantes

  1. Le pulse symétrique f(t)=\left\{\begin{array}{ll} 1,&\mbox{ si }-t_c<t<t_c\\ 0,&\mbox{ sinon.} \end{array}\right.

  2. Le pulse asymétrique f(t)=\left\{\begin{array}{ll} 1,&\mbox{ si } 0<t<2t_c\\ 0,&\mbox{ sinon.} \end{array}\right.

  3. L’exponentielle décroissante f(t)=\left\{\begin{array}{ll} e^{-at},&\mbox{ si } t>0\\ 0,&\mbox{ sinon.} \end{array}\right.

Calculer les transformées de Fourier inverse de la fonction suivante

  1. Le pulse symétrique f(\omega)=\left\{\begin{array}{ll} 1,&\mbox{ si }-\omega_c<\omega<\omega_c\\ 0,&\mbox{ sinon.} \end{array}\right.

Propriétés des transformées de Fourier

La transformée de Fourier possède plusieurs propriétés intéressantes.

[Propriétés]{}

  1. Linéarité. Soit une fonction h(t)=af(t)+bg(t), alors sa transformée de Fourier est donnée par {\hat{h}}(\omega)=a{\hat{f}}(\omega)+b{\hat{g}}(\omega).

  2. Translation temporelle. Soit une fonction g(t)=f(t+t_0), alors sa transformée de Fourier est donnée par {\hat{g}}(\omega)={\hat{f}}(\omega)e^{i\omega t_0}.

  3. Modulation en fréquence. Soit \omega_0\in{\mathbb{R}} et une fonction g(t)=e^{-i\omega_0 t}f(t), alors sa transformée de Fourier est donnée par {\hat{g}}(\omega)={\hat{f}}(\omega+\omega_0).

  4. Contraction temporelle. Soit a\in{\mathbb{R}}^\ast et g(t)=f(at) alors sa transformée de Fourier est donnée par {\hat{g}}(\omega)=\frac{1}{|a|}{\hat{f}}(\omega/a). En particulier, on a la propriété d’inversion du temps quand a=-1, on a h(t)=f(-t)\Rightarrow{\hat{h}}(\omega)={\hat{f}}(-\omega).

  5. Spectres de fonctions paires/impaires. Soit f(t) une fonction paire (impaire), alors {\hat{f}}(\omega) sera une fonction paire (impaire).

La transformée de Fourier à temps discret (TFTD)

Nous allons maintenant plus considérer une fonction continue, mais une série de valeurs discrètes. Notons f[n] une série de nombres, avec n\in{\mathbb{Z}}. Nous voulons définir l’équivalent de la transformée de Fourier de l'@eq:fourier_transform pour ce genre de séries de points. Une façon naturelle de définir l’équivalent à temps discret de cette équation est {\hat{f}}(\omega)=\sum_{n=-\infty}^\infty f[n] e^{-i\omega n}.{#eq:tftd} Pour les transformées de Fourier à temps continu et non périodique, nous avons que la transformée de Fourier est continue et en général non périodique. Pour le cas de la transformée de Fourier à temps discret la transformée de Fourier sera périodique, soit {\hat{f}}(\omega+2\pi)={\hat{f}}(\omega). Nous démontrons cette relation par la définition de la TFTD {\hat{f}}(\omega+2\pi)=\sum_{n=-\infty}^\infty f[n] e^{-i(\omega+2\pi) n}=\underbrace{e^{-i2\pi}}_{=1}\sum_{n=-\infty}^\infty f[n] e^{-i\omega n}={\hat{f}}(\omega). D’une certaine façon nous voyons que nous avons une similarité entre la transformée de Fourier à temps discret et les séries de Fourier. Cette similarité va devenir plus claire dans ce qui suit.

Pour définir la transformée de Fourier en temps discret inverse, nous nous inspirons de la version en temps continu (voir l’équation @eq:inverse_fourier_transform) et on a f[n]=\frac{1}{2\pi}\int_{-\pi}^\pi{\hat{f}}(\omega)e^{i\omega n}{\mathrm{d}}\omega. {#eq:tftdi} Pour prouver cette relation, il suffit de remplacer l’équation @eq:tftd dans cette relation, et il vient f[n]=\frac{1}{2\pi}\int_{-\pi}^\pi \left(\sum_{m=-\infty}^\infty f[m] e^{-i\omega m}\right) e^{i\omega n}{\mathrm{d}}\omega. En supposant que la somme converge, nous pouvons intervertir la somme et l’intégrale et on a \begin{aligned} f[n]&=\frac{1}{2\pi}\left(\sum_{m=-\infty}^\infty f[m] \int_{-\pi}^\pi e^{-i\omega (m-n)} {\mathrm{d}}\omega\right),\nonumber\\ &=\frac{1}{2\pi}\left(\sum_{m=-\infty}^\infty f[m] \delta_{mn} 2\pi\right),\nonumber\\ &=f[n].\nonumber\end{aligned}

Calculer les transformées de Fourier (inverses quand c’est approprié) en temps discret des fonctions suivantes

  1. Le pulse symétrique {\hat{f}}(\omega)=\left\{\begin{array}{ll} 1,&\mbox{ si }-\omega_c<\omega<\omega_c\\ 0,&\mbox{ sinon.} \end{array}\right.

  2. Le pulse discret f[n]=\left\{\begin{array}{ll} 1,&\mbox{ si }n=0\\ 0,&\mbox{ sinon.} \end{array}\right.

Il est intéressant de noter qu’on peut représenter une suite discrète et infinie de points par une fonction continue et périodique.

La transformée de Fourier discrète

Motivation

Pourquoi avons-nous besoin d’encore une transformée de Fourier? Nous avons déjà vu la transformée de Fourier de fonctions périodiques, de fonctions non-périodiques, ainsi que de fonctions à temps discret. Néanmoins, même dans le cas de la transformée de Fourier à temps discret, la transformée de Fourier est une fonction continue. Cela n’est évidemment pas pratique ni même utilisable dans un ordinateur. C’est pourquoi il est nécessaire de définir une transformée de Fourier discrète qui aura les propriétés suivantes

  1. Elle transformera un signal discret de longueur finie.

  2. La transformée de Fourier sera discrète et de longueur finie.

Applications

Avant de voir en détail comment on calcule la transformée de Fourier discrète, on peut discuter quelle est son application. La TFD est utilisée tout le temps en traitement du signal. En gros c’est une approximation de la transformée de Fourier à temps discret. A chaque fois qu’on désire connaître le comportement d’une fonction dans l’espace spectral, on utilisera la TFD. Un exemple typique est l’application pour téléphones portables Shazam que vous connaissez sans doute. Le but de cette application est l’identification de chansons. Elle fonctionne de la façon suivante. Dans un premier temps elle enregistre un signal sonore. Puis avec ce signal sonore elle crée un spectrogramme (une sorte d’emprunte digitale de la chanson) qui est obtenu à l’aide de TFD. Finalement le spectrogramme est comparé avec une base de donnée de spectrogrammes et la chanson peut ainsi être identifiée. Une autre application est le filtrage de signaux. Comme vous l’avez vu (ou verrez) dans les travaux pratiques, la TFD rend très simple le filtrage de fréquences (ou de bande de fréquences). En effet, il suffit d’ôter de la TFD d’un signal les amplitudes voulues et d’effectuer la transformée de Fourier discrète inverse (TFDI) du signal filtré. Ce genre d’applications est très utilisé dans le domaine de la compression de données (jpg, mp3, ...).

La transformée de Fourier discrète à proprement parler

Soit f[n] un séquence de points N points, n=0..N-1. Pour se ramener au cas de la transformée de Fourier à temps discret, on peut aussi se dire qu’on a une séquence infinie de points, mais où f[n]=0, pour n\geq N. On dit qu’on a N échantillons de f.

Avec cette définition il est simple de calculer la transformée de Fourier à temps discret {\hat{f}}(\omega)=\sum_{n=0}^{N-1} f[n] e^{-i\omega n}.{#eq:tftd_fini} On note que la somme à présent ne se fait plus dans l’intervalle (-\infty,\infty), mais uniquement entre [0,N-1], car le signal est de longueur finie.

On représente donc un signal de longueur finie f[n] (n=0,..,N-1) par une fonction continue de la pulsation, {\hat{f}}(\omega). Les deux représentations sont équivalentes. On en déduit que l’information contenue dans un nombre fini de points, est la même que dans une fonction continue (et donc contenant une infinité de points). Une partie de l’information contenue dans la fonction continue doit être redondante...

L’idée à présent va être d’enlever toute l’information redondante de {\hat{f}}(\omega) en échantillonnant {\hat{f}} et en gardant uniquement N échantillons de {\hat{f}}. La fréquence d’échantillonage sera de 2\pi/N et le domaine d’échantillonage sera [-\pi,\pi).

Nous pouvons à présent définir mathématiquement cet échantillonage de {\hat{f}}(\omega) comme étant une suite de points, notée \{{\hat{f}}(\omega_k)\}_{k=0}^{N-1}, où \omega_k=k/(2\pi). Cette suite sera notée {\hat{f}}[k] et appelée la transformée de Fourier discrète de f[n].

On a donc que la transformée de Fourier discrète de f[n] est donnée par {\hat{f}}[k]=\sum_{n=0}^{N-1} f[n] e^{-i\omega_k n} =\sum_{n=0}^{N-1} f[n] e^{-\frac{2\pi i n k}{N}}.{#eq:tfd} En s’inspirant de définition de la transformée de Fourier inverse à temps discret de {\hat{f}}(\omega) (voir l’équation @eq:tftdi), on a que la transformée de Fourier discrète inverse est donnée par f[n]=\frac{1}{N}\sum_{k=0}^{N-1} {\hat{f}}[k] e^{i\omega_k n} =\frac{1}{N}\sum_{k=0}^{N-1} {\hat{f}}[k] e^{\frac{2\pi i k n}{N}}. Montrons à présent que la transformée inverse discrète de la transformée de Fourier discrète donne bien la suite de départ \begin{aligned} f[n]&=\frac{1}{N}\sum_{k=0}^{N-1} {\hat{f}}[k] e^{\frac{2\pi i k n}{N}},\nonumber\\ &=\frac{1}{N}\sum_{k=0}^{N-1} \sum_{m=0}^{N-1} f[m] e^{-\frac{2\pi i k m}{N}} e^{\frac{2\pi i k n}{N}},\nonumber\\ &=\frac{1}{N}\sum_{k=0}^{N-1} \sum_{m=0}^{N-1} f[m] e^{\frac{2\pi i k (n-m)}{N}},\nonumber\\ &=\frac{1}{N}\sum_{m=0}^{N-1} f[m] \sum_{k=0}^{N-1} e^{\frac{2\pi i k (n-m)}{N}},\nonumber\\ &=\frac{1}{N}\sum_{m=0}^{N-1} f[m] N \delta_{nm},\nonumber\\ &=f[n].\end{aligned} Cette relation montre qu’on a bien la même information dans la suite de longueur finie {\hat{f}}[k] que dans f[n]. On a donc enlevé avec succès toute information redondante contenue dans {\hat{f}}(\omega).

On peut maintenant de façon simple implanter la transformée de Fourier discrète sur un ordinateur car on a discrétisé toutes les étapes du calcul. Néanmoins les formules ci-dessus ne sont pas d’une grande efficacité. En effet, on peut montrer que la complexité de l’équation @eq:tfd est de l’ordre N^2.

On peut écrire l'@eq:tfd comme un produit matrice-vecteur sous la forme suivante \begin{array}{l} \underbrace{\begin{pmatrix} {\hat{f}}[0] \\ {\hat{f}}[1] \\ f[2] \\ \vdots \\ {\hat{f}}[N-1] \end{pmatrix}}_{\hat{\bm{f}}} = \underbrace{\begin{pmatrix} 1 & 1 & 1 & \cdots & 1\\ 1 & w & w^2 & \cdots & w^{N-1}\\ 1 & w^2 & w^4 & \cdots & w^{2(N-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots &\\ 1 & w^{N-1} & w^{2(N-1)} & \cdots & w^{(N-1)^2} \end{pmatrix}}_{\bm{W}}\cdot \end{array} \underbrace{\begin{pmatrix} f[0] \\ f[1] \\ f[2] \\ \vdots \\ f[N-1] \end{pmatrix}}_{\bm{f}},w = e^{-\frac{2 \pi i}{N}}. On peut donc de façon plus compacte l’écrire \hat{\bm{f}}=\bm{W}\cdot \bm{f}. Les éléments de la matrice \bm{W} peuvent être précalculés et il reste donc à calculer uniquement le produit matrice vecteur \bm{W}\cdot\bm{f}. Pour ce faire il faut pour chaque ligne de \hat{\bm{f}} induit le calcul de N produit et N sommes (donc une complexité N). Comme il y a N lignes à \hat{\bm{f}}, il y a donc N\cdot N de complexité.

Il existe des algorithmes beaucoup plus efficaces pour effectuer de genre de calculs que nous allons brièvement discuter maintenant. Ils réduisent la complexité algorithmique à N\log(N) en général. Nous allons brièvement discuter un de ces algorithmes dans la sous-section @sec:tfr.

La transformée de Fourier discrète étant un échantillonage de la transformée de Fourier à temps discret, toutes les propriétés discutées pour la transformée de Fourier à temps discret restent valides. En particulier la transformée de Fourier discrète est périodique, de période N {\hat{f}}[k]={\hat{f}}[k+N]. A démontrer en exercice.

La transformée de Fourier rapide {#sec:tfr}

L’algorithme présenté ici est une version “simplifiée” de l’algorithme de Cooley-Tukey (publié en 1965). Cet algorithme a en fait été “inventé” par Gauss en 1805 quand il essayait d’interpoler la trajectoires d’astéroides dans le système solaire.

L’idée de l’algorithme radix-2 est d’abord de séparer le signal en deux parties. D’une part les indices pairs et d’autres part les indices impairs \begin{aligned} &\left\{f[2m]\right\}_{m=0}^{N/2-1}=\left\{f[0],f[2],...,f[N-2]\right\},\\ &\left\{f[2m+1]\right\}_{m=0}^{N/2-1}=\left\{f[1],f[3],...,f[N-1]\right\}.\end{aligned} Puis les transformées de Fourier discrètes de chacune de ces sous-suites sont calculées et combinées pour avoir la transformée de Fourier du signal en entier. En fait on va appliquer cette décomposition de façon récursive sur chacune des deux parties. On fait donc l’hypothèse que la longueur du signal est une puissance de 2. Ce n’est en pratique pas un problème, car on peut facilement rajouter des “zéros” dans notre signal pour avoir un signal d’une longueur d’une puissance de 2.

Commençons donc par réécrire la transformée de Fourier {\hat{f}}[k] lorsqu’on a décomposé le signal en deux sous-signaux \begin{aligned} f[k]&=\sum_{m=0}^{N/2-1} f[2m]e^{-\frac{2\pi i (2m) k}{N}}+\sum_{m=0}^{N/2-1} f[2m+1]e^{-\frac{2\pi i (2m+1) k}{N}},\nonumber\\ &=\sum_{m=0}^{N/2-1} f[2m]e^{-\frac{2\pi i m k}{N/2}}+e^{-\frac{2\pi i k}{N}}\sum_{m=0}^{N/2-1} f[2m+1]e^{-\frac{2\pi i m k}{N/2}},\nonumber\\ &=\hat{p}[k]+e^{-\frac{2\pi i k}{N}}\hat{j}[k],\end{aligned} où nous avons défini les transformées de Fourier discrètes des parties paires et impaires p[k] et \hat{j}[k] \begin{aligned} \hat{p}[k]&=\sum_{m=0}^{N/2-1} f[2m]e^{-\frac{2\pi i m k}{N/2}},\\ \hat{j}[k]&=\sum_{m=0}^{N/2-1} f[2m+1]e^{-\frac{2\pi i m k}{N/2}}.\end{aligned} La transformée de Fourier discrète étant périodique (comme l’est la transformée de Fourier à temps discret), nous avons les propriétés suivantes \begin{aligned} \hat{p}[k]&=\hat{p}[k+N/2],\\ \hat{j}[k]&=\hat{j}[k+N/2].\end{aligned} De plus, nous avons que e^{-\frac{2\pi i (k+N/2)}{N}}=e^{-\pi i}e^{-\frac{2\pi i k}{N}}=-e^{-\frac{2\pi i k}{N}}. Avec ces propriétés il est aisé de réécrire {\hat{f}}[k]=\left\{\begin{array}{ll} \hat{p}[k]+e^{-\frac{2\pi i k}{N}} \hat{j}[k],&\mbox{ si }0\leq k<N/2\\ \hat{p}[k]-e^{-\frac{2\pi i k}{N}} \hat{j}[k],&\mbox{ si }N/2\leq k<N \end{array}\right. On a donc réduit le nombre de calculs nécessaires pour calculer {\hat{f}}[k] d’un facteur 2. En continuant cette procédure jusqu’à N=2 on peut montrer qu’on réduit la complexité algorithmique à N\log N (mais on ne le démontrera pas dans ce cours).

Fréquence d’échantillonage

Une question primordiale dans le calcul des transformée de Fourier (ou de l’analyse spectrale plus généralement) est la question de l’échantillonage du signal que nous souhaitons analyser. Dans le monde réel un signal sonore, une image,... est considéré comme une quantité continue (il est représentée par une infinité de valeur). Lorsque nous souhaitons faire une analyse spectrale sur un ordinateur de ce signal, il est nécessaire de le digitaliser: de le rendre discret. Dès lors une question très importante est de savoir quelle est la fréquence à laquelle on va enregistrer les valeurs de notre suite temporelle afin de garder toute l’information contenue dans le signal original.

En termes mathématiques, nous avons un signal f(t) que nous enregistrons entre t_0 et t_{N-1}. Nous voulons le transformer en un signal de longueur N finie, f(t_n) avec 0\leq n \leq N-1 afin de pouvoir le représenter sur un support numérique. Pour simplifier on va supposer que l’enregistrement se fait à intervalle régulier, \delta t=\frac{t_{N-1}-t_0}{N-1}. On a donc que t_n=t_0+\delta t n. La question qu’on se pose est quelle doit être la valeur de N pour ne pas perdre d’information sur f(t) quand on échantillonne. En d’autres termes à partir de quel nombre N d’échantillons la transformée de Fourier discrète de f[n] ne change plus.

Le théorème de Shannon-Nyquist nous dit que pour pouvoir représenter exactement un signal avec une fréquence maximale F_c=1/\delta t_c, alors on doit l’échantillonner avec une fréquence 1/\delta t_e=F_e\geq 2F_c. De façon similaire, si on choisit un signal et qu’on peut l’échantillonner avec une certaine précision (on détermine la fréquence maximale, F_c qu’on veut pouvoir représenter dans le signal) on a simplement besoin de choisir une fréquence d’échantillonage F_e\geq 2F_c. Nous notons F_N=2F_c la fréquence de Nyquist. En prenant F_e=F_N on a que N=1/F_e=1/F_N et que l’échantillonage permet de représenter les fréquences plus petites que F_N/2. Si la fréquence d’échantillonage est plus petite que la fréquence de Nyquist de notre signal, on verra apparaître le phénomène de repliement de spectre (aliasing en anglais).

Probabilités et statistiques

Introduction à la statistique descriptive

En statistique, une population est un ensemble d’objets (d’individus) possédant un ou plusieurs caractères communs. L’étude des caractères d’une population a pour but de révéler des tendances au sein de la population. Ces études sont particulièrement intéressantes quand le nombre d’individus de notre population est trop élevé pour pouvoir être analysé en entier. On prélève alors un échantillon représentatif de notre population au hasard et on mène l’analyse statistique sur ce sous ensemble. Les éventuelles conclusions tirées de l’étude statistique sur le sous ensemble seront ensuite appliquée à l’ensemble de la population. Grâce au calcul de probabilité nous pourrons alors avoir une confiance plus ou moins grande dans les conclusions tirées en fonction de la taille de l’échantillon. En effet plus celui-ci sera grand, plus la confiance dans les résultats sera élevée.

Un exemple de ce genre d’étude qui est très à la mode ces temps est le sondage (concernant le résultat d’élections ou de votations). Les sondeurs tentent en questionnant un sous-ensemble d’environ 1000 d’électeurs d’un pays (citoyens de plus de 18, moitié d’hommes et de femmes plus ou moins, ...) de déterminer les résultats d’élections ou de votations où participeront des millions d’électeurs potentiels. Il faut avouer que la tâche semble pour le moins complexe. Et la plus grande difficulté tient dans le “représentatif de la population”.

Représentations

Il existe différentes façon de représenter les caractères d’une population selon que sa nature est discrète ou continue. Dans le cas discret d’un caractère pouvant prendre k\in{\mathbb{N}} valeur différentes \{x_i\}_{i=0}^{k-1}, on représente le nombre d’individus pouvant prendre la valeur x_i par le nombre n_i. On a donc un ensemble \{n_i\}_{i=0}^{k-1} d’individus pour les k valeurs des caractères de la population. Dans le cas continu le nombre d’individus d’un caractère correspondrait à une subdivision en k parties de l’ensemble des valeurs possibles pour le dit caractère.

  1. Cas discret: On étudie la distribution de salaires annuels dans une entreprise. Les salaires possibles sont 40'000, 50'000, 60'000 et 1'000'000 de CHF.

    • Il y a 35 personnes payées 40'000 CHF.

    • Il y a 20 personnes payées 50'000 CHF.

    • Il y a 5 personnes payées 60'000 CHF.

    • Il y a 1 personne payée 1'000'000 CHF.

  2. Cas continu: Lors du benchmark d’une application, A, nous effectuons plusieurs mesures (la population) du temps d’exécution (le caractère) de l’application. Les résultats obtenus sont les suivants:

    • 7 exécutions ont pris entre 50 et 51 secondes.

    • 12 exécutions ont pris entre 51 et 52 secondes.

    • 8 exécutions ont pris entre 52 et 53 secondes.

    • 23 exécutions ont pris entre 53 et 54 secondes.

Pour représenter de façon un peu plus parlante ces valeurs, deux méthodes principales existent: le tableau ou le graphique. Pour illustrer les exemples précédents sous forme de tableau on obtient pour le cas des salaires (voir Tabl. @fig:salaires)

Salaire Nombre de salariés


40000            35
50000            20
60000            5
1000000 1

Tableau du nombre de salariés par salaire. {#tbl:salaires}

et du benchmark de l’application (voir Tabl. @fig:exec)

Temps d’exécution     Nombre

    \[50,51)             7
    \[51,52)            12
    \[52,53)             8
    \[53,54)            23

: Tableau du temps d'exécution et du nombre d'exécutions. {#tbl:exec}

Sous forme de graphique on peut représenter le tableau des salaires sous la forme d’un graphique bâton (voir Fig. @fig:salaires)

Nombre salariés en fonction du salaire.

ou d’un histogramme pour le temps d’exécution de l’application (voir Fig. @fig:exec).

Nombre d’exécutions en fonction du temps d’exécution.

Fréquences

Plutôt que de faire apparaître le nombre d’individus d’une population possédant un caractère, il peut être plus intéressant et parlant de faire intervenir la fréquence ou le nombre relatif à la place. En effet, la fréquence donne immédiatement la proportion d’individu plutôt qu’un nombre absolu qui n’est pas forcément très interprétable tout seul.

La population totale, n, est donnée par n=\sum_{i=0}^{k-1}n_i. On peut donc définir la fréquence d’un caractère i, f_i comme f_i=\frac{n_i}{n}.

[Fréquence]{}

Les tableaux de fréquence des deux exemples précédents sont donnés par

  1. Cas discret: la population totale est de n=35+20+5+1=61.

    Salaire Nombre de salariés Fréquence


    40000            35           $35/61\cong0.573770$
    50000            20           $20/61\cong0.327869$
    60000            5            $5/61\cong0.081967$
    1000000 1 1/61\cong0.016393

    Tableau des salaires, du nombre de salariés et la fréquence.

  2. Cas continu: la population totale est de n=7+12+8+23=50. Le tableau @tbl:exec_freq affiche les différentes fréquences des temps d’exécution.

    Temps d’exécution Nombre Fréquence


       \[50,51)          7      $7/50=0.14$
       \[51,52)         12      $12/50=0.24$
       \[52,53)          8      $8/50=0.16$
       \[53,54)         23      $23/50=0.46$

    : Tableau des temps d'exécution et la fréquence des temps d'exécution. {#tbl:exec_freq}

La fréquence possède un certain nombre de propriétés que nous retrouverons dans les sections suivantes qui sont assez intuitives

[Propriétés de la fréquence]{}

  1. Les fréquences sont toujours dans l’intervalle [0,1] 0\leq f_i\leq 1.

  2. La somme de toutes les fréquences donne toujours 1 \sum_{i=0}^{k-1} f_i = 1.

Relié avec la propriété 2 ci-dessus, il peut également être intéressant d’obtenir la fréquence cumulée, notée F(x), d’un caractère qui se définit comme la fréquence des individus qui présentent une valeur de caractère x_i\leq x. Les tableaux correspondants aux tableaux @tbl:salaires et @tbl:exec (voir le @tbl:salaires_freqcum et le @tbl:exec_freqcum)

Salaire Nombre de salariés Fréquence Fréquence cumulée