Externa funktioner och numeriska metoder

Externa funktioner

Vi har tidigare sett hur man kan hantera enklare funktioner, med ett direkt samband mellan indata och utdata, genom att använda def kommandot.

Ifall man arbetar med större projekt vill man gärna samla alla funktioner som projektets olika delar har gemensamma i en separat fil och importerar dem när de behövs (tänkt på hur vi exempelvis importerat numpy när vi behövt funktionaliteten, i stället för att kopiera och klistra in all kod som behövs för att köra sin, cos, tan, ... osv.)

För att importera en funktion från en given fil skriver vi

from filnamn import funktionsnamn

Vill vi importera flera funktioner från filen skriver vi

from filnamn import funktionsnamn1, funktionsnamn2, funktionsnamn3, ...

Vill vi importera alla funktioner i filen skriver vi

from filnamn import *

Vi skapar nu ett illustrerande exempel. Skapa en fil som heter minafunktioner.py i samma mapp som denna Jupyter notebook. Klistra in följande kod i filen:

def addition(a,b):
    return a+b

def subtraktion(a,b):
    return a-b

def multiplikation(a,b):
    return a*b

def division(a,b):
    return a/b

Kör nu följande kod i rutan nedan (notera att man inte skriver .py när man anger filens namn)

from minafunktioner import addition
print(addition(4,6))

Ifall du kör subtraktion(4,6) i rutan ovan kommer du att få ett felmeddelande, eftersom funktionen inte är definierad. Vi har ju bara importerat första funktionen från filen.

För att importera alla funktioner skriver vi alltså import * i stället:

from minafunktioner import *
tal1 = 4
tal2 = 6
print(addition(tal1,tal2))
print(subtraktion(tal1,tal2))
print(multiplikation(tal1,tal2))
print(division(tal1,tal2))

Prova detta i rutan nedan.

Man vill oftast importera så lite som möjligt eftersom det påverkar tiden det tar för programmet att ladda samt kan många stora importer göra programmet långsammare vid körning.

Låt oss kombinera det vi lärt oss från tidigare pass och skapa en fil innehållande en funktion som beräknar det största tal i en lista en användare anger:

def maximum(lista):
    maxi = lista[0]
    for tal in lista[1:len(lista)]:
        if tal > maxi:
           maxi = tal

    return maxi

Denna kod tar emot en lista, sätter först maxi till första värdet i listan och stegar sedan igenom den. Ifall ett tal i listan är större än maxi så skrivs maxi över med det talet i stället. Först när loopen är färdig skickas maxi ut.

Spara nu denna funktion i filen minamaxmin.py och skriv sedan en egen funktion minimum som tar ut det minsta värdet i stället för det största och kör slutligen följande kod i rutan nedan:

from minamaxmin import *
print(maximum([2,55,1,0]))
print(minimum([2,55,1,0]))

Rekursiva funktioner

I detta avsnitt skall vi studera en teknik som påminner om hur vi tänker i induktionsbevis. Grundidén är att dela upp ett problem i liknande problem av mindre storlek.

Antag att vi vill beräkna och spara samtliga summor $$S(n) = \sum_{k=1}^{n} \frac{1}{k} $$ för alla heltal $1\leq n \leq 1000$.

Istället för att för varje $n$ summera alla tal i motsvarande summa så kan vi använda värdet på $S(n-1)$ för att beräkna $S(n)$ genom sambandet

$$S(n) = S(n-1) + \frac{1}{n}$$

dvs. Vi delar upp problemet med att beräkna $S(n)$ genom att lösadet analoga men mindre problemet att beräkna $S(n-1)$ samt utnyttja ett samband dem emellan. Detta kan implementeras i Python genom att t.ex. skapa följande externa funktion:

def harmtal(n):
    if n > 1:
        S = harmtal(n-1) + 1/n
    else:
        S = 1
    return S

Det som är speciellt med ovanstående funktion (mot de vi studerade i föregående avsnitt) är att funktionen harmtal anropar sig själv inuti funktionen. Den använder harmtal(n-1) för att beräkna harmtal(n). T.ex. behöver harmtal(374) värdet på harmtal(373) för sin beräkning, som i sin tur behöver värdet på harmtal(372). Under kodrutan nedan finns en beskrivning av vad som sker när koden körs för harmtal(374), rulla gärna ner till den och läs innan du går vidare.

Spara ovanstående kod i en fil vid namn harmtal.py och kör följande nedan:

import numpy as np
from harmtal import harmtal
# np.zeros(1000) skapar en vektor med tusen nollor
S_vektor = np.zeros(1000)

for tal in np.arange(1,1000,1):
    S_vektor[tal] = harmtal(tal)

print(S_vektor)

Med detta har vi skapat en vektor innehållande alla summor för $1\leq n \leq 1000$

