Commit 2fccdddb authored by Guillaume Garrigos's avatar Guillaume Garrigos
Browse files

TP1 premier jet

parent ecbe3a44
MIT License
Copyright (c) 2020 Guillaume Garrigos
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
# L3Optim TP
# Travaux Pratiques pour le cours d'Optimisation
test
Cetta page regroupe tous les documents nécéssaires pour réaliser les TPs du cours d'Optimisation (Université de Paris - L3 MFA/MI/IngeM).
Pour lancer le TP, il suffit de suivres les étapes suivantes :
1. Cliquer sur ce badge : [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/Guillaume-Garrigos/teaching_online_notebooks/main) afin de lancer une instance de Binder, qui vous permettra de coder dans votre navigateur, sans aucune installation prérquise. Cela va prendre pas mal de temps (entre1 minute) pour se lancer, soyez patients!
2. Une fois l'instance lancée, vous arriverez sur une page ressemblant à l'image ci-dessous, qui liste des fichiers.
3. Cliquez sur le fichier TPX.ipynb de votre choix, et cela ouvira le Notebook Jupyter correspondant!
Quelques informations importantes:
* Si vous restez inactifs pendant plus de 5 minutes, le serveur s'arrêtera :(
* Tout ce que vous ferez ne sera pas sauvegardé en ligne. Une fois votre TP terminé, si vous souhaitez sauvegarder votre travail, il vous faudra télécharger votre notebook (REF). Notez que parfois le document se télécharge avec l'extension ".json". Si c'est le cas, pensez à la remplacer par l'extension ".ipynb"
%% Cell type:markdown id: tags:
# TP1 : Méthode de Newton
%% Cell type:markdown id: tags:
# Avant de commencer : Introduction à Python et Jupyter
Python est un langage de programmation complet : il permet de définir des fonctions, faire du calcul numérique, traiter des bases de données, etc. Dans cette introduction, nous allons (re)visiter rapidement les caractéristiques principales de ce langage.
Jupyter est un environnement qui permet de coder dans plusieurs langages de programmation, et en particulier Python. Il permet d'écrire à la fois du **code** dans des cellules dédiées au code, et du **texte** (comme celui que vous lisez actuellement) dans des cellules dédiées au texte.
### Executer du code et des fonctions en Python
Pour exécuter le contenu d'une cellule, il suffit de se placer dans celle-ci et de taper Ctrl+Entrée ou Maj+Entrée. Vous pouvez essayer dans les cellules suivantes:
%% Cell type:code id: tags:
``` python
2+2
```
%% Cell type:code id: tags:
``` python
2*3
```
%% Cell type:code id: tags:
``` python
4**2
```
%% Cell type:code id: tags:
``` python
1+1
2+3
4*8
2+2
```
%% Cell type:markdown id: tags:
Comme vous pouvez le voir, par défaut Jupyter n'affiche que le résultat du *dernier* contenu de la cellule.
%% Cell type:markdown id: tags:
Python étant très versatile, et s'appliquant à des domaines différents, il a été choisi par ses concepteurs de ne pas proposer par défaut trop de fonctionalités. Si l'on veut avoir accès à des fonctions particulières, il faut au préalable **importer** les librairies qui leur correspondent. Ici on va importer les librairies `numpy` (calcul numérique et matriciel)
%% Cell type:code id: tags:
``` python
import numpy
```
%% Cell type:markdown id: tags:
Maintenant nous avons accès à de nouvelles fonctions! Prenons l'exemple de la fonction `sqrt` (*SQuare RooT*, pour racine carrée) : elle n'est pas implémentée par défault en python! Donc la cellule suivante va renvoyer un message d'erreur:
%% Cell type:code id: tags:
``` python
sqrt(4)
```
%% Cell type:markdown id: tags:
Par contre sqrt est une fonction de la librairie `numpy`. Pour dire à Python que l'on veut appeler une `fonction` depuis une certaine `librairie`, on utilise une commande de la forme `librairie.fonction()`. Dans notre cas on va donc appeler `numpy.sqrt()`, essayez vous même:
%% Cell type:code id: tags:
``` python
numpy.sqrt(4)
```
%% Cell type:markdown id: tags:
Par contre écrire `numpy` à chaque fois va être un peu laborieux. En général on choisit d'importer une librairie avec un surnom dont on peut se servir comme raccourci. Traditionellement `numpy` est raccourcie en `np`
%% Cell type:code id: tags:
``` python
import numpy as np
```
%% Cell type:markdown id: tags:
Cela permet d'écrire de manière plus concise:
%% Cell type:code id: tags:
``` python
np.sqrt(4)
```
%% Cell type:markdown id: tags:
### Boucles et logique en Python
La boucle **FOR** en Python s'écrit de la façon suivante:
%% Cell type:code id: tags:
``` python
for k in range(5):
print(k)
```
%% Cell type:markdown id: tags:
Il y a plusieurs commentaires à faire ici:
- la boucle affiche les numéros de 0 à 4, alors que l'on aurait pu s'attendre aux numéros de 1 à 5! Ceci est dû au fait (très important!) qu'en Python, les indices des listes, tableaux et autres "itérables" commencent à 0.
- Grosso modo, `range(5)` renvoie une liste `[0,1,2,3,4]`. Donc la première ligne veut simplement dire "pour k dans [0,1,2,3,4]".
- la première ligne se termine par `:` , il faut le mettre après chaque expression de ce type (voir plus bas)
- la deuxième ligne (intérieur de la boucle) doit être décalée! Cela peut être un décalage d'une tabulation, 2 espaces, 4 espaces, etc. Choisissez ce que vous voulez, l'important c'est d'être consistant. La communauté recommande d'utiliser une tabulation = 4 espaces.
%% Cell type:markdown id: tags:
La condition **IF** en Python s'écrit de la façon suivante :
%% Cell type:code id: tags:
``` python
if 2 > 0 and 6 > 4:
print("Ces deux propriétés sont vraies")
```
%% Cell type:code id: tags:
``` python
if 2 != 0 and 6 == 4:
print("Ces deux propriétés sont vraies")
else:
print("Une de ces propriétés au moins est fausse")
```
%% Cell type:markdown id: tags:
- Tout comme la boucle for, chaque "phrase" doit se terminer par un `:`. Cela concerne la ligne du `if`, mais aussi celle du `else`, si il y en a.
- On indente toujours les instructions à l'intérieur du `if` ou `else`
- On peut combiner plusieurs propositions syntaxiques avec `and` et `or`. On peut aussi utiliser des parenthèses pour les regrouper si besoin
On peut évidemment combiner tout ceci:
%% Cell type:code id: tags:
``` python
for k in range(5):
if k in [0,2,4]:
print("k est pair")
else:
print("k est impair")
```
%% Cell type:markdown id: tags:
### Fonctions en Python
La syntaxe pour définir des variables en python est assez simple:
%% Cell type:code id: tags:
``` python
a = 3 # on définit une variable a que l'on affecte de la valeur 3
```
%% Cell type:code id: tags:
``` python
b = a # on définit une nouvelle variable qui va prendre la même valeur que la première
```
%% Cell type:code id: tags:
``` python
b + 2 # on renvoie 3 + 2
```
%% Cell type:markdown id: tags:
Voyons maintenant comment définir une fonction en Python. Ici on veut une fonction qui renvoie le carré d'un nombre donné :
%% Cell type:code id: tags:
``` python
def carré(n):
c = n**2 # on élève n au carré
return c # on renvoie c
```
%% Cell type:code id: tags:
``` python
carré(4)
```
%% Cell type:markdown id: tags:
Comme vous pouvez le voir, la structure est assez simple:
- On commence la définition d'une fonction avec `def`, suivi du nom de la fonction, puis entre parethèses les variables d'entrée. Et on n'oublie pas le `:` !
- toutes les commandes à l'intérieur de la fonction sont indentées
- A la fin de la fonction, on renvoie une ou des variables de sortie avec `return`
Voyons un exemple avec plusieurs variables d'entrée et de sortie:
%% Cell type:code id: tags:
``` python
def division_euclidienne(a,b):
# Renvoie le résultat de la division de a par b
# càd un quotient q et un reste r qui vérifient
# a = b*q + r, avec 0 <= r < b
r = a % b # a modulo b
q = (a - r)/b
return q, r
```
%% Cell type:code id: tags:
``` python
division_euclidienne(13,3) # on peut afficher toutes les variables de sortie
```
%% Cell type:code id: tags:
``` python
q, r = division_euclidienne(13,3) # on peut mettre les variables de sortie directement dans des variables
```
%% Cell type:code id: tags:
``` python
q # on peut ensuite avoir accès à chaque variable de sortie
```
%% Cell type:markdown id: tags:
Vous savez maintenant à peu près tout ce dont vous avez besoin pour commencer le TP.
Voici quelques commandes python dont vous pourrez avoir besoin par la suite:
| | |
|-|-|
|`a*b`| Calcule le produit $a\times b$ |
| `a**b` | Calcule $a^b$ |
| `x=2` | Affecte $2$ à la variable `x` |
| `np.random.randn()` | Génère un nombre aléatoire selon la loi normale |
S'il s'avère que vous vous retrouviez coincés par la suite, n'hésitez pas à demander de l'aide.
D'autre part, si vous voulez faire quelque chose, mais ne savez pas *comment* le faire, Google est votre ami! Je dis ça de manière non-ironique : la documentation sur Python en ligne est énorme, tant en français qu'en anglais, et taper simplement votre question dans un moteur de recherche vous donnera immédiatement la réponse que vous chrechez!
%% Cell type:markdown id: tags:
# Le TP : Introduction
Dans ce TP, on s'intéresse au problème suivant : *comment calculer des racines carrées rapidement?* Plus précisément, étant donné un nombre réel $a>0$, comment calculer rapidement (une bonne approximation de) $\sqrt{a}$?
Cette question peut paraître étrange au premier abord, puisque on vient de voir que Python est parfaitement capable de calculer des racines carrées avec la fonction `sqrt` de `numpy`. Mais ici la question n'est pas tant de savoir comment obtenir une racine carrée, que de le faire **rapidement**.
En effet il existe des problèmes spécifiques qui requièrent de calculer des racines carrées **plusieurs millions de fois par seconde**, en temps réel. C'est par exemple le cas pour l'animation de synthèse (films Pixar), ou pour l'affichage de graphiques 3D dans les jeux vidéos, où une opération de base consiste à *normaliser* des vecteurs. Autrement dit, étant donné un vecteur $v \in \R^3$, à calculer sa direction normalisée $\frac{v}{\Vert v \Vert}$.
Or la norme euclidienne $\Vert v \Vert$ est égale à $\sqrt{\Vert v \Vert^2}$, où $\Vert v \Vert^2 = \sum_i v_i^2$ est facile à calculer.
%% Cell type:markdown id: tags:
Considérons un exemple concret : en 1999 est sorti un jeu vidéo en 3D (Quake III Arena), dont le code source a été rendu public en 2004. Des personnes ont alors exploré le code source du jeu, et se sont alors rendues compte que les développeurs du jeu avaient codé une fonction `Quake_isqrt` (*isqrt* pour *inverse square root*) dont le but est de calculer l'inverse de la racine carrée d'un nombre $a$, c'est-à-dire $\frac{1}{\sqrt{a}}$ :
| Quake III Arena | Le code original du jeu (traduit en Python) |
| --- | --- |
|<img src="https://miro.medium.com/max/600/1*vBFusX31Iwr5LYULQGolpQ.gif" width=200px> | <img src="https://i.imgur.com/GJLMmvD.png" width=600px> |
%% Cell type:markdown id: tags:
Au premier coup d'oeil, ce code parait un peu mystérieux (voire mystique, de l'aveu même de la personne qui a programmé, au vu des commentaires laissés). Essayer de comprendre ce qu'il se passe avec cette fonction va nous servir d'excuse pour étudier lors de ce TP un algorithme très important : la méthode de Newton.
%% Cell type:markdown id: tags:
## I. Méthode de Newton pour calculer $\sqrt{a}$
Etant donné un nombre positif $a>0$, on cherche à calculer ${\sqrt{a}}$ de manière approchée. Pour cela, on introduit la fonction $$f_a(x):= {x^2} - a,$$ qui vérifie la propriété: $$f_a(x)=0 \ \Leftrightarrow x = {\pm\sqrt{a}}.$$ Donc calculer cette racine carrée $\sqrt{a}$ est équivalent à trouver un **zéro** de la fonction $f_a$.
Pour cela, nous allons utiliser un [algorithme](https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Newton) que l'on doit à [Newton](https://fr.wikipedia.org/wiki/Isaac_Newton), dont le but est de trouver un zéro d'une fonction dérivable $f: \mathbb{R} \to \mathbb{R}$.
Cet algorithme s'énonce de la manière suivante: étant donné un point de départ $x_0 \in \mathbb{R}$, construire la suite $(x_n)_{n \in \mathbb{N}}$ définie par
\begin{equation}\label{D:Newton_algo}
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}. \tag{1}
\end{equation}
Alors, sous *certaines hypothèses*, la suite $x_n$ va converger, lorsque $n\to +\infty$, vers un point $x_\infty$ qui est un zéro de $f$. Autrement dit, un point qui vérifie $f(x_\infty)=0$.
Dans cette première partie du TP nous allons:
- Implémenter cette méthode,
- Essayer de comprendre (expérimentalement) sous quelles conditions elle fonctionne.
%% Cell type:markdown id: tags:
**0)** On souhaite appliquer la méthode de Newton à la fonction $f_a$ introduite plus haut. pour cela, vous devez calculer, à la main, l'expression de récurrence (1).
%% Cell type:markdown id: tags:
**1)** Programmer une fonction `Newton_sqrt` sous la forme `Newton_sqrt(a,x0,nmax)` telle que:
- les variables d'entrée sont
- `a`, le nombre dont on veut calculer la racine carrée
- `x0`, le point de départ de l'algorithme
- `nmax`, le nombre d'itérations maximal que l'on souhaite faire
- la variable de sortie est l'itéré $x_{nmax}$ de la méthode de Newton.
- la fonction implémente la méthode de Newton avec une boucle for sur le nombre des itérés.
Vous pourrez au besoin vous référer à la documentation en ligne pour connaitre la syntaxe d'une [fonction](https://python.doctor/page-apprendre-creer-fonction-en-python) et d'une [boucle FOR](https://python.doctor/page-apprendre-boucles-python-loop).
%% Cell type:code id: tags:
``` python
def Newton_sqrt(a, x0, nmax):
return x
```
%% Cell type:markdown id: tags:
**2)** Supposons qu'on veuille tester la méthode avec $a=4$. Si tout marche bien, quelle est (au signe près) la valeur que la fonction est supposée renvoyer?
Tester que votre méthode marche bien pour `a=4`, avec `nmax=100`, et divers points initiaux de votre choix. Est-ce que cela marche bien? Est-ce toujours le cas, quelque soit le point initial?
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
**3)** Tester la méthode pour `a=4`, avec `x0=10`, et cette fois-ci vous testerez diverses valeurs pour `nmax`. Faut-il faire beaucoup d'itérations pour obtenir un résultat raisonnable?
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
**4)** Au vu des questions précédentes, il semble que lorsque l'algorithme converge, il converge **vite**. Donc il n'y a pas besoin de faire beaucoup d'itérations pour obtenir un résultat satisfaisant.
Cependant, ce nombre d'itérations peut fluctuer, et il n'est pas non plus une bonne idée de fixer une fois pour toutes un certain nombre d'itérations : si `nmax` rest trop grand on perd notre temps, si il est trop petit on n'obtient pas une très bonne solution. Il serait donc préférable d'avoir un critère qui *arrête* l'algorithme lorsqu'on est *proche de la solution*, afin d'économiser du temps de calcul.
Vous allez donc reprendre le code de la fonction `Newton_sqrt`, et le modifier de telle façon:
- la fonction va prendre un argument d'entrée supplémentaire `tol`. Ceci désigne une certaine tolérance que l'on exige pour la solution. Plus exactement, on veut que $\vert f_a(x_n)\vert < tol$.
- La boucle qui fait tourner l'algorithme ne sera plus une boucle FOR mais une boucle WHILE. Vous pourrez au besoin vous référer à la documentation en ligne pour connaitre la syntaxe d'une [boucle WHILE](https://python.doctor/page-apprendre-boucles-python-loop).
- On **arrêtera** l'algorithme dès lors que $n>nmax$ **ou** $\vert f(x_n)\vert < tol$. C'est à vous de déterminer comment écrire proprement la boucle WHILE, notamment quelle est la condition à vérifier pour pouvoir continuer à entrer dans la boucle.
- Juste avant de renvoyer la solution, la fonction devra afficher un message, à l'aide de la fonction `print`, pour indiquer le nombre d'itérations qui ont été nécéssaires. Si vous ne savez pas comment faire cela, je vous conseille de tester à part comment faire. Vous pourrez également consulter la [doc](https://python.sdv.univ-paris-diderot.fr/03_affichage/).
%% Cell type:code id: tags:
``` python
def Newton_sqrt(a, x0, nmax, tol):
return x
```
%% Cell type:markdown id: tags:
**5)** Tester que la fonction fonctionne, avec `a=4`, `nmax=100`, `x0=4` et `tol=1e-2`. On fera varier la valeur de `tol` et on regardera son impact sur le nombre d'itérations nécéssaires.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
Pas mal non?
%% Cell type:markdown id: tags:
## II. Méthode de Newton pour $1/\sqrt{a}$
Etant donné un nombre positif $a>0$, on veut calculer $\frac{1}{\sqrt{a}}$ de manière approchée. La raison pour cela est que dans le problème d'images 3D mentionné précédemment, on voulait normaliser des vecteurs $v$, et donc calculer $\frac{v}{\sqrt{a}}$ avec $a = \sum_i v_i^2$.
On pourrait légitimement se demander **pourquoi** faire ça alors qu'on a déjà trouvé une méthode qui calcule $\sqrt{a}$ !? Il suffit de faire `1/Newton_sqrt(a,x0,nmax,tol)` non? En fait, c'est un peu plus compliqué que ça:
1. La division est une opération qui coûte *cher* en comparaison avec l'addition et la multiplication. Bon, si vous faites une ou deux divisions vous n'allez pas vraiement vous en rendre compte, mais quand vous devez en faire plusieurs millions par seconde, la différence commence à se faire sentir...
2. Calculer $\frac{1}{\sqrt{a}}$ à partir de $\sqrt{a}$ requiert une division *évidente*. Mais cela ne s'arrête pas là! Par exemple, regardez votre code pour la fonction `Newton_sqrt`. Il y a une division à chque étape de la boucle!
3. Maintenant, si on repart de zéro et qu'on cherche à calculer directement $\frac{1}{\sqrt{a}}$ avec la méthode de Newton, on va voir qu'il n'y a **pas besoin de divisions** (c'est magique).
%% Cell type:markdown id: tags:
**0)** Revenons donc à notre problème, qui est de calculer $\frac{1}{\sqrt{a}}$ avec la méthode de Newton. Pour cela il nous faut d'abord trouver une fonction qui s'annule seulement en $\frac{1}{\sqrt{a}}$.
On introduit donc la fonction $$f_a(x):= \frac{1}{x^2} - a,$$ qui vérifie la propriété $$f_a(x)=0 \Leftrightarrow x = \pm \frac{1}{\sqrt{a}}.$$ Donc calculer cette inverse de racine carrée est équivalent à trouver un zéro de la fonction $f_a$.
- Ecrire sur papier la relation de récurence qui définit l'algorithme de Newton appliqué à cette fonction. Vérifier qu'il n'y a pas besoin de faire de divisions.
- Retournez voir la fonction `Quake_isqrt` mentionnée au début du TP, et regardez ses deux dernières lignes. Comprenez-vous maintenant ce que font ces lignes?
%% Cell type:markdown id: tags:
**1)** En se basant sur le squelette de `Newton_sqrt`, programmer une fonction `Newton_isqrt` (*isqrt* pour *inverse square root*) sous la forme `Newton_isqrt(a,x0,nmax,tol)` qui :
- prend en entrée une valeur `a`, un point de départ `x0`, un nombre d'itérations maximal `nmax` et une tolérance `tol`
- applique la méthode de Newton calculée à la question précédente, avec une boucle WHILE
- arrête l'algorithme lorsque $n>nmax$ **ou** $\vert f_a(x_n)\vert < tol$
- indique combien d'itérations ont été effectuées
- renvoie l'approximation `x` calculée par la méthode de Newton
%% Cell type:code id: tags:
``` python
def Newton_isqrt(a, x0, nmax, tol):
return x
```
%% Cell type:markdown id: tags:
**2.1)** Supposons qu'on veuille tester la méthode avec $a=4$. Si tout marche bien, quelle est (au signe près) la valeur que la fonction est supposée renvoyer? Soyez attentifs!
Vérifiez que la méthode marche avec `a=2`, `tol = 10**-8`, `nmax = 100` et `x0 = 1`.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
**2.2)** On s'intéresse maintenant au comportement de la méthode en fonction de l'initialisation `x0`. En gardant les mêmes paramètres que précédemment, on testera des valeurs comme 0.1, 0.2, ... 1.5. Que constatez-vous?
Si vous voulez comprendre ce qui se passe, il vous faudrait savoir comment évolue $x_n$ lorsque l'algorithme tourne. Vous pouvez par exemple temporairement ajouter une ligne `print(x)` dans la boucle WHILE (à adapter selon vos notations).
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
**3)** Au vu de la question précédente, on voit que pour cette fonction l'initialisation est *très* importante. En fait la théorie garantit, de manière générale, que l'algorithme de Newton ne converge vers une solution que s'il est initialisé *suffisament proche* d'une solution.
- Evidemment, il y a des cas favorables où tout se passe bien. Regardez par exemple la partie I du TP où il n'y avait aucun problème!
- Il est impossible de savoir à l'avance exactement quels sont les `x0` qui vont permettre de converger, et ceux qui vont diverger. Et même dans des cas simples, la répartition entre "bons" points initiaux `x0` et "mauvais" peut faire apparaitre des [fractales](https://fr.wikipedia.org/wiki/Fractale_de_Newton) ...
Du coup, retournons voir ce code issu du jeu vidéo Quake:
| |
| --- |
|<img src="https://i.imgur.com/GJLMmvD.png" width=600px> |
On comprend bien que la dernière ligne (plus celle commentée) correspondent à une étape de l'algorithme de Newton appliquée au nombre `x`. Donc on peut supposer que toutes les lignes bizarres qui précèdent servent à calculer une initialisation de l'algorithme à partir de `a`. Une initialisation tellement bonne qu'on a besoin que d'une seule itération de la méthode de Newton! Deux commentaires à ce propos:
- Nous n'avons pas pour objectif ici de comprendre comment se fait cette initialisation. Ce n'est pas franchement très important, ni facile à comprendre. Si vous êtes curieux/curieuse, vous pourrez aller jeter un oeil plus tard à [ceci](https://fr.wikipedia.org/wiki/Racine_carr%C3%A9e_inverse_rapide) ou [cela](https://www.blog-libre.org/2015/12/15/la-mysterieuse-fonction-racine-carree-inverse-de-quake-3/)).
- Par contre il serait intéressant de vérifier que cela marche bien! C'est ce que l'on fait dans la prochaine question:
%% Cell type:markdown id: tags:
**4)** Vous trouverez ci-dessous une fonction implémentant l'initialisation *magique* utilisée par les développeurs. Utiliser cette fonction, combinée à votre fonction `Newton_isqrt`, dans les mêmes conditions qu'à la question 2). Que constatez-vous?
En restreignant le nombre d'itérations maximal à 1, que pensez-vous de la solution ainsi obtenue?
%% Cell type:code id: tags:
``` python
def Quake_initial(a):
x0 = np.float32(a)
i = x0.view(np.int32)
i = np.int32(0x5f3759df) - np.int32(i >> 1)
x0 = i.view(np.float32)
return float(x0)
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
## III. Une approche plus simple : une bonne vieille méthode de dichotomie
%% Cell type:markdown id: tags:
Revenons au problème de calculer $\sqrt{a}$. Au lieu de partir sur des choses aussi compliquées que la méthode de Newton, on pourrait se dire qu'il serait bien plus simple de chercher la solution dans un certain intervalle, en y mettant une grille, ou plus efficacement, en faisant une recherche de *dichotomie*.
Ici encore, on va rechercher les zéros de la fonction $f_a(x)=x^2-a$.
La méthode de dichotomie consiste à se donner un intervalle $[m,M]$ dans lequel on recherche un $x$ tel que $f(x)=0$. Afin d'être sûrs qu'il existe une solution sur cet intervalle, on demandera que $f(m)$ et $f(M)$ soient de signe opposé (le théorème des valeurs intérmédiaires nous garantissant alors l'existence d'une solution). La méthode peut s'énoncer comme suit, en fonction d'un paramètre de tolérance:
```
Vérifier que f(m)f(M) < 0
Tant que tol < M-m # tant que notre intervalle est trop gros
x = (M+m)/2 # on prend le point milieu
Si f(m)f(x) < 0, # si il y a encore une solution dans l'intervalle [m,x]
alors M <-- x, # on rétrécit l'intervalle à [m,x]
sinon m <-- x # sinon forcément la solution est de l'autre coté donc on rétrécit à l'intervalle [x,M]
```
%% Cell type:markdown id: tags:
**1)** Définir une fonction `dichotomie(a,m,M,tol)` pour résoudre $f_a(x)=0$ sur un intervalle [m,M], et qui renvoie une solution approchée `x`. La fonction affichera également le nombre d'itérations nécessaires pour le trouver.
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
def dichotomie(a, m, M, tol):
return x
```
%% Cell type:markdown id: tags:
**2)** Tester cette fonction `dichotomie` avec $a=4$, $m=0,\ M=4$ et avec différentes valeurs de la tolérance $10^{-2},10^{-3},10^{-4},10^{-5}$... Noter le nombre d'itérations $n$ et contrôler le résultat.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
**3)** Comparer les méthodes de Newton et Dichotomie, et déterminer laquelle est la plus rapide?
%% Cell type:code id: tags:
``` python
```
name: example-environment
channels:
- conda-forge
dependencies:
- numpy
- matplotlib
- addict
- ipympl
- ipywidgets
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import rcParams
from mpl_toolkits.mplot3d import Axes3D
import ipywidgets as widgets
from ipywidgets import AppLayout, FloatSlider
from ipywidgets import GridspecLayout
from matplotlib import rc
from copy import deepcopy
try:
from addict import Dict
except: