Note
Le schéma d’approximation de Newton est un algorithme pour calculer numériquement les solutions d’équations de type \(f(x)=0\).
Le principe est de partir d’une solution approchée \(x_0\), et de l’améliorer en résolvant l’équation linéarisée \(f(x_0)+df(x_0).(x-x_0) = 0\), et de répéter le procédé avec la solution \(x_1\) obtenue.
En appliquant la formule de Taylor si \(f\) est de class \(\cC^2\) au voisinage de \(x\) et telle que \(f'\) ne s’annule pas, on démontre que la convergence obtenue est doublement exponentielle, c’est-à-dire que la précision obtenue (le nombre de décimales correctes) double à chaque étape : c’est exponentiellement plus rapide que la méthode du point fixe.
Mots-clefs xcas: evalf, diff, Digits, normal.
On considère la fonction \(f(x):=x^2-2\). Faites calculer la dérivée de \(f\), puis déterminer l’itération de Newton associée à \(f\)
f(x):=x^2-2
diff(f)
Nf(x) := normal(x-f(x)/diff(f)(x)); Nf(x)
On considère la suite récurrente \((u_n)\) donnée par l’itération ci-dessus et
l’initialiation \(u_0=1\).
Entrer u:=1.; (avec un . pour avoir le nombre flottant)
et appliquer l’itération à u plusieurs fois :
vérifier que l’on obtient bien une valeur approchée de \(\sqrt{2}\).
u:=1.
u:=Nf(u)
u:=Nf(u)
u:=Nf(u)
u:=Nf(u) // le resultat ne change plus
evalf(sqrt(2))
Écrire une fonction Newton_iteration(f) qui calcule l’expression
d’itération de Newton associée à une fonction \(f\) quelconque.
Newton_iteration(f):=normal(x-f(x)/diff(f)(x))
Comment travailler avec \(k\) chiffres significatifs ? Fixer la précision à 150 chiffres significatifs.
Digits:=150;
Écrire la suite récurrente donnée par la méthode de Newton qui permet de calculer \(\sqrt{a}\).
On calcule une racine de \(f(x)=x^2-a=0\). L’itération est donnée par
Newton_iteration(x->x^2-a)
Vérifier la convergence de la suite obtenue, en prenant toujours \(1\) comme valeur initiale.
On pourra par exemple écrire une procédure
sqrt_newton(a,n) qui effectue \(n\) étapes du calcul et affiche à chaque étape la
valeur de \(x_{n+1}-x_n\) (ou mieux, \(\log_{2}(\abs{x_{n+1}-x_n})\).
1 2 3 4 5 6 7 8 9 10 11 | sqrt_newton(a,n):={
local(x); local(x1); local(ii);
x1 := 1;
for(ii:=0;ii<n;ii++) {
x := x1;
x1 := (x+a/x)/2;
print(evalf(logb(abs(x-x1),10),10));
}
}
sqrt_newton(5.,7)
|
Prédire le nombre d’itérations nécessaires pour obtenir \(10^5\) décimales correctes de \(\sqrt{5}\) (noter le nombre d’itérations initiales en partant de \(x_0=1\) pour avoir un premier chiffre correct). Le vérifier expérimentalement.
On a besoin de deux itérations initiales pour avoir un premier chiffre correct. La précision double à chaque itération, la suite converge exponentiellement vite en nombre de chiffres corrects.
Il faut donc \(2+\log_2(10^5)\approx18\) itérations pour obtenir \(10^5\) décimales de \(\sqrt{5}\).
evalf(logb(10^5,2),2)+2
Digits:=10^5
evalf(log10(sqrt(5.) - sqrt_newton(5.,18)),10)
evalf(log10(sqrt(5.) - sqrt_newton(5.,17)),10)
Écrire une fonction sqrt_init(a,x0) qui calculer \(\sqrt{a}\) à l’aide de la
méthode de Newton en partant de la valeur x0. Pour quelles valeurs initiales
a-t-on convergence ? (repasser à \(150\) chiffres significatifs).
[option C] Avec l’algorithme ci-dessus, quelle est la complexité du calcul de \(\sqrt{a}\) à précision \(b\) ? Proposer une amélioration pour éviter de calculer à chaque étape des décimales inutiles, et en donner la complexité.
On suppose \(f\) de classe \(\cC^2\) sur un intervalle \(I\) contenant une solution \(x\) de l’équation \(f(x)=0\), et telle que \(0<\inf\abs{f'(t)}=m\) et \(M = \sup\abs{f''(t)}<\infty\) sur \(I\).
Montrer que pour toute valeur \(x_0\in I\), on a
On écrit la formule de Taylor exacte sur l’intervalle \(]x_0,x[\subset I\)
ainsi que l’équation définissant \(x_1\)
par différence on obtient exactement
En déduire un critère de convergence de la méthode, ainsi qu’une majoration de \(-\log\abs{x_n-x}\).
Notons \(C=\frac{M}{2m}\). Tant que chaque \(x_k\in I\), on a par récurrence
En particulier, si \(\abs{x_0-x}\leq 1/C\), l’itération est contractante donc la suite converge. Dans ce cas le nombre de chiffres corrects à l’itération \(n\) satisfait
On appelle fonction de Lambert la fonction \(W:y\mapsto W(y)\) définie par l’équation implicite
Donner un schéma de Newton pour calculer \(W(y)\).
g(x) := x*exp(x)-y
Newton_iteration(g)
On trouve le schéma
L’implanter. On pourra effectuer un développement asymptotique de la solution pour choisir une valeur initiale.
Déterminer un schéma de Newton permettant de calculer l’inverse \(\frac 1a\) d’un nombre \(a\) sans effectuer de division.
Pour inverser \(a\), \(f(x)=ax\) ne donne rien, en revanche \(f(x)=a-1/x\) donne l’itération \(x_{n+1} = 2*x_n-a*x_n^2\) qui ne contient pas de division.
Note
C’est ce schéma qui permet de réduire le problème de la division à celui de la multiplication. C’est ainsi que calculent tous les ordinateurs.
Écrire une fonction inverse_newton(a,a0,k) qui effectue \(k\) itération
du schéma en partant de \(a_0\). Examiner la convergence sur quelques exemples
ainsi que l’influence de l’initialisation.
1 2 3 4 5 6 7 8 9 10 11 12 | inverse_newton(a,a0,k):={
local(x); local(ii);
x := a0;
for(ii:=0;ii<k;ii++){
x := 2*x-x**2*a;
}
return(x);
}
inverse_newton(evalf(Pi),1,3); // diverge
inverse_newton(evalf(Pi),0.3,3); // converge
inverse_newton(evalf(Pi),0.5,3); // converge doucement au debut
|
On considère un nombre réel \(a\). Quel est le domaine de valeurs initiales pour lequel le schéma converge ?
On pourra utiliser la fonction plotseq et étudier analytiquement la suite récurrente.
common_digits(a,b):={if(a==b){Digits}else{max(0,min(Digits,-log10(abs(a-b))))}};
R := evalf(1/Pi):;
plot( common_digits(inverse_newton(evalf(Pi),x,7),R),x=0..3,nstep=100)
Si l’on fait l’étude de la suite \(u_{n+1}=f(u_n)=2u_n-au_n^2\), on remarque que \(f\) croise la droite \(y=x\) en \(0\) et \(1/a\), avec des dérivées valant respectivement \(2\) et \(0\) : le premier point fixe est répulsif vers le second, attractif. En regardant les variations, le domaine de convergence est \([0,2/a]\).
À partir de la fonction \(-\log_{10}(\abs{a-b})\), définir une fonction dist(a,b)
définie sur \(\C\times\C\) à valeurs dans
\([0,\mathtt{Digits}]\).
dist(a,b):={if(a==b){Digits}else{max(0,min(Digits,-log10(abs(a-b))))}};
On considère le polynôme \(X^5+1\). Écrire une fonction
newton5(x0,k) qui effectue \(k\) itérations de Newton en partant de x0,
et renvoie dist(x_k,x_{k-1}).
L’itération est \(x\mapsto \frac{4x-x^{-4}}5\).
newton5(x0,k):={
local(x); local(y); local(ii);
y:=x0;
for(ii:=0;ii<k;ii++) {
x := y;
y := (4*x-x^(-4))/5;
};
return(dist(x,y));
}
En utilisant la fonction plotdensity, représenter la convergence en
fonction de la valeur initiale \(x_0\) dans le plan complexe (utiliser par
exemple \(7\) itérations).
Note
si l’on représente graphiquement la convergence d’un schéma de Newton en fonction de l’initialisation dans \(\C\), on obtient la plupart du temps de très jolies fractales. Les ensembles de Mandelbrot ou de Julia peuvent être obtenus de cette manière.