Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download

Jupyter notebook Semaine 2 - Methode de Newton/newton-polynomes.ipynb

Views: 163
Kernel: SageMath 7.3

Comprendre la méthode de Newton sur l'exemple des polynômes.

Le but de ce bloc-notes est de comprendre la méthode de Newton-Raphson d'abord telle qu'elle a été expliquée par Newton, puis avec les notations de Raphson, mais toujours sur l'exemple des polynômes.

La méthode de Newton-Raphson est une méthode pour approcher les racines d'une fonctions dérivables - c'est à dire les arguments sur lesquels la fonction s'annule. Vous en saurez plus très vite.

Crédits: Ce bloc-notes s'inspire de l'article «La méthode de Newton et son histoire» d'André Bonnet, paru dans un bulletin de l'APMEP et disponible en libre téléchargement ici : http://www.apmep.fr/IMG/pdf/APMEP_article_BV_7_C_A.pdf Cependant les seuls extraits de texte recopiés tels quel depuis cet article, sont ceux des textes historiques. (Il n'est pas clair par quelle licence de droit d'auteur est régi cet article.)

1. La méthode telle qu'elle est décrite par Newton

Dans les deux paragraphes suivants, vous avez des extraits du texte de Newton (en anglais) où ce qu'on appelle aujourd'hui la méthode de Newton apparaît pour la première fois. Ce texte a été publié après sa mort par Colson sous le titre "The method of fluxions and infinite series". Vous pouvez le télécharger ici : http://www.e-rara.ch/zut/doi/10.3931/e-rara-10833 Ce même texte traduit en français par Buffon est lisible ici: http://gallica.bnf.fr/ark:/12148/bpt6k62411f.r=La méthode des fluxions et des suites infinies?rk=21459;2

«§20. Let this equation y³−2y−5=0 be proposed to be resolved, and let 2 be a number (any how found) which differs from the true Root less than by a tenth part of itself. Then I make 2+p=y, and substitute 2+p for the given Equation, by which is produced a new Equation p³+6p²+10p−1=0, whose Root is to be sought for, that it maybe added to the Quote. Thus rejecting p³+6p² because of its smallness, the remaining Equation 10p−1=0, or p = 0.1 will approach very near to the truth. Therefore I write this in the Quote and suppose 0,1+q=p, and substitute this fictious value of P as before, which produces q³+6,3q²+11,23q+0,061=0. And since 11,23q+0,061=0 is near the truth, or q = −0.054 nearly, (that is dividing 0,061 by 11,23, till so many Figures arise as there are places between the first figures of this, and of the principal Quote exclusively, as there are two places between 2 and 0,005) I write −0,0054 in the lower part of the quote, as being negative; and supposing −0,0054+r=q, I substitute this as before. And thus I continue the Operation as far as I please, in the manner of the following diagram: ┌────────────────┬─────────────────────────────┐ │y³−2y−5=0 │+2,10000000 │ │ │−0,00544852 │ │ │─────────── │ │ │+2,09455148, &c=y │ ├────────────────┼─────────────────────────────┤ │2+p=y. +y³│+ 8 + 12p + 6p² + p³ │ │ −2y │− 4 −2p │ │ −5 │− 5 │ ├────────────────┼─────────────────────────────┤ │ The sum │− 1 + 10p + 6p² + p³ │ ├────────────────┼─────────────────────────────┤ │0,1+q=p. +p³│+ 0,001 + 0,03q + 0,3q² + q³│ │ +6p²│+ 0,06 + 1,2 q + 6 q² │ │ +10p │+ 1 + 10 q │ │ −1 │− 1 │ ├────────────────┼─────────────────────────────┤ │ The sum │+ 0,061 + 11,23q + 6,3q² + q³│ ├────────────────┼─────────────────────────────┤ │−0,0054+r=q. +q³│[…] │ │ +6,3 q²│[…] │ │ +11,23 q │[…] │ │ +0,061 │[…] │ ├────────────────┼─────────────────────────────┤ │ The sum │0,0005416+11,162r […] │ ├────────────────┼─────────────────────────────┤ │−0,00004852+s=r.│[…] │ └────────────────┴─────────────────────────────┘ »

