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 quadratique, c’est-à-dire que la précision obtenue (le nombre de décimales correctes) double à chaque étape. A comparer avec la convergence linéaire de 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)
def sqrt_newton(a):
x1, x0 = 1, 0
while x0 != x1:
x0 = x1
x1 = (x0+a/x0)/2
return x0
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).
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
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]\).
On considère désormais l’anneau \(\Z/m\Z\), et même la suite d’anneaux \(\Z/m^k\Z\).
Soit \(P\in\Z[X]\) un polynôme, la méthode de Newton permet de «remonter» une racine de \(P\) modulo \(m\) en une racine modulo \(m^{2^k}\) (sous réserve de pouvoir calculer l’inverse de \(P'(a)\) modulo \(m\)).
On considère le polynôme \(P(X):=X^3-3*X^2+3*X-10^9-1\). Vérifier que \(P(0)\equiv 0\bmod 7\), et que \(P'(0)\) est inversible modulo \(7\).
P(x):=x^3-3*x^2+3*x-10^9-1;
P(0)%7;
P'(0)%7;
Calculer \(a_1\in\Z\) tel que \(a_1\equiv a-P'(a)^{-1}*P(a) \bmod 49\) (utiliser
la fonction invmod pour calculer un inverse modulaire).
a := 0;
a1 := a - invmod(P'(a),49)*P(a) % 49;
P(a1);
Continuant l’itération modulo \(7^4\) et \(7^8\). Que se passe-t-il ? Trouver une racine de \(P\) dans \(\Z\).
a2 := a1 - invmod(P'(a1),7^4)*P(a1) % 7^4;
P(a2);
a3 := a2 - invmod(P'(a2),7^8)*P(a2) % 7^8
La suite semble stationner : en fait \(1001\) est une racine dans \(\Z\):
P(1001)
Écrire une procédure hensel(P,a,m,k) qui à partir d’un polynôme
\(P\in\Z[X]\) et d’entiers \(a,m\) tels que \(P(a)\equiv 0\bmod m\) et
\(\pgcd(P'(a),m)=1\), renvoie un entier \(b\) tel que \(P(b)\equiv 0\bmod m^{2^k}\).
1 2 3 4 5 6 7 8 9 10 11 12 13 | hensel(P,a,m,k):= {
local Pa, Pam, dPa, ah;
if(k==0) return a;
Pa := P(a);
print(Pa,a,m,k);
if(irem(Pa,m)!=0) return( "mauvais argument : P(a)!=0 mod m");
Pam := iquo(Pa,m);
dPa := diff(P)(a);
if(gcd(dPa,m)!=1) return( "mauvais argument : P'(a) non inv mod m");
Pam := iquo(Pa,m);
ah := irem(a - m*Pam*((dPa mod m)^(-1) %0),m^2);
return hensel(P,ah,m^2,k-1);
}
|
Montrer que \(-1\) et \(4\) sont les deux autres racines de \(P\) modulo \(7\), et que l’on peut leur appliquer la méthode de Hensel. Que se passe-t-il à présent ?
En utilisant le schéma de calcul d’inverse (et non plus la fonction
invmod), calculer l’inverse de \(94\) modulo \(6561\). Combien
d’itérations sont nécessaires ?
6561 vaut 3^8 = 3^(2^3) et 94 = 1 mod 3. On calcule l’itération d’inversion 3 fois:
a := 1;
for ii from 1 to 3 do
a := irem(2*a - 94*a**2,6561) ;
end_for;
Essayer d’améliorer la fonction hensel pour que l’étape d’inversion
modulo \(m^{2^k}\) soit elle-même réalisée par le schéma de Newton.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | Hensel := proc(P,a,m,k)
local Pa, inv_dPa, dP;
dP := diff(P);
// il faut calculer le premier inverse
inv_dPa := bezout_entiers(dP(a),m)[0];
l := 0;
pow2l := 1;
while(pow2l<k) {
m := m^2;
a := irem(a - P(a) * inv_dPa,m);
inv_dPa := irem(2*inv_dPa-dP(a)*inv_dPa^2,m);
};
return a;
end proc;
|
Utiliser la méthode de Newton pour inverser une matrice \(A\in M_{5,5}(\R)\).
L’itération est \(X\mapsto 2X-XAX\), et on utilisera comme valeur initiale \(A^T/(\norm{A}_1\norm{A}_{\infty})\).
Note
À l’inverse de tous les cas précédents, ce n’est pas du tout une bonne méthode d’inversion matricielle, mais elle permet d’augmenter la précision sur une solution initiale.
Implanter aussi une fonction inverse_ser qui calcul l’inverse d’une
série entière de terme constant 1 via l’itération de Newton.
À 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.