Een rappere Newton-Raphson Edward Omey EHSAL (Stormstraat 2, 1000 Brussel) [
[email protected]]
1. Inleiding Bij vele kwantitatieve problemen is het nodig om nulpunten te bepalen van functies. Soms kunnen deze nulpunten expliciet gevonden worden omdat men (door toeval) de vergelijking f(x) = 0 exact kan oplossen. In andere situaties is het zoeken van de nulpunten niet zo evident en moeten er numerieke procedures gevolgd worden om tot een benaderende oplossing te komen. In vele toestellen is nu een “SOLVER” of een “OPLOSSER” aanwezig. Deze tool is in feite de voordeur van een ganse waaier aan technieken en wiskundige analyse! Tijd om daar even bij stil te staan.
2. Newton-Raphson We vertrekken van een functie f(x) en zoeken een punt (de punten) t waarvoor f(t) = 0. Een populaire aanpak bestaat er in deze nulpunten numeriek te benaderen via een éénstaps iteratief schema van de vorm (1) x( n + 1) = g ( x (n)), n = 0,1,2,... Hierbij is g(x) een geschikte functie waarvoor g(t) = t en waarbij x(0) oordeelkundig moet worden gekozen. Omdat g(t) = t noemt men t een vast punt van de functie g(x). Indien men x(0) = t kiest, volgt uit (1) dat x(n) = t voor alle waarden van n =1, 2, ... Bij de werkwijze van Newton-Raphson wordt de functie g(x) als volgt bepaald: Stap 1. Kies een startpunt x(0) Stap 2. Bepaal f(x(0)) en bepaal de vergelijking van de raaklijn aan f(x) door het koppel (x(0), f(x(0)). De vergelijking van de raaklijn is gelijk aan: y = f ( x (0)) + f ' ( x (0))( x − x (0))
(2)
Stap 3. Zoek het snijpunt van de raaklijn met de horizontale as. Via formule (2) vinden we hier f ( x(0)) x(1) = snijpunt = x(0) − f ' ( x(0)) We noteren dit snijpunt als x(1) . Zie figuur 1. Stap 4. We vervangen x(0) door x(1) en gaan terug naar stap 2. Stap 5. We bepalen achtereenvolgens x(1), x(2), ...
In figuur 1 tekenden we de raaklijn aan f(x) = 15 – x² in x(0) = 2. De raaklijn heeft als vergelijking y = 11 – 4(x – 2). Het snijpunt met de horizontale as is gelijk aan x(1) = 19/4=4,75.
Figuur 1 De bovenstaande procedure leidt tot een rij x(0), x(1), x(2),... waarbij f ( x(n)) x(n + 1) = x(n) − f ' ( x(n)) Dit betekent dat (1) geldt met f ( x) g ( x) = x − f ' ( x)
(3)
Bemerk dat f(t) = 0 enkel en alleen als g(t) = t maar dat er problemen kunnen zijn indien ook f’(t) = 0. In de meeste gevallen kunnen we op met deze werkwijze de nulpunten vinden van afleidbare functies f(x). Deze werkwijze noemt men de werkwijze van Newton-Raphson. Isaac Newton formuleerde dit procédé voor het eerst in 1669. In 1690 vereenvoudigde Joseph Raphson de aanpak van Newton.
Voorbeeld 1. We zoeken de nulpunten van f(x) = x² – 3. We vinden hier f’(x) = 2x en g(x) = x – (x² – 3)/(2x) = x/2 + 3/(2x). De procedure van Newton-Raphson leidt tot de volgende tabel en figuur: n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
x(n) 5,0000000 2,8000000 1,9357143 1,7427649 1,7320837 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508 1,7320508
x(n) -0,5000000 -3,2500000 -2,0865385 -1,7621632 -1,7323081 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508 -1,7320508
6,0000000 5,0000000 4,0000000 3,0000000 2,0000000 1,0000000 0,0000000 -1,0000000
1
3
5
7
9
11
13
15
-2,0000000 -3,0000000 -4,0000000
Wanneer we als startwaarde x(0) = 5 nemen, dan vinden we het eerste nulpunt van f(x); bij de startwaarde x(0) = – 0,5 vinden we de negatieve wortel. Uit de tabel en de grafiek blijkt dat de convergentiesnelheid best meevalt.
Voorbeeld 2. Voor de functie f ( x) = n 0,000 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000 11,000 12,000 13,000 14,000 15,000
x(n) 100,000 170,469 257,863 328,880 355,262 357,679 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696
3620 + ln( x) − 16 , x > 0 vinden we de volgende tabel: x
x(n) 200,000 286,530 342,974 357,059 357,695 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696
x(n) 500,000 296,042 346,626 357,336 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696 357,696
x(n) 1000,000 -1088,643 #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM! #NUM!
Voor de startwaarden x(0) = 100, 200, 500 vinden we de benadering t ≈ 357,696. Bij de startwaarde x(0) = 1000 leidt het procédé niet tot een oplossing. Oefening. 1) Neem andere startwaarden x(0) (zoals bijvoorbeeld x(0) = 5 000; x(0) = 100 000; x(0) = 8 000 000) en ga na of x(n) convergeert of niet. 2) Schets een grafiek van f(x). 3) Verklaar waarom x(0) = 1000 niet werkt. 4) Heeft f(x) nog (een) ander (e) nulpunt(en)? Zo ja, zoek een goede benadering. Oefening. Werk de volgende voorbeelden uit en ga na wat er verkeerd loopt. 1) f(x) = (x – 3)³ (nulpunt voor t = 3, maar f’(3) = f”(3) = 0) 2) f(x) = sin(x) (vele nulpunten en je weet (bijna) niet naar welk nulpunt x(n) convergeert). 3) f(x) = x² + 2 (de rij x(n) oscilleert rond het minimul van f(x); f(x) heeft géén nulpunten). Opmerking. Er bestaan ook twee-staps interatieve (en uiteraard ook drie-staps enzovoort procedures) schema’s. In dit geval gebruikt men in de plaats van (1) een formule van de vorm x( n + 1) = g ( x ( n), x( n − 1)) waarbij de functie g(x,y) en x(0),x(1) op een verstandige manier moeten worden gekozen. In §6 komen een aantal andere werkwijzen aan bod.
3. MacLeod De wiskundige MacLeod paste het procédé van Newton-Raphson aan. In de plaats van de functie g(x) zoals in (3), gebruikte MacLeod de volgende functie: f ( x) g ( x) = x − (4) f ' ( x) + cf ( x) Hierbij is c een reële constante die we oordeelkundig moeten kiezen. De oorsprong van dit voorstel komt verder in § 4 aan bod. Bij de keuze c = 0 vinden we uiteraard formule (3) terug. Voorbeeld 1. (vervolg) We bestuderen f(x) = x² − 3 en gebruiken (4) met c = -1/3. We vinden de volgende tabel en figuur:
n 0,000 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000 11,000 12,000 13,000 14,000 15,000 16,000 17,000
x(n) 4,000 0,455 1,973 1,728 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732 1,732
x(n) -6,000 -4,565 -3,382 -2,501 -1,966 -1,761 -1,733 -1,732 -1,732 -1,732 -1,732 -1,732 -1,732 -1,732 -1,732 -1,732 -1,732 -1,732
c = -1/3 6,000 4,000 2,000 0,000 -2,000
1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16 17
-4,000 -6,000 -8,000
Oefening. Werk voorbeeld 1 uit voor andere keuzes van c en x(0).
4. Convergentiesnelheid Uit de vorige oefening zal blijken dat c een belangrijke rol speelt en dat de convergentie niet steeds even vlot verloopt. Hoe kunnen we nu op een verantwoorde manier het getal c kiezen? Om deze vraag te beantwoorden stellen we e(n) gelijk aan de benaderingsfout die we maken bij de n-de stap: e(n) = x(n) – t. We vinden dan dat x( n + 1) = g ( x ( n)) = g (e( n) + t ) .
Via een reeksontwikkeling(zie NOOT verderop) vinden we nu dat 1 1 g (e(n) + t ) = g (t ) + e(n) g ' (t ) + e²(n) g ' ' (t ) + e³( n) g ' ' ' (t ) + ... 2 6 Omdat f(t) = 0 en dus ook g(t) = t volgt dat 1 1 x(n + 1) = t + e(n) g ' (t ) + e²(n) g ' ' (t ) + e³(n) g ' ' ' (t ) + ... 2 6 1 1 en e(n + 1) = e(n) g ' (t ) + e²(n) g ' ' (t ) + e³(n) g ' ' ' (t ) + ... 2 6 We berekenen nu g’(t) en vinden f '²( x) − f ( x) f ' ' ( x) g ' ( x) = 1 − ( f ' ( x) + cf ( x))² g’(t) = 0. 1 Uit het voorgaande volgt nu dat e(n + 1) ≈ e²(n) g ' ' (t ) . 2 Als g”(t) klein is, dan verwachten we dat de convergentiesnelheid best meevalt. Bij elke stap is de volgende fout evenredig met het kwadraat van de vorige fout. Men noemt deze methode daarom een tweede-orde methode. f ' ' (t ) Na (veel) rekenwerk vinden we eveneens dat g ' ' (t ) = 2c + . Indien we c kunnen f ' (t ) 1 kiezen zo dat g”(t) = 0, dan volgt dat e(n + 1) ≈ e³(n) g ' ' ' (t ) . In dit geval zal de 6 convergentiesnelheid nog beter zal zijn! Bij deze keuze is de methode een derde-orde methode. Het probleem hier is wel dat we t en bijgevolg ook c niet kennen.
5. Veralgemening In Omey (1988) werd de aanpak van Macleod veralgemeend. Hierbij kiezen we een functie h(x) zodanig dat h(t) NIET gelijk is aan nul. Als we nu de functie F(x) = h(x)f(x) bekijken, dan is f(t) = 0 enkel en alleen als F(t) = 0. Mogelijke keuzes voor h(x) zijn: 1 h( x) = exp(cx) = e cx ; h( x ) = exp( cx ²) ; h( x ) = ; h( x) = 1 + f ²( x ) ; enz. 1 + x² Wanneer we de methode van Newton-Raphson toepassen voor F(x) = h(x)f(x), dan vinden we F ( x) f ( x ) h( x ) g ( x) = x − = x− (5) F ' ( x) f ' ( x ) h( x ) + f ( x ) h' ( x ) Wanneer h(x) = 1, dan vinden we formule (3) terug. Wanneer we de keuze h(x) = exp(cx) maken, dan vinden we (4) terug.
Voorbeeld 3 We bestuderen f(x) = x³ – 0,165x² +3,993/10000. Via een grafiek zien we dat de functie drie nulpunten heeft (in de buurt van 0, rond 0,05 en rond 0,15). We zoeken het grootste nulpunt. In de tabel hieronder staan de opeenvolgende benaderingen waaarbij we eerst werkten met (3) (volle lijn) en daarna met (5) (stippellijn) met als keuze h(x) = 1/x.
n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
formule (3) x(n) 1,0000000 6,6852026 4,4754387 3,0024143 2,0206248 1,3664375 0,9308192 0,6411622 0,4491807 0,3228592 0,2411102 0,1902445 0,1615759 0,1491797 0,1464864 0,1463598 0,1463595 0,1463595 0,1463595 0,1463595 0,1463595
formule (5) x(n) 1,0000000 5,0415811 2,5626990 1,3239013 0,7057454 0,3989899 0,2496805 0,1813364 0,1549162 0,1477785 0,1465502 0,1463840 0,1463626 0,1463599 0,1463596 0,1463595 0,1463595 0,1463595 0,1463595 0,1463595 0,1463595
8,0000000 7,0000000 6,0000000 5,0000000 4,0000000 3,0000000 2,0000000 1,0000000 0,0000000 1
3
5
7
9
11
13
15
17
19
21
We vinden als nulpunt bij benadering t ≈ 0,1464595. We merken ook dat (5) vlugger convergeert. Voor formule (5) vinden we, net zoals in de vorige paragraaf, dat g(t) = t g’(t) = 0 h' (t ) f ' ' (t ) g ' ' (t ) = 2 + h(t ) f ' (t ) 1 en e( n + 1) ≈ e ²(n) g ' ' (t ) . 2 Als g”(t) klein is dan verwachten we dat de convergentiesnelheid goed zal zijn. Bij de keuze h( x) = 1 / f ' ( x) vinden we dat g”(t) = 0. In dit geval vinden we dat 1 (6) e( n + 1) ≈ e³( n) g ' ' ' (t ) 6 De enige (belangrijke!) voorwaarde is dat h(x) minstens in de buurt van x = t niet gelijk is aan 0. Voorbeeld 4 De hoeveelheid (in miljoenen) datapaketten q(t) die op een bepaald netwerk worden verstuurd neemt toe met de tijd t. Voor een bepaald netwerk vonden Waner en Costenoble (1996) de volgende formule: 2 exp(0,69t ) q(t ) = ,t≥0 3 + 1,5 exp(−0,4t ) De grafiek vertoont een sterk stijgend patroon. q(t) 180000 160000 140000 120000 100000 80000 60000 40000 20000 18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0
Men is nu geïnteresseerd in het tijdstip waarop q(t) precies gelijk is aan 20 000. Op de grafiek zien we dat er juist één dergelijk tijdstip is. Met de voorgaande methodes zoeken we nu het nulpunt van f(x) = q(x) – 20000 In de volgende tabel en grafiek gebruikten we formule (3) (volle lijn) en formule (5) (stippellijn) met h(x) gelijk aan h(x) = exp(-0,345x). We startten met x(0) = 4.
n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
formule (3) x(n) 4,00000 31,00905 29,55980 28,11058 26,66147 25,21264 23,76458 22,31859 20,87823 19,45305 18,06820 16,78647 15,74304 15,12775 14,95371 14,94239 14,94235 14,94235 14,94235 14,94235 14,94235 14,94235 14,94235 14,94235 14,94235
formule (5) x(n) 4,00000 6,37750 9,03697 11,76066 14,05766 14,91203 14,94047 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051 14,94051
35 30 25 20 15 10 5 0 1
3
5
7
9
11
13
15
17
19
21
23 25
We merken dat de convergentiesnelheid bij formule (5) beduidend beter is dan bij formule (3). Oefening. Het verspreiden van nieuwe technologie kan bijvoorbeeld gemodelleerd worden als volgt (zie [6], p. 295): 3 exp(−0,1t ) p (t ) = ,t ≥ 0 1 + exp(4,50 − 0,477t ) waarbij t de tijd is en p(t) het percentage firma’s is die op tijdstip t een bepaalde technologie toepassen. In de grafiek zien we de opkomst en het verval van de nieuwe techniek.
p(t) 0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1 0 1
4
7 10 13 16 19 22 25 28 31 34 37 40 43
1) Een techniek is goed ingeburgerd indien 50% van de bedrijven de techniek gebruiken. Op welk tijdstip gebeurt dit hier? 2) Een techniek dooft uit als nog hoogstens 10% van de bedrijven de techniek gebruikt. Wanneer gebeurt dit in dit voorbeeld?
6. Andere werkwijzen Naast de formules van de vorige paragrafen zijn er nog tal van varianten. Ook zijn er procedures beschikbaar om de nulpunten te vinden van functies in meerdere variabelen. Hier gaan we niet dieper op in.
6.1. Bissectie In deze procedure gaat men er van uit dat men ongeveer weet waar het nulpunt van f(x) ligt. Men weet bijvoorbeeld dat t in het interval [a, b] ligt en we veronderstellen hier dat f(a) < 0 < f(b). Men bepaalt nu een rij x(n) als volgt:
1. [a(0), b(0)] = [a, b] en x(0) = (a + b)/2 2. Î als f(x(0)) < 0, dan is [a(1), b(1)] = [x(0), b(0)] en x(1) = (x(0) + b(0))/2 Î als f(x(0) > 0, dan is [a(1), b(1)] = [a(0) ,x(0)] en x(1) = (a(0) + x(0))/2 3. ... 4. Î als f(x(n)) < 0, dan is [a(n+1), b(n+1)] = [x(n), b(n)] en x(n+1) = (x(n) + b(n))/2 Î als f(x(n) > 0, dan is [a(n+1), b(n+1)] = [a(n) ,x(n)] en x(n+1) = (a(n) + x(n))/2
6.2. Halley’s methode Dit is de methode beschreven door formule (5) met de keuze h( x) = 1 / f ' ( x) . Bij deze keuze is g(t) = g’(t) = g’’(t) = 0 en voldoen de fouten aan (6).
6.3. Schröder’s methode Dit is de methode beschreven door formule (5) met de keuze h( x ) = 1 / f ' ( x ) .
6.4. Secant methode Bij deze werkwijze vertrekt men van de formules (1) en (3): f ( x(n)) x(n + 1) = x(n) − . f ' ( x(n)) en via de middelwaardestelling benadert men f’(x(n)) als volgt: f ( x(n)) − f ( x(n − 1)) f ' ( x(n)) ≈ x(n) − x(n − 1) Door deze twee formules te combineren vinden we: f ( x(n))( x(n) − x(n − 1)) x(n + 1) = x(n) − ( f ( x(n)) − f ( x(n − 1)) Dit is een twee-staps formule van de vorm x(n+1)=g(x(n),x(n-1)).
Noot De reeksontwikkeling van een functie f(x) rond x = a is gelijk aan 1 1 f ( x ) = f (a ) + f ' ( a )( x − a ) + f ' ' ( a )( x − a )² + f ' ' ' ( a )( x − a )³ + ... 2 3! De functie moet wel aan een aantal vereisten voldoen om over een zincolle uitdrukking te beschikken. Welbekende voorbeelden zijn (telkens met a = 0) 1 = 1 + x + x ² + x ³ + ... 1− x 1 1 e x = 1 + x + x ² + x 3 + ... 2 3!
Bibliografie [1] A..K. Kaw (2006). An interactive e-book for illustrating Newton-Raphson method of solving nonlinear equations. http://numericalmethods.eng.usf.edu/ebooks/newtonraphson_03nle_ebook.htm [2] Macleod, A.J. (1984). How to accelerate convergence in Newton-Raphson? Int. Journal of Math. Education, Science and Technology. Vol. 15, p. 117. [3] Omey, E. (1988). Note on a paper of MacLeod. Int. Journal of Math. Education, Science and Technology. Vol. 19, p. 341-342. [4] T.R. Scavo and J.B. Thoo (1995). On the geometry of Halley’s method. Amer. Math. Monthly 102, 417-426. [5] E. Schröder (1870). Uber unendlich viele Algorithmen zur Auflösung der Gleichungen. Math. Ann. 2, 317-365. [6] Waner, S. en Costenoble, S. (1996). Calculus applied to the real world. Harper Collins College Publisher, New York. [7] E.W. Weinstein. Secant method. Mathworld-A Wolfram Web Resource. http://mathworld.com/SecantMethod.html.