«§21. But the Work may be much abbreviated towards the end by this Method, especially in Equations of many Dimensions. Having first determined how far you intend to extract the Root, count so many places after the first Figure of the Coefficient of the last Term but one, of the Equations that result on the right side of the Diagram, as there remain places to be filled up in the Quote, and reject the Decimals that follow. But in the last Term the Decimals may be neglected, after so many more places as are the decimal places that are filled up in the Quote. And in the antepenultimate Term reject all that are after fo many fewer places. And so on, by proceeding Arithmetically, according to that Interval of places: Or, which is the same thing, you may cut off every where so many Figures as in the penultimate Term, so that their lowest places may be in Arithmetical Progression, according to the Series of the Terms, or are to be supposed to be supplied with Cyphers, when it happens otherwise. Thus in the present Example, if I desired to continue the Quote no farther than to the eighth place of Decimals, when I substituted 0,0054+r for q, where four decimal places are compleated in the Quote, and as many remain to be compleated, I might have omitted the Figures in the five inferior places, which therefore I have marked or cancelled by little Lines drawn through them ; and indeed I might also have omitted the first Term r³, although its Coefficient be 0,99999. Those Figures therefore being expunged, for the following Operation there arises the Sum 0,0005416+11,162r, which by Division, continued as far as the Term prescribed, gives 0,00004852 for r, which compleats the Quote to the Period required. Then subtracting the negative part of the Quote from the affirmative part, there arises 2,09455148 for the Root of the propofed Equation. »

Vous allez maintenant refaire l'exemple de Newton.

Déclarations et affectations

Commencons par déclarer les objets utiles.

Nous allons chercher une racine approchée du polynôme y³−2y−5. Il nous faudra donc déclarer ce polynôme. Pour mettre en application les concepts appris cette semaine, je vais vous demander de bien préciser son Anneau de coefficients.