Ifall vi skriver ut vad som sker när vi exempelvis kör harmtal(374) set det alltså ut så här:

Vi kör harmtal(374)

if 374 > 1: (=True)
   S = harmtal(373) + 1/374


som i sin tur kör harmtal(373)

if 373 > 1: (=True)
   S = harmtal(372) + 1/373


som i sin tur kör harmtal(372)

if 372 > 1: (=True)
   S = harmtal(371) + 1/372

...

som i sin tur kör harmtal(3)

if 3 > 1: (=True)
   S = harmtal(2) + 1/3


som i sin tur kör harmtal(2)

if 2 > 1: (=True)
   S = harmtal(1) + 1/2


som i sin tur kör harmtal(1)

if 1 > 1: (=False)
else:
    S = 1
return S (dvs. 1)


Nu returneras 1 ett steg upp och vi får

   S = 1 + 1/2
   return S (dvs. 1.5)


Nu returneras 1.5 ett steg upp och vi får

   S = 1.5 + 1/3
   return S (dvs. 1.8333333333333333)

...

Nu returneras 6.494764830964528 ett steg upp och vi får

   S = 6.494764830964528 + 1/372
   return S (dvs. 6.497453003007539)

Nu returneras 6.497453003007539 ett steg upp och vi får

   S = 6.497453003007539 + 1/373
   return S (dvs. 6.500133968154993)

Nu returneras 6.500133968154993 ett steg upp och vi får

   S = 6.500133968154993 + 1/374
   return S (dvs. 6.502807764946437)

Skälet till att vi här döpt funktionen till harmtal är för att summan i detta exempel har en viktig roll i vissa delar av matematiken och kallas då för harmoniska tal (Harmoniska tal är viktiga inom talteori och är nära relaterade till Riemanns zetafunktion).

Numeriska metoder

Många matematiska problem saknar lösning som kan beskrivas med exakta värden eller explicita uttryck med elementära funktioner t.ex. lösningen på en ekvation, värdet på en integral, lösningen på en differentialekvationer mm. I någon mening är det bara undantagsvis man kan lösa sådana problem "exakt". Ofta är det då konstruerade exempel i matematikkurserna för att bättre förstå olika begrepp, tekniker och samband.

I stället får man ofta nöja sig med försöka hitta approximativa lösningar med s.k. numeriska metoder. Ofta bygger dessa metoder på att man startar med en grov gissning på vad lösningen kan vara, och utifrån denna gissning sedan successivt komma närmare den exakta lösningen på problemet. Det blir en s.k. iterativ process där approx1 -> approx2 -> approx3 -> approx4 osv.

I en sådan iteration måste man också ha ett lämpligt villkor för att avbryta processen när man kommit tillräckligt nära den exakta lösningen. Eftersom man inte känner den exakta lösningen och därför inte kan avgöra (exakt) hur bra de approximativa lösningarna är så väljer man ofta att avbryta processen när skillnaden mellan två successiva approximationer är tillräckligt liten (mätt på lämpligt sätt).

Som exempel på numeriska metoder skall vi i nedanstående delavsnitt studera två olika sätt att hitta approximativa lösningar till ekvationer. Under ingenjörsstudierna kommer ni (sannolikt) stöta på många andra numeriska metoder t.ex. för att beräkna integraler och lösa differentialekvationer.

Fixpunktsiteration

Antag att vi vill lösa en ekvation på formen $x = F(x)$ för någon funktion $F(x)$. Antag också att vi på något sätt insett att det finns en lösning på ekvationen nära en viss punkt $x_1$ t.ex. genom att studera grafen till funktionen. Då kan vi bilda en talföljd genom:

$$x_n = F(x_{n-1}), ~~~ n = 2,3,4,...$$

Som schematiskt skulle kunna illustreras på följande sätt:

Låt oss se vad som händer med talföljden för några olika funktioner $F(x)$ och startvärden $x_1$.

$F(x) = \frac{x}{2} + \frac{1}{x}$

Med denna funktion och startpunkt $x_1=1$ får vi talföljden:

def fun(x):
    return x/2 + 1/x

x1 = 1 # startgissning
print(x1)
x2 = fun(x1)
print(x2)
x3 = fun(x2)
print(x3)
x4 = fun(x3)
print(x4)
x5 = fun(x4)
print(x5)

Detta ser ut att närma sig $\sqrt{2}$, eller hur? Lös nu ekvationen $ x = x/2 + 1/x $ med papper och penna för att se vad som händer.

Vi skulle också kunna använda en for-loop för att spara plats vid många iterationer:

import numpy as np
def fun(x):
    return x/2 + 1/x

x = 1 # startgissning
for n in np.arange(1,10,1):
    x = fun(x) # skicka in x, beräkna och skriv över
    print(x)

