Note
Le schéma d’approximation de Newton est un algorithme «miraculeux» pour approcher numériquement les solutions d’équations de type \(f(x)=0\).
L’idée 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.
Sous des hypothèses minimales, 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.
Initialement conçu par Newton pour déterminer des racines de polynômes, ce schéma s’applique à approximer les zéros de fonctions définies sur un grand nombre d’espaces : \(\R\) ou \(\C\), \(\Z/p^n\Z\) ou des corps \(p\)-adiques, des anneaux de polynômes ou de séries formelles, des espaces de Banach…
Un nombre incalculable d’algorithmes doivent leur faible complexité à l’itération de Newton (en particulier l’algorithme de division), mais à l’instar du théorème du point fixe, ce schéma trouve aussi une place clé dans un certain nombre de démonstrations théoriques.
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).
Améliorer le code pour de grandes précisions en commençant l’algorithme avec une faible précision que l’on double à chaque étape. Comparer les temps de calcul pour calculer \(10^5\) chiffres de \(\sqrt{5}\).
sqrt_newton_incrprec(a):= {
}
def sqrt_newton_incrprec(a):
pa = a.precision()
x = sqrt(a*1.) # precision machine
px = x.precision()
while px < pa:
px <<= 1
x = RealNumber(x,min_prec = px)
x = (x+a/x)/2
return x
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.
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 \(anm\) 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.
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.