Table des matières

Ipython notebook : cours INPROS ¶

Auteur: Marc BUFFAT, Pr dpt de Mécanique, UCB Lyon 1

Contributeurs: Violaine Louvet, Michel Kern, Loic Gouarin, Laurence Viry </h5>

Mise à disposition selon les termes de la Licence Creative Commons
Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France.

Algorithmes numériques de base¶

Contenu de la leçon¶

1. Algorithme recursif
2. Tri d’un tableau
3. Recherche de racine
4. Polynômes
5. Interpolation
6. Intégration
7. Équations différentielles
8. Exercices

Bibliothèques scientifiques sous Python¶

• numpy [http://www.numpy.org/]
• scipy [http://docs.scipy.org/doc/scipy/reference]

Algorithme récursif¶

récursion méthode algorithmique où la solution du problème s’exprime en fonction de solutions du même problème, mais sur des données plus petites.

• Très utilisé en informatique, car fournit des solutions élégantes.

• Attention, peut être inefficace (comparé à l’itération)

Exemple 1: calcul du PGCD de 2 entiers¶

• PGCD(a,b)=PGCD(a-b,b) si a>b
• PGCD(a,b)=PGCD(b-a,a) si a<b
• PGCD(a,a)=a

Algorithme PGCD recursif¶

Algorithme PGCD(a,b)
  si a == b 
     retour a
  sinon si a > b
     retour PGCD(a-b,b)
  sinon
     retour PGCD(b-a,a)

Programme Python¶

Exemple 2: calcul du déterminant d’une matrice¶

soit $A$ une matrice d’ordre $n$, la méthode de Cramer permet le calcul récursif du déterminant par développement par rapport à la 1ere colonne

$$ det(A) = \sum_{k=1}^n (-1)^{k+1} A_{k,1} \times det(A^k) $$

où $A^k$ est la sous-matrice obtenu par élimination de la colonne 1 et de la ligne $k$ de $A$

Pour une matrice de dimension 1 $$ det(a) = a $$

Algorithme déterminant¶

 Algorithme determinant(A)
     n = dim(A)
     si n==1 alors
         retour A[0,0]
     fin si
     # A1 sous matrice d'ordre n-1
     A1 = tableau(n-1,n-1)
     det = 0
     signe = 1
     # développement par rapport a la 1ere colonne
     pour k de 0 a n-1
         # sous matrice Ak
         A1[0:k,:] = A[0:k ,1:n]
         A1[k:,:]  = A[k+1:,1:n]
         det = det + signe*A[k,0]*determinant(A1)
         signe = -signe
         fin pour
     retour det

Programme Python¶

Exercise: tour de Hanoi¶

Ecrire un algorithme récursif pour résoudre le problème des tours de Hanoi

Tri d’un tableau ¶

tri par insertion¶

algorithme de tri classique: on compare successivement chaque élément du tableau en l’insérant au non endroit dans le tableau (méthode classique de tri d’un jeu de carte)

principe¶

algorithme¶

# Tri d'un tableau X de n éléments
Algorithme Tri_insertion(X)
    n=dimension de X
    pour i de 1 a n-1 faire  # trie des élts
        x = T[i]  # élément a trier
        j = i     # position
        tant que (j>0) et (T[j-1]>x) faire
            T[j] = T[j-1] # permute j et j-1
            j = j -1  # nouvelle position de x
        fin tant que
        T[j] = x  # position de l’élément 
    fin i
    retour # le tableau X est trie

Visualisation de l’algorithme en Python

Efficacité $O(n^2)$

• algorithme assez efficace pour $n$ petit ou des tableaux presque triés
• tri sur place

Algorithme de tri de la bibliothéque numpy¶

sort: méthode de trie plus efficace (tri rapide ou quicksort), renvoi le tableau trié

Recherche de racine $F(x)=0$ ¶

• méthode itérative à partir d’une estimation $x_0$
• methode de Newton: $x_{n+1}= x_n - \frac{F(x_n)}{F’(x_n)}$
• méthode de la sécante: Newton avec $F’(x_n) \approx \frac{F(x_n)-F(x_{n-1})}{x_n - x_{n-1}}$

fonction fsolve dans scipy.optimize

• généralisation dans $\mathbf{R}^p$, on remplace $F’(x_n)$ par la matrice Jacobienne $J = \{\frac{\partial F_i}{\partial x_j}\} $

Exercice¶

Un projectile de masse $M$ est tiré avec un angle $\alpha$ et une vitesse initiale $U_0$. En tenant compte de la résistance de l’air qui exerce une force proportionnelle à la vitesse $F=K U$, la trajectoire du projectile est donnée par: \begin{eqnarray} x(t)&=& C\cos{\alpha} (1-e^{-t/C})\ y(t)&=& (C\sin{\alpha}+32C^2)(1-e^{-t/C}) - 32 C t \end{eqnarray} avec $C=M/K$. Écrire un programme python permettant de calculer le temps et la distance du point d’impact.

AN: $\alpha=45°$, $U_0=40 m/s$

Polynômes¶

soit le polynôme $P(x) = \sum_{i=0}^n a_i x^i $

Pour évaluer ce polynôme, on peut programmer direction la relation précédente:

Pour évaluer $P(x)$ on peut aussi utiliser la forme factorisée de Hörner, en récrivant $P(x)$ $$ P(x) = a_0 + x* (a_1 + x*(a_2 + x*( \cdots x*(a_{n_1} + a_n =*x)$$

La formule factorisée de Hörner est donc beaucoup plus efficace!! C’est celle utilisée dans la bibliothéque numpy

Interpolation¶

Interpolation de Lagrange¶

Soient $N+1$ points $\left\{ X_i,Y_i \right\}_{i=0..N}$, le polynôme d’interpolation de Lagrange est le polynôme de degré $N$ passant par ces $N+1$ points.

Formule de Lagrange (instable si $N>10$)¶

$$P(x) = \sum_{i=0}^N Y_i L_i(x) \mbox{ avec } L_i(x) = \frac{\prod_{j=0, j\neq i}^N (X-X_j)}{\prod_{j=0, j\neq i}^N (X_i-X_j)} $$

    # bibliotheque scipy
    scipy.interpolate.lagrange

Pour $N$ grand: interpolation locale (linéaire, spline, ..)¶

    # bibliotheque scipy
    scipy.interpolate.interp1d

Droite des moindres carrées¶

Soient $n$ couples de points expérimentaux: $(X_i,Y_i)_{i=0,n-1}$, on cherche à lisser ces points par la droite des moindres carrés: $y = a x +b $

**Principe**

on détermine la droite $y = a x + b$ qui minimise l’erreur $Err$ entre les points de mesure et les points de la droite: $$ Err(a,b) = \sum_{i=0}^{n-1} \Delta_i ^2 $$ où $\Delta_i ^2 = (a X_i +b - Y_i)^2$ représente le carré de l’erreur entre un point de mesure $(X_i,Y_i)$ et le point correspondant de la droite $(X_i, a X_i + b)$

Les valeurs de $a$ et $b$ sont déterminées en minimisant l’erreur $Err(a,b)$ par rapport à $a$ et $b$, i.e. tel que $\frac{\partial Err}{\partial a} = 0 $ et $\frac{\partial Err}{\partial b} = 0$

On obtient ainsi le système linéaire de 2 équations: $$ a \overline{X} + b = \overline{Y} \mbox{ avec } \overline{X}=\frac{1}{n}\sum_{i=0}^{n-1}{X_i} \mbox{ et } \overline{Y}=\frac{1}{n}\sum_{i=0}^{n-1}{Y_i}$$ $$ a \overline{X^2} + b \overline{X} = \overline{XY} \mbox{ avec } \overline{X^2}=\frac{1}{n}\sum_{i=0}^{n-1}{X_i^2} \mbox{ et } \overline{XY}=\frac{1}{n}\sum_{i=0}^{n-1}{X_i Y_i}$$ dont la solution s’écrit:

$$ a = \frac{\overline{XY}- \overline{X}\ \overline{Y}}{\overline{X^2}- (\overline{X})^2} \mbox{ et } b = \overline{Y} - a \overline{X} $$

Intégration¶

formule de quadrature par approximation de la fonction $f(x)$ par un polynôme d’interpolation $P(x)$ $$ \int_a^b{f(x)\,dx} \approx \int_a^b{P(x)\,dx} = \sum_{i=0}^N w_i f(x_i) $$

méthode des trapèzes (P(x) polynôme de degré 1 par morceaux)¶

$$ \int_a^b{f(x)\,dx} \approx \sum_{i=1}^N h\frac{f(a+ih)+f(a+(i-1)h)}{2} \mbox{ avec } h=\frac{b-a}{N}$$

méthode d’ordre 2: erreur $O(h^2)$

# librairie scipy
scipy.integrate.cumtrapz

méthode de Simpson (P(x) polynôme de degré 2 par morceaux)¶

$$ \int_a^b{f(x)\,dx} \approx \sum_{i=1}^N h\frac{f(a+ih)+4f(a+(i-\frac{1}{2})h)+f(a+(i-1)h)}{6} \mbox{ avec } h=\frac{b-a}{N}$$

méthode d’ordre 4: erreur $O(h^4)$

    # librairie scipy
    scipy.integrate.simps

méthode de Gauss¶

• on utilise des points de Gauss comme points d’interpolation
• avec $N$ points de Gauss, formule exacte pour des polynômes de degré $\leq 2n+1$

• très précise pour $f(x)$ analytique

    # bibliotheque scipy
    scipy.integrate.quadrature

Comparaison des méthodes¶

Calcul de $I = \int_0^\pi \sin x\, dx = 2$

Equations différentielles¶

Equation différentielle du 1er ordre¶

$$ \frac{d Y}{d t} = F(t,Y) \mbox{ avec la C.I.} Y(0)=Y_0$$

• $Y(t)$ est une courbe dont l’équation donne la tangente
• calcul de la trajectoire à partir de la vitesse
• approximation de $\frac{d Y}{d t}$ par différences finies $$ \frac{d Y}{d t} = \lim_{\Delta t \rightarrow \infty}{\frac{\Delta Y}{\Delta t}} \approx \frac{Y_{n+1}-Y_n}{\Delta t}$$

Méthode de Runge Kutta d’ordre 2¶

$$ Y_{n+1} = Y_n + \Delta t F(Y_{n+1/2},t_n + \frac{\Delta t}{2}) \mbox{ avec } Y_{n+1/2} = Y_n + \frac{\Delta t}{2} F(Y_n,t_n) $$

• approximation de la tangente au milieu de l’intervalle
• méthode d’ordre 2

# Bibliothéque scipy.integrate.ode 
methode = scipy.integrate.ode(F).set_integrator(nom)
# nom = methode d'integration
# 'dopri5' methode Runge Kutta d'ordre 4(5)
# 'lsoda'  methode implicite (Adams Bashford et BDF)
methode.set_initial_value(Y0,t0)
methode.integrade(t)

exemple: intégration de dY/dt = Y avec Y(0)=1¶

solution $Y(t)=e^t$

Equation différentielle du second ordre¶

$$\frac{d^2 y}{d t^2} = f(t,y,\frac{dy}{dt}) \mbox{ avec les C.I.: } y(0)=y_0, \frac{dy}{dt}(0)=u_0$$

• transformation en un système EDO d’ordre 1: $$ \frac{d Y}{dt} = F(t,Y) $$
• on pose $Y_1 = y$ et $Y_2 = \frac{dy}{dt}$ $$ \frac{d}{dt} \left[\begin{array}{c} Y_1 \ Y_2 \ \end{array}\right] = \left[\begin{array}{c} Y_2 \ f(t,Y_1,Y_2) \ \end{array}\right] \mbox{ avec } Y_1(0)=y_0 \mbox{ et } Y_2(0)=u_0$$

Méthode de Verlet (1967)¶

• Integration de l’ équation du mouvement en astrophysique et dynamique moléculaire
• Equation de Newton avec des forces $F$ dépendants de la position $X(t)$

$$ M \frac{d^2 X}{dt^2} = F(X(t)) \mbox{ avec } X(0)=X_0, \frac{dX}{dt}(0)=V_0$$

• Approximation centrée d’ordre 2 de la dérivée seconde $$ \frac{d^2 X}{dt^2} = \lim_{\Delta t \rightarrow \infty}{\frac{\Delta}{\Delta t}(\frac{\Delta X}{\Delta t})}\approx \frac{X_{n+1}-2X_n+X_{n-1}}{\Delta t^2}$$

Algorithme de Verlet

$$ X_{n+1} = 2X_n - X_{n-1} + \Delta t^2 A(X_n) \mbox{ avec } A(X)=M^{-1} F(X) $$

• erreur d’ordre 4, simplectique pour les systèmes Hamiltonien (conservatif)

• démarrage des itérations (ordre 3)

$$ X_1 = X_0 + \Delta t V_0 + \frac{\Delta t^2}{2} A(X_0) $$

exemple: trajectoire de la terre autour du soleil¶

loi de la gravitation universelle: planete masse $m$ distante de $(x,y)$ d’un soleil masse $M$: $$ \begin{cases} m \frac{d^2 x}{dt^2} = - \frac{G M m }{r^2} \frac{x}{r} \ m \frac{d^2 y}{dt^2} = - \frac{G M m }{r^2} \frac{y}{r} \end{cases}$$

Exercices¶

Swinging pendulum¶

On considère le mouvement de $N$ pendules simples de longueurs $L$ variables.

Écrire un programme python qui calcule le mouvement de ces pendules. Déterminez la répartition des longueurs $L$ permettant d’obtenir le type de mouvement ci-dessous avec $N=15$ pendules.

On comparera aussi le résultat avec les vidéos d’expériences trouvées sur le net.

pendule simple avec une corde souple¶

On considère les oscillations d’un pendule simple fixé à une corde souple de longueur maximum L et de masse négligeable par rapport au pendule. Lorsque la longueur de la corde est maximum, elle peut exercer uniquement une force de tension centripète sur la masse m du pendule. Le pendule est lancé avec une vitesse initiale, et commence à osciller. Si la tension dans le fil devient nulle, le pendule est alors soumis uniquement à la force de gravité et sa trajectoire devient parabolique. Lors que cette trajectoire intercepte le cercle de longueur l, un choc se produit puisque le fil atteint alors sa longueur maximum et exerce à nouveau une tension sur le pendule. Ce choc est un choc mou avec perte d’énergie, et l’on suppose que ce choc est telle que la vitesse tangentielle se conserve et la vitesse normale s’annule. Le pendule se met alors à osciller, mais avec des oscillations de plus faibles amplitudes.

Ecrire un programme python qui calcule la trajectoire du pendule.

problème à N corps¶

On considère un système planétaire de $N$ masses en interaction gravitationnelle. La force exercée par la masse $m_i$ sur la masse $m_j$ est donnée par la relation de Newton: $$ F_{i,j} = \frac{G m_i m_j}{r_{i,j}^2} \mbox{ avec } r_{i,j}=\mbox{distance entre les 2 masses}$$ $G$ est la constante de gravitation universelle. Pour $N=3$, on a le problème classique des 3 corps qui intrigue les physiciens depuis Newton.

Écrire un programme python calculant la trajectoire de $N$ corps se trouvant dans un même plan, modélisant un système planétaire où la première masse représente le soleil et les autres masses les planètes.

Quelle est la complexité de l’algorithme en fonction de N ?

Dans le cas du système à trois corps Soleil-Terre-Jupiter, retrouver le mouvement du système solaire.

unités: Dans le système d’unités astronomiques, i.e. la distance $L_t=1.49\,10^{11}\,m$ de la terre au soleil comme unité de longueur, la masse $M_t=6\,10^{24}\,kg$ de la terre comme unité de masse, et l’année terrestre moyenne $T_t=365,25\,jours$ comme unité de temps, on a $G=0.00012$. On a aussi la masse du Soleil $m_0=3.33\,10^5$, la masse de Jupiter $m_2=316.7$ et la distance Soleil Jupiter $r_2=5.2$.

Étudier le comportement du système lorsque la masse du soleil diminue.

Essayer de retrouver la solution périodique en huit découverte en 1993 par Christopher Moore.

pendule double: mouvement chaotic¶

On veut étudier le mouvement d’un double pendule par intégration des équations du mouvement.

Écrire un programme python pour calculer la position des 2 masses $m_1$ et $m_2$ au cours du temps.

On comparera la simulation aux vidéos d’expérience disponible sur le net.

<img src=”Double_Pendulum.png”, style=”width:200px;”/>

Évolution du volume d’une particule fluide¶

En suivant la trajectoire d’une particule fluide dans un champ de vitesse $U=[u(x,y),v(x,y)]$, on détermine la variation du volume $V(t) de la particule.

On veut vérifier la relation suivante: $$ \frac{1}{V}\frac{D V}{Dt} = div\,U = \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$$

Ecrire un programme python qui calcule la trajectoire $(x(t),y(t)$ d’une particule fluide: $$ \frac{d x}{dt}= u(x,y) \mbox{ et } \frac{d y}{dt}= v(x,y)$$ et détermine la variation de son volume $V(t)$

FIN¶