Vi kanske vill iterera tills skillnaden mellan två successiva tal är tillräckligt liten, t.ex. mindre än 0.00000001. Då passar det bättre med en while-loop, men får då vara försiktig så att inte snurran fastar i en evig loop, utifall talföljden inte skulle konvergera. Man kan då lägga in en räknare som försäkrar att snurran avbryts om den håller på för länge t.ex. efter 1000 varv.

import numpy as np
def fun(x):
    return x/2 + 1/x

x = 1 # startgissning
n = 1 # antal steg
tol = 0.00000001 # tolerans
while np.abs(x-fun(x)) > tol and n < 1000:
    print(n,x)
    x = fun(x) # skicka in x, beräkna och skriv över
    n = n + 1

$F(x) = \cos(x)$

Skriv själv kod som tar fram början av talföljden här och kör den i rutan nedan. Använd en while-loop för att iterera.

Närmar den sig något tal?

Tryck här för att se svaret

Här är det inte lika tydligt att följden verkligen närmar sig ett specifikt tal, men om man fortsätter iterationen ett tag så upptäcker man att den trots allt verkar närma sig ett värde runt ca. 0.7391. Även här kan man visa att talföljden kommer konvergera mot lösningen på ekvationen $ x = \cos(x) $, men vi har inte något analytiskt sätt att lösa denna ekvation exakt!

$F(x) = x/2-1/x$

Skriv själv kod som tar fram början av talföljden här och kör den i rutan nedan. Använd en for-loop för att iterera.

Närmar den sig något tal?

Tryck här för att se svaret

Här känns det inte som följden närmar sig något specifikt tal, oavsett hur länge man fortsätter iterationen. En talföljd av typen $x_n=F(x_{n-1})$ konvergerar alltså inte alltid mot något specifikt tal. Konvergensen beror både på vilken funktion som används och på startvärdet.

Låt oss nu se varför vissa talföljder konvergerar just mot lösningen på ekvationen $x = F(x)$.

Antag att vi har talföljden:

$$x_n = F(x_{n-1}), ~~~ n = 2,3,4,...$$

samt att talföljden $x_n$ närmar sig ett tal $x^*$ då $n$ går mot oändligheten.

Då får vi ju att:

$x_n\rightarrow x^*$ då $n\rightarrow \infty$

men även att

$F(x_{n-1}) \rightarrow F(x^*)$ då $n\rightarrow \infty$

Men det innebär ju i så fall att:

$$x^* = F(x^*)$$

Alltså att $x^*$ är en lösning till ekvationen $x = F(x)$. Ett sådant tal kallas en fixpunkt eftersom man får tillbaka samma tal som man stoppar in i funktionen!

Metoden vi lärt oss att använda här kallas därför fixpunktsiteration.

Här är en bild som illusterar vad som händer i ett steg med fixpunktmetoden

Intervallhalveringsmetoden

En annan metod för att lösa ekvationer är den s.k. intervallhalveringsmetoden. Här studerar vi i stället en ekvation på formen $f(x)=0$.

Antag att funktionen $f(x)$ är kontinuerlig och växlar tecken på intervallet $[a,b]$, dvs. antingen är $f(a)$ positiv och $f(b)$ negativ, eller vice versa.

Notera att $f(a)$ och $f(b)$ har olika tecken om och endast om $f(a)\cdot f(b) < 0$. Övertyga dig om detta innan du läser vidare.

Enligt satsen om mellanliggande värden vet vi att det måste finnas minst ett tal i intervallet sådant att $f(x)=0$ där. Övertyga dig om även om detta innan du läser vidare.

För att komma lite närmare detta nollställe kan vi undersöka funktionsvärdet i mittpunkten på intervallet dvs. $f(m)$, där $m=\frac{a+b}{2}$. Om nu $f(a)\cdot f(m) < 0$, så betyder det att ett nollställe också måste finnas mellan a och m. I annat fall måste ett nollställe finnas mellan m och b (eftersom det i så fall måste gälla att $f(b)\cdot f(m) < 0$).

Oavsett vilket så har vi på detta sätt halverat intervallet för vilket vi vet att det finns en rot. Vi kan sedan upprepa denna process tills intervallet är smalt nog så att vi är nöjda och väljer en punkt på det intervallet som lösning till ekvationen (differensen mellan intervallets ändpunkter kan ju exempelvis bara vara 0.000001 när vi är nöjda).

Här är en bild som illusterar vad som händer i ett steg med intervallhalveringsmetoden