Pour rappel, lorsqu'on veut que x represente le polynôme x à coefficients dans l'anneau RR = RealField(), qui est une des implementations à virgule flottante du corps des nombres réels dans Sage, on écrit x = polygen(RR,'x'). (C'est à dire qu'ensuite x est vu comme un élément de l'anneau de polynômes RR.)

y = polygen(QQ,'y')
p = polygen(QQ,'p')
q = polygen(QQ,'q')
r = polygen(QQ,'r')
s = polygen(QQ,'s')
type (y)
<type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>
P = y^3 - 2*y -5
type (P)
<type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>
S = P.change_ring(QQbar)
R = P.change_ring(AA)
P.factor()
y^3 - 2*y - 5
S.factor()
(y - 2.094551481542327?) * (y + 1.047275740771164? - 1.135939889088929?*I) * (y + 1.047275740771164? + 1.135939889088929?*I)
R.factor()
(y - 2.094551481542327?) * (y^2 + 2.094551481542327?*y + 2.387145908831156?)
//--Remarque: - P ne sest pas factorisé --// - S a été factorisé avec des coefficients rationnels --// - R a été factorisé avec des coefficients réels --//
P.roots()
[]
S.roots()
[(2.094551481542327?, 1), (-1.047275740771164? - 1.135939889088929?*I, 1), (-1.047275740771164? + 1.135939889088929?*I, 1)]
R.roots()
[(2.094551481542327?, 1)]
//-- Remarque: - P na pas donné de racines --// //-- - S a donné des racines rationnelles --// //-- - R a donné des racines réelles --//
y0=2
P(y0+1/10)
61/1000
P(y0-1/10)
-1941/1000
//-- Effectivement les nombres sont de signes différents --//
P1= P(y0+p)
P1
p^3 + 6*p^2 + 10*p - 1
type (P)
<type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>
type (P1)
<type 'sage.rings.polynomial.polynomial_rational_flint.Polynomial_rational_flint'>
//-- type (P) est égal à type (P1) --//
P1.degree()
3
def termes_grand_degre(polynome, d=2): pp=polynome.parent().0 pd=polynome.degree() res=0 for i in range (d, pd+1): res = res + polynome[i]*pp^i return res
termes_grand_degre(P)
y^3
P1_diff=termes_grand_degre(P1)
P1_diff
p^3 + 6*p^2
P1_diff(0.1)
0.0610000000000000
P1_lin=P1-P1_diff
P1_lin
10*p - 1
p0=P1_lin.roots()[0][0]
p0
1/10
P1(p0-1/100)
-50671/1000000
P1(p0+1/100)
173931/1000000
P2=P1(p0+q)
P2
q^3 + 63/10*q^2 + 1123/100*q + 61/1000
P2_diff=termes_grand_degre(P2)
P2_diff
q^3 + 63/10*q^2
P2_diff(1/100)
631/1000000
P2_lin=P2-P2_diff
P2_lin
1123/100*q + 61/1000
q0=P2_lin.roots()[0][0]
q0
-61/11230
P2(q0-1/10000)
-1317646813783233867/1416247867000000000000
P2(q0+1/10000)
1843884895441261867/1416247867000000000000
P3=P2(q0+r)
P3
r^3 + 35283/5615*r^2 + 351906913/31528225*r + 32878756/177030983375
P3_diff=termes_grand_degre(P3)
P3_diff
r^3 + 35283/5615*r^2
P3_diff(1/10000)
70567123/1123000000000000
P3_lin=P3-P3_diff
P3_lin
351906913/31528225*r + 32878756/177030983375
r0=P3_lin.roots()[0][0]
r0
-32878756/1975957316495
P3(r0-1/100000000)
-6781410317536312951892697005138134752911223750719899/61719537630659037720669473970719899000000000000000000000000
P3(r0+1/100000000)
6996165200333154319824600523126839922827724190719899/61719537630659037720669473970719899000000000000000000000000
P4=P3(r0+s)
P4
s^3 + 12416232975111/1975957316495*s^2 + 43578799130857674916984057/3904407316610121599085025*s + 13422175326999498219819346928/7714942203832379715083684246339987375
P4_diff=termes_grand_degre(P4)
P4_diff
s^3 + 12416232975111/1975957316495*s^2
P4_diff(1/100000000)
248324659897411463299/395191463299000000000000000000000000
P4_lin=P4-P4_diff
P4_lin
43578799130857674916984057/3904407316610121599085025*s + 13422175326999498219819346928/7714942203832379715083684246339987375
s0=P4_lin.roots()[0][0]
s0
-13422175326999498219819346928/86109846986684169676738889168418120215
//-- Conclusion --//
app=y0+p0+q0+r0+s0
app
180361507581342374686204847776335588181/86109846986684169676738889168418120215
numerical_approx(app)
2.09455148154233
n(P(app))
1.52669567034175e-19
def raphson(g, b=2, c=5): while abs(g^3-b*g-c)>(1/10^15): if 3*g^2-b<>0: x=(c+b*g-g^3)/(3*g^2-b) g=g+x return g
raphson(2,2,5)
180361507581342374686204847776335588181/86109846986684169676738889168418120215
n(raphson(2,2,5))
2.09455148154233
numerical_approx(app) == n(raphson(2,2,5))
True
//-- on obtient la même approximation à la solution --//
P.derivative()
3*y^2 - 2
//-- Remarque: la dérivée de P est équivalente au dénominateur de Δx soit Δx= (c+b*g-g^3)/P.derivative --//
from sage.geometry.newton_polygon import NewtonPolygon NP = NewtonPolygon([ (0,0), (1,1), (2,6) ]) polygon = NP.plot() polygon.show()
Image in a Jupyter notebook