TOPIC 6: Numerieke Benaderingen¶
Tot nu toe zijn we ervan uitgegaan dat iedere (wiskundige) definitie
rechttoe rechtaan in Python te vertalen valt. Als het gemiddelde van 2
getallen gedefinieerd is als de som van beide getallen gedeeld door 2,
dan is een Python functie schrijven die het gemiddelde uitrekent
kinderspel. Ook functies zoals fac en fib die een inductieve
definitie hebben kunnen rechtstreeks vertaald worden naar Python syntax.
Maar dit is lang niet altijd het geval. Beschouw bijvoorbeeld de
doodgewone vierkantswortel. Een wiskundige zou deze kunnen definiëren
als \(\sqrt{x} = y \iff x = y^2\). Voor een computerwetenschapper is dit
echter niet voldoende. De definitie zegt wel wat een wortel is, maar
zegt niks over hoe men die wortel uitrekent. Dit verschil tussen “wat”
en “hoe” is bepalend voor de computerwetenschapper. Het zet
computerwetenschappen op de kaart als een combinatie van wetenschap en
ingenieurswetenschap. Het komt er niet alleen op aan om dingen te
bestuderen en definiëren maar ook om dingen te maken die technologisch
werken.
De zojuist vernoemde definitie van de vierkantswortel illustreert het verschil tussen twee soorten kennis die voor computerwetenschappers van zeer groot belang zijn:
Declaratieve kennis
is kennis die beschrijft wat iets is. De wiskundige definitie van de wortel is een zeer goed voorbeeld van declaratieve kennis. Declaratieve kennis is heel aantrekkelijk voor mensen aangezien ze dikwijls op heel elegante en abstracte manier iets definieert. Ze is echter dikwijls onbruikbaar voor computers.
Operationele kennis
is kennis die op operationele wijze beschrijft hoe je iets kan aanreiken dat door declaratieve kennis beschreven wordt. Operationele kennis wordt door mensen dikwijls als “technisch” en detaillistisch ervaren. Maar het is dit soort kennis die computers nodig hebben om een vraagstuk tot een goed einde te brengen.
Het omzetten van declaratieve kennis naar operationele kennis is dikwijls het kernprobleem dat bij het schrijven van een programma moet worden opgelost. De reden waarom deze omzetting niet triviaal is, is dat we van een operationele definitie verwachten dat ze een oplossing genereert binnen een aanvaardbaar tijdsbestek.
Veel kennis die opgebouwd werd in de wiskunde is declaratieve kennis. In dit hoofdstuk bestuderen we enkele algoritmen die de operationele variant vormen van enkele zeer gekende wiskundige vraagstukken. We bespreken achtereenvolgens algoritmen om bepaalde integralen, vierkantswortels en nulpunten van functies uit te rekenen. Deze algoritmen leveren geen precies antwoord maar trachten het theoretische antwoord zo goed mogelijk te benaderen. Zo’n benaderingsmethode of approximatiemethode zal typisch in opeenvolgende stappen potentiële oplossingen aanreiken die steeds “nauwer bij” of “korter bij” de theoretische oplossing liggen.
Numeriek Uitrekenen van Integralen¶
Integreren ken je van de cursus wiskunde. Ter herinnering vermelden we dat \(F(x) = \int f(x)dx\) de onbepaalde integraal heet van de functie \(f(x)\). \(F(x)\) wordt dan ook wel de stamfunctie van \(f(x)\) genoemd. De onbepaalde integraal wordt ook gebruikt om een bepaalde integraal uit te rekenen over een interval \([a,b]\). We hebben dan \(\int_a^b f(x)dx = F(b)- F(a)\).
Integreren heel dikwijls wat meer voeten in aarde dan afleiden. Er zijn immers geen eenvoudige rechttoe-rechtaan regels om de stamfunctie van een functie te bepalen. Het uitrekenen van stamfuncties is dan ook een vaardigheid, een kunst, waarin we onszelf wel kunnen trainen (en er dus steeds bedrevener in worden) maar die we nooit kunnen opschrijven d.m.v. een stelselmatige procedure die altijd werkt. Sommige geavanceerde rekenmachines gaan redelijk ver in het zoeken van stamfuncties maar het zal altijd mogelijk zijn om één of andere functie op te schrijven waar gelijk welke rekenmachine niet uit raakt.
Voor het berekenen van onbepaalde integralen \(\int f(x)dx\) is dat effectief een probleem. Voor het berekenen van bepaalde integralen \(\int_{a}^{b}f(x)dx\) is dat eigenlijk niet zo’n probleem omdat we immers niet écht geïnteresseerd zijn in de stamfunctie \(F\) maar wel in het getal dat we uiteindelijk verkrijgen na het uitrekenen van \(F(b)-F(a)\). Daarom hebben numerieke analisten een aantal numerieke integratiemethoden ontwikkeld die integralen kunnen uitrekenen zonder daarbij via een stamfunctie te gaan. In deze sectie bespreken we twee zulke methoden: de rechthoeksregel en de trapeziumregel. Deze twee methoden vormen slechts een klein onderdeeltje van het arsenaal aan numerieke integratietechnieken die bekend zijn in de numerieke analyse.
De rechthoeksregel¶
Eén van de meest eenvoudige regels om integralen numeriek uit te rekenen is de rechthoeksregel. Deze regel wordt geïllustreerd inonderstaande figuur. Hij correspondeert hoogstwaarschijnlijk met de intuïtie die je van integralen hebt. We benaderen de integraal door het oppervlak bepaald door de functie op te delen in rechthoeken. De reden waarom we de integratieoppervlakte opdelen in rechthoeken is dat de oppervlakte van iedere rechthoek zeer eenvoudig te berekenen is. Dat kunnen we immers al van op de lagere school. De breedte van iedere rechthoek wordt bepaald door de gewenste precisie: hoe meer rechthoekjes we kiezen, hoe minder breed iedere rechthoek is. De breedte noteren we met \(\Delta\). Het aantal rechthoekjes noteren we met \(n\). We hebben dus \(\Delta = \frac{b-a}{n}\). De hoogte van iedere rechthoek wordt bepaald door de functiewaarde te berekenen in het middelpunt van ieder deelinterval.
Met wiskundige notatie wordt de rechthoeksregel als volgt opgeschreven:
Hoe komen we hieraan? De breedte van iedere rechthoek is \(\Delta\) dus moeten we zeker iedere keer vermenigvuldigen met \(\Delta\). De hoogte van iedere rechthoek correspondeert met de functiewaarde genomen in het midden van de rechthoek. Het eerste midden ligt op \(a+\frac{\Delta}{2}\), het tweede midden op \(a+2.\frac{\Delta}{2}\) en het \(k\)de midden dus op \(a+k.\frac{\Delta}{2}\). De totale integraal is de som van alle oppervlaktes van al die rechthoekjes. Het volgende stukje code laat zien hoe dat er in Python uitziet:
def rect(f,a,b,n):
delta = (b - a)/n
start = a + delta/2
area = 0
for k in range(n):
area = area + f(start+k*delta)*delta
return area
Om de code uit te testen zullen we de sin functie als voorbeeld gebruiken. pi en sin zitten in de math bibliotheek. We ‘weten’ uit de wiskundeles dat bv. op het interval [0,pi/2]de integraal 1 is. De resultaten in onze experimenten benaderen 1. Een grotere n geeft een betere benadering. Ook het experiment op het interval [-pi,+pi] geeft zoals verwacht een extreem klein getal terug, i.e. een benadering van 0 (de teruggegeven waarde is een heel klein getal in de wetenschappelijke notatie voor heel kleine floats; 1.2e-16 staat voor \(1.2 \times 10^{-16}\)).
from math import sin, pi
rect(sin,0,pi/2,10)
1.0010288241427083
rect(sin,0,pi/2,100)
1.0000102809119051
rect(sin,-pi,+pi,100)
1.196959198423997e-16
De functie rect neemt vier parameters: de functie f die we willen
integreren, \(a\) en \(b\), en het aantal rechthoekjes dat we in beschouwing
willen nemen. Hoe groter \(n\), hoe preciezer het resultaat zal zijn.
delta is de breedte van iedere rechthoek en het eerste punt waarop we
de functiewaarde moeten uitrekenen is a+delta/2 aangezien dat precies
het midden is van de eerste rechthoek. Vanaf dan laten we k variëren
over het totaal aantal rechthoekjes en tellen we telkens de oppervlakte
van ieder rechthoekje bij de variabele area dewelke de totale
oppervlakte van de integraal vergaart. Op het eind geven we deze
variabele terug als resultaat van het integratieprocedé.
De trapeziumregel¶
De trapeziumregel is een verbetering van de rechthoeksregel. Hij wordt geïllustreerd in onderstaande figuur. Intuïtief voelen we aan dat het stukje onder de functie benaderen met een trapezium een preciezer resultaat zal opleveren dan de rechthoeksregel. De schuine zijde van het trapezium volgt immers veel beter de helling van de functie indien die erg grillig is. We verwijzen door naar cursussen wiskunde om te bepalen hoeveel preciezer deze methode juist is.
De volgende Python functie naive_trap is een eerste versie van onze
trapeziumregel. Zoals je aan de naam naive_trap van de functie kan
zien volgt er zo meteen een verbeterde oplossing. De structuur van de
code is heel gelijklopend met die van de rechthoeksregel. Het enige wat
verschilt is de oppervlakte die we bij area tellen in iedere slag van
de iteratie. In de trapeziumregel is dat de oppervlakte van een
trapezium, wat gegeven wordt als de breedte maal het gemiddelde van
beide basissen.
def naive_trap(f,a,b,n):
delta = (b - a)/n
area = 0
for k in range(n):
area = area + (f(a+k*delta)+f(a+(k+1)*delta))/2*delta
return area
Merk vooreerst nog eens op dat range(n) alle mogelijke gehele getallen
tussen \(0\) en \(n-1\) bevat. De eerste keer is \(k\) dus gelijk aan \(0\). De
oppervlakte van het eerste trapezium is dus
\(\Delta \frac{f(a)+f(a+\Delta)}{2}\), die van het tweede trapezium is
\(\Delta \frac{f(a+\Delta)+f(a+2 \Delta)}{2}\) enzovoort. In het algemeen
is dat dus \(\Delta \frac{f(a+k \Delta)+f(a+(k+1)\Delta)}{2}\) wat meteen
ook de body van de functie naive_trap verklaart. We herhalen dezelfde experimenten met sin om onze code te testen.
naive_trap(sin,0,pi/2,100)
0.9999794382396072
naive_trap(sin,0,pi/2,10)
0.9979429863543573
naive_trap(sin,-pi,pi,100)
-1.8041124150158794e-16
Trapeziumregel Versie 2¶
Vanuit computationeel standpunt is bovenstaande methode nogal verspillend. Voor elke \(k\) met \(0<k<n\) wordt \(f(a+k \Delta)\) immers 2 keer uitgerekend (in opeenvolgende slagen van de iteratie) en in iedere slag wordt echter slechts \(\frac{f(a+k \Delta)}{2}\) bij de totale oppervlakte geteld; de helft dus. We kunnen dus eigenlijk minder rekenwerk doen door ervoor te zorgen dat we slechts \(1\) keer \(f(a+k \Delta)\) bij de oppervlakte tellen maar dat getal dan ook nooit delen door \(2\). We moeten echter wel opletten voor de allereerste waarde \(f(a+0 . \Delta)\) en \(f(a+(n+1) \Delta)\). Die worden er in het originele algoritme slechts \(1\) keer voor de helft bijgeteld. Een tweede oorzaak van verspilling is dat de functiewaarden in iedere slag van de iteratie vermenigvuldigd worden met \(\Delta\). Hier kunnen we ook veel rekenwerk uitsparen wegens de distributiviteit van de vermenigvuldiging ten opzichte van de optelling: \(\sum ( f(...) \Delta) = (\sum f(...) ) \Delta\). Het volstaat dus van de vermenigvuldiging met \(\Delta\) helemaal op het einde te doen.
Deze twee computationele besparingen worden in de functie trap
hieronder doorgevoerd.
def trap(f,a,b,n):
delta = (b - a)/n
area = 0
for k in range(1,n):
area = area + f(a+k*delta)
return (area + (f(a) + f(b))/2)*delta
Het algoritme telt nu (vanaf \(k=1\)) iedere functiewaarde \(f(a+k \Delta)\) één keer volledig (dus niet delen door \(2\)) bij de oppervlakte. We moeten dan wel opletten van helemaal op het einde nog eens de helft van de allereerste en de allerlaatste functiewaarde bij het resultaat bij te tellen. De vermenigvuldiging met \(\Delta\) doen we slechts \(1\) keer helemaal op het einde.
We kunnen zeggen dat deze versie van het algoritme grosso modo de helft
minder werk zal verrichten wat het uitrekenen van \(f\) betreft. Bovendien
worden \(n\) vermenigvuldigingen met \(\Delta\) vervangen door slechts één
vermenigvuldiging met \(\Delta\). De verbeterde functie zal dus een heel
pak sneller werken dan de naïeve versie. Voor kleine \(n\) is dit verschil
niet voelbaar maar indien we \(n\) moeten opdrijven om zeer grote
precisies te bereiken kan het verschil enorm groot worden. Ook hier herhalen we de experimenten met sin om de code te testen. De resultaten wijken licht af van deze van de naive-trap fubnctie omdat door de verschillen in bewerkingen ook andree afrondingsfouten gebeuren.
trap(sin,0,pi/2,10)
0.9979429863543572
trap(sin,0,pi/2,100)
0.9999794382396076
trap(sin,-pi,+pi,100)
-6.801343571116833e-17
Merk ten slotte op dat de laatste versie van de trapeziumregel correspondeert met de formule die men typisch in wiskundeboeken zal terugvinden:
Snelheid van de integratiealgoritmen¶
Merk ten slotte op dat het vergroten van \(n\) de enige manier is om onze drie integratiemethoden te verbeteren: hoe groter we \(n\) kiezen, hoe kleiner \(\Delta\) zal zijn en hoe nauwkeuriger het resultaat dus zal zijn. De prijs die we daarvoor betalen zijn meer iteraties en dus een langere looptijd van het programma. Merk ook op dat de iteratie in alle drie de gevallen \(n\) (of \(n-1\)) keer wordt uitgevoerd. De drie gevallen genereren dus een lineair computationeel proces. Het zijn m.a.w. algoritmen die \(O(n)\) zijn.
Wortels Berekenen¶
Het tweede probleem dat we bestuderen bestaat uit het uitrekenen van vierkantswortels. Declaratief weten we dat \(\sqrt{x} = y \iff x = y^2\) maar dat levert ons zoals gezegd geen methode op om de wortel effectief ook uit te rekenen.
Wortels van Perfecte Kwadraten¶
Een eerste oplossing voor het probleem bestaat uit het schrijven van een Python functie die wortels van natuurlijke getallen kan berekenen zodat het resultaat eveneens een natuurlijk getal is. De bedoeling is dus om wortels binnen \(\mathbb{N}\) te bestuderen. Zulke wortels bestaan natuurlijk enkel voor natuurlijke getallen die perfecte kwadraten \(1\), \(2\), \(4\), \(16\), …zijn. Onze Python functie kan er als volgt uitzien:
def sqrt_int(n):
for i in range(n):
if i*i == n: return i
print (n, " is not a perfect square")
return None
sqrt_int(16)
4
sqrt_int(10)
10 is not a perfect square
De functie sqrt_int krijgt een natuurlijk getal \(n\) als parameter en
geeft de wortel ervan terug als die bestaat. Indien de wortel niet
bestaat (bvb. als we de functie zouden uitproberen op \(17\) of \(23\)) dan
wordt dit meegedeeld en wordt er None als waarde teruggegeven.
De techniek die we hier gebruikt hebben om het verschil tussen
declaratieve en operationele kennis te overbruggen draagt de naam
brute-kracht methode of generate-and-test methode. Het idee van deze
methode bestaat erin om voldoende dikwijls te “raden” naar een bepaalde
oplossing. Dat raden gebeurt natuurlijk niet zomaar lukraak maar bestaat
erin van systematisch alle mogelijke oplossingen te proberen. M.a.w. we
genereren telkens een mogelijke oplossing en testen meteen of de
gegenereerde oplossing ook daadwerkelijk een oplossing is. Indien dat
niet het geval is doen we verder door de volgende mogelijke oplossing te
genereren. Deze oplossingsmethode is een brute-kracht oplossing omdat de
functie domweg alle mogelijke oplossingen beschouwt en erop vertrouwt
dat de computer snel en krachtig genoeg is om alle tussenliggende
non-oplossingen te doorlopen en te verwerpen alvorens bij de echte
oplossing uit te komen. De for lus van Python is hiervoor extreem
geschikt aangezien deze ons precies toelaat om te zeggen dat we alle
mogelijke oplossingen aflopen.
Toevoegen van Gezond Verstand¶
Eigenlijk is de laatste zin van de vorige paragraaf niet helemaal
correct. We lopen immers niet echt alle getallen van \(\mathbb{N}\) af
maar beperken ons tot de getallen uit range(n). Dat is omdat we weten
dat een wortel van een getal steeds kleiner is of gelijk aan dat getal.
Het heeft geen zin om grotere getallen te genereren aangezien de test
voor deze getallen nooit zal slagen.
We kunnen echter nog iets slimmer zijn en stoppen van zodra we detecteren dat de laatst gegenereerde waarde in het kwadraat reeds groter is dan \(n\). Ook dan is het zinloos om nog volgende waarden te beschouwen aangezien hun kwadraat telkens groter dan \(n\) zal zijn. Deze kennis is opgenomen in de volgende versie van onze Python functie:
def smart_sqrt_int(n):
for i in range(n):
if i*i == n:
return i
elif i*i > n:
print (n, " is not a perfect square")
return None
smart_sqrt_int(16)
4
smart_sqrt_int(10)
10 is not a perfect square
Deze functie genereert dus waarden vanaf \(0\) tot en met \(n-1\). Indien de
gegenereerde waarde in het kwadraat gelijk is aan onze invoer (t.t.z. de
test slaagt) zijn we klaar en hebben we de wortel gevonden. Indien de
genereerde waarde in het kwadraat groter is aan de invoer zullen we
zeker geen wortel meer vinden aangezien het kwadraat blijft groter
worden. Dan kunnen we dus meteen stoppen. Zo’n stukje extra domeinkennis
(in dit geval dus wiskundekennis) om een brute krachtmethode wat te
versnellen wordt doorgaans een heuristiek genoemd. Merk op dat
Python’s for lus ook hiervoor nog steeds goed werkt: we lopen gewoon
alle mogelijke oplossingen af en de heuristiek is opgenomen als een if
test die de lus vroegtijd kan verlaten.
Wortels in float¶
Laat ons nu overstappen van numerieke vraagstukken m.b.t. het type int
naar vraagstukken in de wereld van floats. We zullen gewoon verder
doen met het uitrekenen van wortels maar dan nu voor reële getallen
i.p.v. natuurlijke getallen. Twee problemen zorgen ervoor dat de brute
krachtmethode niet meer toepasbaar is:
De brute krachtmethode is enkel toepasbaar als er voor iedere gegenereerde oplossing duidelijk een “volgende” potentiële oplossing gegenereerd kan worden. Voor het vraagstuk dat bestaat uit het trekken van natuurlijke wortels is dit inderdaad het geval: we dienen gewoon de natuurlijke getallen één voor één af te lopen om te testen of hun kwadraat al dan niet de invoerwaarde is. Reële getallen zijn echter overaftelbaar waardoor het niet langer mogelijk is om alle getallen op te sommen die potentieel de wortel van een gegeven getal zijn en die telkens te testen. In een computer zijn de reële getallen benaderd door het type
floatdat wel aftelbaar is aangezien een computer een eindig groot geheugen heeft. Alle opeenvolgende floats opsommen is in de praktijk echter onbegonnen werk omdat het onredelijk veel tijd in beslag zou nemen.Een tweede probleem bestaat uit het feit dat de oplossing voor een gegeven vraagstuk dikwijls niet eens bestaat in de wereld van de
floats. In theorie heeft elk positief reëel getal een wortel maar die wortel is niet noodzakelijk een getal dat met een eindig aantal decimalen kan uitgedrukt worden. Bijgevolg is de exacte wortel in een computer niet altijd als eenfloatuit te drukken. We denken bijvoorbeeld aan \(\sqrt{2}\).
We hebben dus een andere methode nodig om oplossingen te genereren. Dit is wat zogenoemde benaderingsmethodes oftewel approximatiemethodes doen. Om het tweede probleem op te lossen gaan we uit van een benadering (oftewel: approximatie) van de echte oplossing en stoppen we de iteratie van zodra de benadering “dicht genoeg” bij de theoretische oplossing ligt. We beschouwen een benadering als “dicht genoeg bij de theoretische oplossing” indien het verschil tussen die benadering en de echte oplossing klein genoeg is, t.t.z. kleiner dan een vooraf gekozen precisie. Zelfs indien een gegenereerde oplossing nog niet de exacte oplossing is stoppen we het procedé indien onze gekozen precisie dit toelaat. Om het eerste probleem aan te pakken gaan we op een intelligente wijze steeds een nieuwe benadering genereren in de hoop dat dit procedé snel genoeg tot een oplossing leidt die voldoende dicht bij de theoretische oplossing ligt.
Iteratie met while¶
In de vorige paragraaf treedt een nieuw gegeven op waar we tot nu toe
nog geen Python oplossing voor hebben: we moeten iets kunnen herhalen
tot een bepaalde voorwaarde voldaan is (nl. voldoende precies). Met
for is dit niet mogelijk aangezien we moeten aangeven hoe dikwijls de
body herhaald dient te worden. Dat is waar het while statement uit
Python voor dient. In zijn algemene vorm ziet het er als volgt uit:
Dit statement zal beginnen met de conditie te testen. Indien de test
True oplevert zullen de statements van het block één voor één worden
uitgevoerd. Nadien wordt de conditie opnieuw getest. Enzovoort. Dat
blijft dus doorgaan tot de conditie False oplevert. Op dat ogenblik
stopt de iteratie. Merk op dat het niet gegarandeerd is dat het block
überhaupt wordt uitgevoerd. Indien de conditie van meet af aan False
is gebeurt er helemaal niks.
Veronderstel dat we willen bestuderen hoe snel bacteriën kweken indien de populatie op tijdstip \(t+1\) gegeven is in functie van de populatie op tijdstip \(t\) met de volgende vergelijking:
waarbij \(r\) een of andere groeifactor is die model staat voor de omgevingsfactoren (bvb. de hoeveelheid voedsel). In dit vereenvoudigde model gaan we er dus van uit dat er geen bacteriën afsterven. We kunnen ons nu de vraag stellen na hoeveel tijdstippen het aantal bacterieën in de kolonie verdubbeld zal zijn. Met andere woorden: hoe dikwijls dienen we bovenstaande formule uit te rekenen alvorens we het dubbele krijgen van \(P(0)\). Dit is eenvoudig te bepalen door volgend Python programma:
def grow(init_pop, r):
time = 0
pop = init_pop
while pop < init_pop * 2:
print (pop)
pop = pop + r * pop
time = time + 1
return time
grow(10,0.2)
10
12.0
14.4
17.28
4
Alvorens we de iteratie in gang steken stellen we de tijd op nul en
initialiseren we de “huidige populatie” pop met de waarde van de
initiële populatie. Dan begint de while loop te werken. In iedere slag
van de iteratie herberekenen we de waarde van pop op basis van de
vorige waarde van pop (volgens bovenstaande formule), verhogen we de
tijd en drukken we de huidige populatie af op het scherm. Indien de
huidige populatie pop minstens het dubbele is van de initiële
populatie stopt de iteratie en wordt time teruggegeven met return.
Een fenomeen dat we bij de iteratie met for niet zijn tegen gekomen is
dat een iteratie kan ontsporen in een oneindig computationeel proces.
Dat is omdat we bij de for bij het begin van de iteratie weten wat de
eerste en laatste waarde van de variabele in de for zal zijn. Bij de
while is dit niet langer waar. Stel dat we bovenstaande functie
oproepen met een negatieve groeifactor. Dan zal de conditie van de
while nooit True worden en zal de iteratie oneindig lang blijven
doorlopen (tot we de iteratie met de
hand stoppen; in deze notebook door de kernel te ‘interupten’). Het is dus opletten geblazen bij het gebruik van
while!
Iteraties die met for worden genoteerd noemt men daarom bounded
loops. Iteraties die met while worden genoteerd noemt men free
loops.
In sectie over lijsten hebben we gezien dat het for statement
gebruikt kan worden om bepaalde recursieve processen eenvoudiger uit te
drukken. Hetzelfde is waar voor while. We kunnen dus ook alles wat met
een while valt uit te drukken recursief schrijven. Hier zien we de
recursieve versie van bovenstaande functie grow:
def grow_aux(init_pop,pop,r,time):
if pop >= init_pop * 2:
return time
else:
print (pop)
return grow_aux(init_pop, pop+pop*r, r, time+1)
def grow2(init_pop,r):
return grow_aux(init_pop, init_pop, r, 0)
grow(10,0.2)
10
12.0
14.4
17.28
4
De Methode van Newton¶
We pikken nu terug de draad op van wortels berekenen. De methode van Newton laat ons toe om via opeenvolgende benaderingen te komen tot een approximatie voor \(\sqrt{x}\) voor een reëel getal \(x\). De methode gaat uit van een gok \(g\) (d.i. de huidige benadering) en zal op basis van die gok en van \(x\) telkens weer een betere gok genereren. Dit blijven we dan herhalen tot we merken dat het verschil tussen \(g^2\) en \(x\) voldoende klein is. Dat laatste doen we door te testen dat \(| (x-g^2) | < \epsilon\) waarbij \(\epsilon\) de gekozen precisie is.
De inbreng van Newton is dat indien \(g\) een gok is voor \(\sqrt{x}\), dan is \(\frac{g+\frac{x}{g}}{2}\) een betere gok. De bedoeling is dus om deze breuk te blijven herberekenen tot \(g\) voldoende dicht bij de wortel van \(x\) zit. In Python doen we dit als volgt:
def sqrt(x):
guess = 1.0
prec = 0.001
while abs(x-(guess*guess)) > prec:
guess = (guess + x/guess)/2
return guess
sqrt(16)
4.000000636692939
sqrt(10)
3.162277665175675
prec is de precisie die hierboven als \(\epsilon\) genoteerd werd.
Hieronder laten we zien wat er gebeurt als we sqrt op
verschillende waarden loslaten. Zoals we kunnen zien is de
uitkomst niet helemaal precies: zo is bijvoorbeeld de wortel van
\(16\) niet exact gelijk aan \(4\). Dit is eigen aan numerieke
benaderingsmethoden. De methode stopt
immers met verder rekenen van zodra een oplossing gevonden
is die nauw genoeg (volgens de geprogrammeerde precisie) aansluit bij de echte oplossing.
Snelheid van dit algoritme¶
Hoe snel loopt dit algoritme? Laten we enkele uitvoeringen ervan tracen d.m.v. een aangepaste versie van het algoritme die we traced_sqrt hebben genoemd. We hebben de functie traced_sqrt twee extra parameters gegeven zodat we de startwaarde en de precisie als argument kunnen meegeven in plaats van ze als constanten te beschouwen. Er zijn verder goed geformateerde print instructies aan de code toegevoegd zodat we de waarde van de gok en van het kwadraat van de gok kunnen opvolgen.
def traced_sqrt(x, prec, guess):
print("current guess:", "%16.8f"% guess, " guess*guess:", "%16.8f"% (guess*guess))
while abs(x-(guess*guess)) > prec:
guess = (guess + x/guess)/2
print("current guess:", "%16.8f"% guess, " guess*guess:", "%16.8f"% (guess*guess))
return guess
Hieronder zien we eerst het resultaat van de functie voor \(x=25\) met precisie 0.0001 en startwaarde 1.
traced_sqrt(25,0.0001,1)
current guess: 1.00000000 guess*guess: 1.00000000
current guess: 13.00000000 guess*guess: 169.00000000
current guess: 7.46153846 guess*guess: 55.67455621
current guess: 5.40602696 guess*guess: 29.22512752
current guess: 5.01524760 guess*guess: 25.15270851
current guess: 5.00002318 guess*guess: 25.00023178
current guess: 5.00000000 guess*guess: 25.00000000
5.000000000053722
Laten we nu eens bij wijze van experiment eens kijken of het kiezen van een andere precisie een grote invloed heeft. We zijn nu al tevreden met een precisie op 2 decimalen. We zijn net 1 iteratie vroeger klaar wat wel strookt met de intuitie dat je langer moet doorgaan als je een grotere preciesie wil.
traced_sqrt(25,0.01,1)
current guess: 1.00000000 guess*guess: 1.00000000
current guess: 13.00000000 guess*guess: 169.00000000
current guess: 7.46153846 guess*guess: 55.67455621
current guess: 5.40602696 guess*guess: 29.22512752
current guess: 5.01524760 guess*guess: 25.15270851
current guess: 5.00002318 guess*guess: 25.00023178
5.000023178253949
Laten we nu eens bij wijze van experiment een veel grotere \(x\) kiezen, bvb. \(x=2500\). Er zijn duidelijk meer iteraties nodig, maar dat is niet noodzakelijk een gevolg van de grotere \(x\). De startgok ligt hier ook een stuk verder van de beoogde uitkomst \(\sqrt{x}\).
traced_sqrt(2500,0.0001,1)
current guess: 1.00000000 guess*guess: 1.00000000
current guess: 1250.50000000 guess*guess: 1563750.25000000
current guess: 626.24960016 guess*guess: 392188.56170048
current guess: 315.12080934 guess*guess: 99301.12447808
current guess: 161.52713731 guess*guess: 26091.01608736
current guess: 88.50220639 guess*guess: 7832.64053583
current guess: 58.37504486 guess*guess: 3407.64586203
current guess: 50.60078221 guess*guess: 2560.43915978
current guess: 50.00356654 guess*guess: 2500.35666655
current guess: 50.00000013 guess*guess: 2500.00001272
50.000000127192884
Als laatste experiment gebruiken we daarom een veel grotere startwaarde voor guess bij ons oorspronkelijk experiment met \(x=25\). Ook hier zien we dat we een aantal iteraties meer nodig hebben, niet omdat \(x\) groter is maar omdat de startgok verder van de beoogde uitkomst \(\sqrt{x}\) ligt.
traced_sqrt(25,0.0001,1000)
current guess: 1000.00000000 guess*guess: 1000000.00000000
current guess: 500.01250000 guess*guess: 250012.50015625
current guess: 250.03124938 guess*guess: 62515.62566403
current guess: 125.06561844 guess*guess: 15641.40891538
current guess: 62.63275675 guess*guess: 3922.86221836
current guess: 31.51595445 guess*guess: 993.25538520
current guess: 16.15460174 guess*guess: 260.97115730
current guess: 8.85107420 guess*guess: 78.34151449
current guess: 5.83779506 guess*guess: 34.07985117
current guess: 5.06011692 guess*guess: 25.60478328
current guess: 5.00035711 guess*guess: 25.00357124
current guess: 5.00000001 guess*guess: 25.00000013
5.0000000127519
Deze traces doen ons vermoeden dat het aantal stappen dat het algoritme uitvoert niet direct afhankelijk is van de grootte van de invoerwaarde \(x\). De absolute fout (t.t.z. het verschil tussen de benadering en de gezochte wortelwaarde) verkleint echter wel in elke iteratie. De convergentie hangt dus niet af van de grootte van de invoerwaarde maar enkel van de gewenste precisie en van de afstand tussen de eerste gok en de werkelijke wortelwaarde. Men kan bewijzen dat de absolute fout als volgt kleiner wordt:
Eens de fout kleiner dan \(1\) is convergeert de methode dus bliksemsnel. Men zegt dat de methode kwadratisch convergent is. Dit wil zeggen dat het aantal significante cijfers verdubbelt in elke slag van de iteratie! Dat betekent dat de fout gedeeld wordt door \(100\) in elke slag van de iteratie. Eens de gok dicht genoeg ligt bij de wortel (en het is niet zo eenvoudig om te bepalen wanneer dit gebeurt) ligt het proces dus in \(O(log_{10^2})\) wat zeer goed is.
Nulpunten van functies¶
Een derde voorbeeld van een numerieke approximatiemethode bestaat erin van nulpunten van willekeurige functies te bepalen. We beschouwen een functie \(f\) van \(\mathbb{R}\) naar \(\mathbb{R}\) die “voldoende glad” is (t.t.z. een voldoend aantal keren continu-differentieerbaar; de details hiervan laten we over aan de cursus wiskunde). Indien we een interval \([a,b]\) beschouwen zodat \(f(a)\) en \(f(b)\) een verschillend teken hebben, dan weten we dat \(f\) minstens één keer de \(x\)-as moet snijden. M.a.w. er moet ergens een \(c\) bestaan tussen \(a\) en \(b\) zodat \(f(c)=0\). Zo’n \(c\) noemen we een nulpunt van \(f\). Indien \(f(x)\) een eenvoudige veelterm is, is het zoeken van een nulpunt niet zo moeilijk. Bijvoorbeeld, voor een veelterm van orde \(2\) kennen we de aloude methode gebaseerd op de discriminant. Voor een willekeurige functie \(f\) is er echter geen eenvoudige analytische oplossing. Dan gaan we over tot een numerieke approximatiemethode.
Opnieuw is de generate-and-test methode hier niet praktisch toepasbaar.
We zouden van \(a\) tot \(b\) alle mogelijke tussenliggende \(c\)-waarden
kunnen aflopen en telkens testen of inderdaad \(f(c)=0\). Net zoals bij
ons wortelprobleem is dit theoretisch onzin omdat er oneindig veel
getallen tussen \(a\) en \(b\) liggen. Ook praktisch is deze methode
onbruikbaar aangezien er enorm veel eindige floats tussen \(a\) en \(b\)
liggen zodat het aantal stappen uitgevoerd door deze methode
onaanvaardbaar groot zou zijn. Bovendien is het nulpunt misschien niet
eens voorstelbaar met floats, bijvoorbeeld omdat het irrationaal is of
een oneindig repetitieve decimale ontwikkeling heeft. We gaan dus
opnieuw op zoek naar een geschikte approximatiemethode die de
opeenvolgende suggesties voor \(c\) min of meer intelligent berekent.
Een zeer elegante methode om dit probleem op te lossen is de zogenoemde bissectiemethode of ook wel de half-interval methode genoemd. Deze wordt grafisch voorgesteld in bovenstaande figuur. Het idee van deze methode is van telkens opnieuw het interval \([a,b]\) te verkleinen om zo dichter en dichter bij de oplossing te komen. Indien het interval klein genoeg is zal zijn middelpunt voldoende dicht bij het nulpunt liggen. We gaan dus telkens het midden van dat interval berekenen (t.t.z. \(\frac{a+b}{2}\)) en dit als voorlopige waarde voor \(c\) kiezen. Indien \(f(c)\) voldoende dicht bij de oplossing zit (t.t.z. \(|f(c)|<\epsilon\)), dan stoppen we. Anders zien we of \(f\) de \(x\)-as snijdt in het interval \([a,c]\) dan wel \([c,b]\) en herhalen we de hele bedoening met dit nieuwe interval. In het meest algemene geval kan de bissectiemethode toegepast worden van zodra \(f(a)\) en \(f(b)\) een verschillend teken hebben. Wij zullen er voor de eenvoud van de code van uitgaan dat \(f(a)<0\) en \(f(b)>0\). We laten het aanpassen van de code tot het algemene geval over als oefening. In Python ziet dat er dan als volgt uit:
def zero(f,a,b):
mid = (a+b)/2
val = f(mid)
while abs(val) > 0.001:
if val < 0:
a = mid
else:
b = mid
mid = (a+b)/2.0
val = f(mid)
return mid
Alvorens we in gaan op de werking van zero, bekijken we hieronder naar
de manier waarop we deze functie kunnen gebruiken. We gebruiken de wiskundige functie
\(test : x \rightarrow x^2-9\) als voorbeeld.
Om de gedachten te vestigen tonen we snel een plot van de voorbeeldfunctie gemaakt met de gekende mathplotlib bibliotheek. We hebben de code om de plot te maken weggestopt in een module plot. Je ziet dat zowel -3 als 3 nulpunten zijn voor deze functie wat in onderstaande plot wordt benadrukt door de rode cirkeltjes.
from plot import plot_voorbeeld
plot_voorbeeld()
We definiëren allereerst de corresponderende Python functie.
Vervolgens roepen we zero op met drie
parameters: de functie test, \(a=-2\) en \(b=4\). Op dat interval verwachten we het nulpunt 3 te vinden.
In het tweede experiment gebruiken we \(a=0\). Ook op dat interval verwachten we het nulpunt 3 te vinden. Zoals we zien is de uitkomst in beide experimenten niet precies dezelfde. Dat komt omdat we in het tweede geval het extreme geluk
hebben dat we na \(2\) iteraties precies bij het nulpunt terecht komen (wat je eventjes uit je hoofd kan controleren maar wat wat verderop ook met de traced versie wordt geïllustreerd).
def test(x): return x**2-9
zero(test,-2,4)
2.9998779296875
zero(test,0,4)
3.0
De werking van zero is niet zo moeilijk te begrijpen. De functie
gebruikt \(2\) lokale variabelen mid (die overeenkomt met de \(c\) uit
onze uitleg) en val (die overeenkomt met \(f(c)\)). Zo lang de afstand
tussen \(f(c)\) en \(0\) (t.t.z. abs(val)) groter is dan de gekozen
precisie, blijft de iteratie doorlopen. We bekijken nu de body van de
iteratie. In iedere slag van de iteratie wordt getest of f(mid)<0.
Indien dat het geval is moeten we verder zoeken tussen \(mid\) en \(b\). We
stellen dus a=mid. In het andere geval zoeken we verder tussen \(a\) en
\(mid\). We stellen dus b=mid. Ten slotte herberekenen we het middelpunt
\(mid\) en de waarde van de functie in \(mid\) en gaan we opnieuw de
iteratie in. De half-interval methode heeft dus haar naam niet gestolen:
de lengte van het interval waartussen we zoeken halveert in iedere
iteratie.
Wat kunnen we zeggen over de convergentie van deze methode naar een oplossing toe? Dat is lang niet altijd zeker. De functie moet vooreerst al een nulpunt hebben. Verder mag de gegeven precisie niet te klein zijn en zijn er bepaalde beperkingen omtrent de gladheid van de functie. We laten het volledige verhaal over aan een eventuele opvolgcursus over numerieke analyse.
Snelheid van dit algoritme¶
Opnieuw kunnen we ons de vraag stellen hoe snel dit algoritme naar een
oplossing convergeert. Laten we enkele uitvoeringen ervan tracen
d.m.v. een aangepaste versie van het algoritme die we traced_zero hebben genoemd. We hebben op strategische plaatsen een print statement toegevoegd aan de zero definitie.
Zoals gezegd halveert de lengte van het interval waarin gezocht in iedere iteratie. Na \(n\) keer itereren is de lengte van dat interval dus \(\frac{|b-a|}{2^n}\). Dus geldt \(\frac{|b-a|}{2^n} <\epsilon\) indien \(n > log_2(\frac{|b-a|}{\epsilon})\). Wat leren we hieruit? Voornamelijk dat door \(\epsilon\) kleiner te maken de breuk groter wordt en dus ook de logaritme groter wordt en dus \(n\) groter moet worden. Hoe preciezer we het resultaat willen, hoe meer we dus moeten itereren. Bovendien is het aantal iteraties logaritmisch afhankelijk van \(|b-a|\). Voor \(a=-3\) en \(b=4\) zien we dat \(log_2(\frac{|b-a|}{0.001}) = 12\) en we zien effectief ook dat het aantal iteraties groter is dan \(12\):
def traced_zero(f,a,b):
mid = (a+b)/2
val = f(mid)
print( "a=", "%9.6f"% a, "b=", "%9.6f"% b, "mid=", "%9.6f"% mid , "f(mid)=" , "%9.6f"% val)
while abs(val) > 0.001:
if val < 0:
a = mid
else:
b = mid
mid = (a+b)/2.0
val = f(mid)
print( "a=", "%9.6f"% a, "b=", "%9.6f"% b, "mid=", "%9.6f"% mid , "f(mid)=" , "%9.6f"% val)
return mid
traced_zero(test,-2,4)
a= -2.000000 b= 4.000000 mid= 1.000000 f(mid)= -8.000000
a= 1.000000 b= 4.000000 mid= 2.500000 f(mid)= -2.750000
a= 2.500000 b= 4.000000 mid= 3.250000 f(mid)= 1.562500
a= 2.500000 b= 3.250000 mid= 2.875000 f(mid)= -0.734375
a= 2.875000 b= 3.250000 mid= 3.062500 f(mid)= 0.378906
a= 2.875000 b= 3.062500 mid= 2.968750 f(mid)= -0.186523
a= 2.968750 b= 3.062500 mid= 3.015625 f(mid)= 0.093994
a= 2.968750 b= 3.015625 mid= 2.992188 f(mid)= -0.046814
a= 2.992188 b= 3.015625 mid= 3.003906 f(mid)= 0.023453
a= 2.992188 b= 3.003906 mid= 2.998047 f(mid)= -0.011715
a= 2.998047 b= 3.003906 mid= 3.000977 f(mid)= 0.005860
a= 2.998047 b= 3.000977 mid= 2.999512 f(mid)= -0.002929
a= 2.999512 b= 3.000977 mid= 3.000244 f(mid)= 0.001465
a= 2.999512 b= 3.000244 mid= 2.999878 f(mid)= -0.000732
2.9998779296875
Door de test in de while lus kunnen we echter ook geluk hebben. We
stoppen immers niet na het volledige doorlopen van \(n\) iteraties maar
wél indien onze \(f(mid)\) voldoende dicht bij nul ligt. Dat kan soms
verbazingwekkend snel zijn als we echt veel geluk hebben:
traced_zero(test,0,4)
a= 0.000000 b= 4.000000 mid= 2.000000 f(mid)= -5.000000
a= 2.000000 b= 4.000000 mid= 3.000000 f(mid)= 0.000000
3.0
Conclusie: zero genereert in het slechtste geval een logaritmisch
computationeel proces maar dat hoeft niet steeds het geval te zijn. We
zeggen dus dat het algoritme in \(O(log_2(\frac{|b-a|}{\epsilon}))\) is.
Herinner echter van hoofdstuk over recursieve functies dat logaritmische processen door
computerwetenschappers felbegeerd zijn omwille van het feit dat de
logaritme zo’n traag stijgende functie is.
Samenvatting¶
In dit hoofdstuk zijn we uitgegaan van declaratieve kennis en operationele kennis. Andere wetenschappen (wiskunde, fysica,…) leveren ons dikwijls uitstekende declaratieve kennis maar daar zijn we in de computerwetenschappen weinig mee. Computerwetenschappers hebben een constructieve operationele methode nodig om een oplossing voor een probleem op een systematische methode op te bouwen. Zó systematisch dat er helemaal geen menselijke intelligentie meer nodig is en het ook door een computer kan uitgevoerd worden.
Een eenvoudige methode bestaat erin om de operationele kennis te beschouwen als een zoekproces dat heel de tijd oplossingen suggereert en de declaratieve kennis gebruikt om te testen of een gesuggereerde oplossing daadwerkelijk een oplossing is. Dit werkt goed zolang de verzameling van mogelijke oplossingen relatief klein is en er voor iedere oplossing makkelijk een volgende oplossing kan gegenereerd worden. Het eist ook dat de oplossing exact voor te stellen is in een computergeheugen.
De tweede methode bestaat erin om de declaratieve kennis te gebruiken om een benaderingsmethode op te stellen die dan als operationele kennis kan dienen. Dit is noodzakelijk van zodra we niet zomaar een “volgende” oplossing kunnen vinden voor een oplossing of indien de oplossing misschien überhaupt niet bestaat in de wereld van eindige computergeheugens.