Vi skriver kod som använder sig av rekursion för att lösa problemet. Vi lägger även till kommentarer omgivna av """ för att skapa texten som visas när man kör help(interhalv).

def interhalv(fun, I_in, tol):

    """
    interhalv - beräknar ett litet intervall innehållande
    nollställe till f(x).

    Syntax:
     I_ut = interhalv(f, I_in, tol)

    Argument:
     f - funktionen vi studerar
     I_in - ett fält med två värden som specificerar intervallet [a,b]
     tol - positivt tal som anger önskad bredd på intervallet
     innehållande nollstället.

    Utargument:
     I_ut - ett fält med två värden som specificerar 
     ett intervall innehållande nollstället.

    Exempel: 
      import numpy as np
      interhalv(np.cos,[1,3], 1e-5)
    """

    start = I_in[0]
    slut  = I_in[1]
    print(start)
    print(slut)
    print(fun)
    medelpunkt = (start+slut)/2

    if (slut-start) < tol:
        I_ut = I_in
    elif fun(start)*fun(medelpunkt) <= 0:
        I_ut = interhalv(fun, [start, medelpunkt], tol)
    else:
        I_ut = interhalv(fun, [medelpunkt, slut], tol)

    return I_ut

Kör ovanstående kod i rutan nedan och se hur det ser ut när vi kör help(interhalv). Lägg sedan till koden som behövs för att använda interhalv för att hitta ett nollställe till $f(x)=x^2-2$.

Algoritmer, flödesdiagram och pseudokod

För att lösa mer komplicerade problem så behöver man en strategi. Det finns olika sätt att angripa ett problem och det kan vara svårt att veta i vilken ände man skall börja. Ofta kan det vara en bra idé att göra förenklingar eller titta på ett närliggande enklare problem för att bättre förstå var svårigheterna finns. Det kan också ofta vara en bra idé att dela upp problemet i mindre delar, och angripa varje del för sig.

Ett vanligt misstag är att börja koda kommandon innan man har en helhetsbild av hur problemet skall lösas. Vi skall i detta avsnitt istället titta på några sätt på vilket man kan angripa problem och beskriva lösningar utan att arbeta med specifika kommandon. Målet är att finna en beskrivning av hur ett visst problem skall lösas (med ett ändligt antal steg), en s.k. algoritm.

Flödesdiagram

Det kan ibland vara användbart med s.k. flödsdiagram för att hitta och beskriva lösningar på problem. Det är en sorts grafisk bild av ett program som kan ge en tydligare överblick av hur ett program är uppbyggt. Vi skall här bara titta på något enkelt exempel på flödesdiagram. För den intresserade finns mycket mer att läsa om detta på internet (sök på flödesdiagram eller flowchart).

Som exempel på ett flödesdiagram, låt oss åter titta på programmet maximum som vi skapade tidigare.

Det finns inga absoluta regler för hur flödesschemat skall se ut, det viktiga är att tydliggöra de olika stegen. Det finns dock konventioner/överenskommelser om vilka former boxarna skall ha för olika aktiviteter. Vi går inte närmare in på det här, det får vänta tills en mer omfattande kurs i programmering.

Gör nu ett flödesdiagram för något av de program vi hittills tittat på i kursen. Det spelar inte så stor roll vilket program du väljer och det är inte så noga med att rita alla boxar på rätt sätt. Vi har ju inte gjort någon djupare analys av detta i ovanstående text så här är det bara tänkt att du skall reflektera lite över det som nämnts.

Pseudokod

Ett annat sätt att arbeta fram ett program är att använda sig av s.k. pseudokod (som direktöversatt betyder falsk eller skenbar kod). Det är också en skiss av programstrukturen, ofta i form av en punktlista, men till skillnad mot en algoritmisk beskrivning med löpande text så är den lite mer kortfattad.

Liksom för flödesprogram, finns det ingen tydlig standard för hur man skall skriva pseudokod, så det finns stor frihet i hur den kan se ut. Pseudokod innehåller normalt inte exekverbar kod men kan innehålla mer allmänna programmeringsbegrepp som INPUT, OUTPUT, IF, ELSE, THEN, FOR, WHILE, BEGIN, END mm.

Tanken är dock att koden inte skall följa syntax i något specifikt programmeringsspråk utan kunna tolkas på ett mer övergripande sätt (och implementeras i vafritt språk). I pseudokod för program av mer matematisk art är det också vanligt med matematiska begrepp, formler och beteckningar.

Pseudokod för maximumfunktionen skulle kunna se ut så här:

Detta kan jämföras med en beskrivning av algoritmen som istället skulle kunna se ut så här;

• Be användaren om en vektor.
• Det största värdet sätts initialt till det första talet i vektorn.
• Använd en snurra för att succesivt jämföra övriga tal i vektorn med det hittills största. Om ett större tal hittas så sätts detta värde till det största.
• Efter att snurran gått igenom alla tal i vektorn så skall det största värdet skrivas ut.