Funksjonsanalyse i Python - Programmering for undervisning

Introduksjon¶

Numerisk funksjonsanalyse

Kompetansemål¶

Denne aktiviteten passer kanskje bedre på VGS, men kan være en god aktivitet i 10. trinn som har disse aktuelle kompetansemålene:

utforske matematiske egenskapar og sammenhenger ved å bruke programmering
lese og forklare tekstbasert programmering

Forkunnskaper¶

For å kunne gjøre denne aktiviteten er det en fordel at man har vært innom følgende

funksjoner, nullpunkter og stigningstall i matematikk
andregradsfunksjoner og/eller andre funksjoner
variabler og operatorer i programmering
Å lage verditabeller i programmering (se egen aktivitet)
sløyfer/løkker i programmering
betingelser i programmering

Didaktisk tilnærming¶

Aktiviteten handler om å omsette kunnskap om funksjoner og verditabeller og betingelser til å kunne finne numeriske løsninger på funksjonenes nullpunkter og ekstremalpunkter. I hovedsak er det aktiviteten tenkt som en USE-MODIFY-CREATE aktivitet hvor man også blir presentert for noen grubleoppgaver underveis. Det er løsningsforslag til alle grubleoppgavene.

TODO¶

Lage flere oppgaver
Dele opp, la derivasjonsdelen være egen aktivitet

Å definere og plotte en funksjon i Python¶

Å definere en funksjon og lage en verditabell¶

I eksempelet under har vi definert en funksjon f som returnerer $4 - x^2$ , laget en verditabell for x-verdier og y-verdier:

# importerer nødvendige biblioteker
import numpy as np

# definerer funksjonen vår
def f(x):
    return 4 - x**2

x = np.linspace(-4,4,9)

# her sender vi x-tabellen inn i funksjonen, og lagrer funksjonsverdiene i en y-tabell
y = f(x)

print(x)
print(y)

[-4. -3. -2. -1.  0.  1.  2.  3.  4.]
[-12.  -5.   0.   3.   4.   3.   0.  -5. -12.]

Å plotte funksjoner¶

I programmet nedenfor har jeg lagt til følgende kode nederst i programmet.

plt.plot(x,y, color="red")
plt.axvline(x=0, color="black")
plt.axhline(y=0, color="black")
plt.show()

# importerer nødvendige biblioteker
import numpy as np
import matplotlib.pyplot as plt

# definerer funksjonen vår
def f(x):
    return 4 - x**2

x = np.linspace(-4,4,50)

# her sender vi x-tabellen inn i funksjonen, og lagrer funksjonsverdiene i en y-tabell
y = f(x)

plt.plot(x,y, color="red")      # lager et x-y plott med rød farge
plt.axvline(x=0, color="black") # tegner x-aksen i svart
plt.axhline(y=0, color="black") # Tegner y-aksen i svart
plt.show()                      # Viser plottet

Å finne nullpunkter i funksjonen¶

Det finnes flere metoder for å finne nullpunkter i en funksjon numerisk, vi skal her bare se på den enkleste måten hvor vi rett og slett går igjennom hele verditabellen og finner nullpunktene der y-tabellen har fortegnsskift. For å kunne gjøre dette må vi først se på en metode hvor vi går igjennom tabellen i en for-sløyfe. Denne metoden kalles for iterasjon.

Iterasjon av Arrays¶

Iterasjon av arrays innebærer å gå igjennom en tabell element for element. Følgende eksempel itererer et array som heter y (tenk y-tabell). Verdien til elementene hentes ut med y[i], der i er elementnummeret som telles opp i for-sløyfa. y.size vil returnere lengden av y-tabellen (her blir det 9, som definert i np.linspace(-4,4,9).

import numpy as np

def f(x):
    return 4 - x**2

x = np.linspace(-4, 4, 9)
y = f(x)
for i in range(y.size):
    print(y[i])

-12.0
-5.0
0.0
3.0
4.0
3.0
0.0
-5.0
-12.0

Løsningsforslag

I eksempelet er det lett å se nullpunktene, siden vi faktisk har med de eksakte løsningene i tabellen.
Forslag til kode:

import numpy as np

def f(x):
    return 4 - x**2

x = np.linspace(-4, 4, 9)
y = f(x)
for i in range(y.size):
    print(y[i])
    # Her trenger vi en betingelse for å sjekke
    # om vi har et nullpunkt i y-tabellen
    if y[i] == 0:
        # elementene i x og y tilsvarer hverandre, og indeksen i kan
        # nå brukes i x-tabellen for å skrive ut x-verdien
        print("nullpunkt for x = " + str(x[i]))

Hvis x-tabellen ikke “treffer” nullpunktene så må vi se på fortegnsskifter i tabellen, se neste avsnitt! Test selv om koden over i punkt 2 klarer å finne nullpunktene når lengden på linspace endres fra 9 til 10!

I de neste to avsnittene skal vi se på kode som kan finne alle nullpunkter, enten om de ligger i mellom to y-verdier, eller om vi faktisk har et reelt nullpunkt i tabellen. I tillegg er det lagt til kode som plotter funksjonen, og det er nyttig for da kan vi visuelt bekrefte at det finnes nullpunkter i definisjonsområdet vårt.

Nullpunkter ved fortegnsskift - metode 1¶

I verditabellen kan vi være heldige, og finne en y-verdi som er akkurat null, men som regel er dette ikke tilfelle. Vi kan imidlertid se på fortegnsskifte: dersom en y-verdi har forskjellig fortegn enn den neste, så har funksjonen krysset x-aksen. Følgende kode bruker np.sign() funksjonen for å finne fortegnet til y-verdien og sjekker for fortegnsskift.

Dersom vi har fortegnsskift må nullpunktet ligge MELLOM de tilhørende x-verdiene, derfor tar vi gjennomsnittet av disse to for å få et mer nøyaktig resultat. Dersom vi er så heldige at vi har et eksakt nullpunkt i tabellsettet, så får vi nødvendigvis ikke fortegnsskift derfor må vi teste for dette også i elif-blokken under.

# Importerer biblioteker
import numpy as np
import matplotlib.pyplot as plt
# Definerer funksjoner
def f(x):
    return 4 - x**2

# Antall elementer i x og y-tabellene
N = 100
# Lager x-tabell
x = np.linspace(-4,4,N)
# Lager y-tabell med tilhørende funksjonsverdier
y = f(x)
# Plotter grafen (ser om det faktisk finnes nullpkt)
plt.plot(x,y, color="red")
plt.axvline(x=0, color="black")
plt.axhline(y=0, color="black")
plt.show()
# Itererer y-tabellen for å finne fortegnsskift eller 0
for i in range(N-1):
    if np.sign(y[i]) != np.sign(y[i+1]):
        # Tar gjennomsnittet av x[i] og x[i+1]
        nullp = (x[i]+x[i+1])/2
        print("nullpunkt x = " + str(nullp))
    # sjekker også etter reelt nullpunkt
    elif y[i] == 0:
        print("nullpunkt x = " + str(x[i]))

nullpunkt x = -2.02020202020202
nullpunkt x = 2.0202020202020208

Nullpunkter ved fortegnsskift - metode 2 (snedig)¶

Man trenger ikke bruke sign(), men vi kan bruke fortegnsreglene for produkter til å lage en enklere test. Dersom y og neste y har samme fortegn, så vil produktet alltid være positivt. Dersom vi har et fortegnsskift mellom y og neste y, så vil produktet være negativt. Koden under viser hvordan dette kan gjøres:

# Importerer biblioteker
import numpy as np
import matplotlib.pyplot as plt
# Definerer funksjoner
def f(x):
    return 4 - x**2

# Antall elementer i x og y-tabellene
N = 100
# Lager x-tabell
x = np.linspace(-4,4,N)
# Lager y-tabell med tilhørende funksjonsverdier
y = f(x)
# Plotter grafen (ser om det faktisk finnes nullpkt)
plt.plot(x,y, color="red")
plt.axvline(x=0, color="black")
plt.axhline(y=0, color="black")
plt.show()
# Itererer y-tabellen for å finne fortegnsskift eller 0
for i in range(N-1):
    if y[i] * y[i+1] < 0:
        # Tar gjennomsnittet av x[i] og x[i+1]
        nullp = (x[i]+x[i+1])/2
        print("nullpunkt x = " + str(nullp))
    # sjekker også etter reelt nullpunkt
    elif y[i] == 0:
        print("nullpunkt x = " + str(x[i]))

nullpunkt x = -2.02020202020202
nullpunkt x = 2.0202020202020208

Nå har vi to versjoner av et verktøy som kan finne nullpunkter på alle funksjoner, i definisjonsområdet til x.

Oppgave

Nå kan du prøve ut andre funksjonsuttrykk, bare husk på at du må sjekke plottet og eventuelt justere definsjonsområdet i linspace() slik at nullpunktene faktisk eksisterer i datasettet. Prøv for eksempel $f(x)=cos⁡(x)+x$ Dette er en såkalt transcendental likning, som ikke er mulig å løse analytisk (ved regning!). Koden for å bruke cosinus i Python er:

y = np.cos(x) + x

Finne ekstremalpunkter i funksjonen¶

Vi kan bruke samme metode for å finne ekstremalpunkter i funksjonen som vi gjorde for nullpunkter, men betingelsen er forskjellig. For nullpunkter undersøkte vi fortegnsskift, for ekstremalverdier må vi undersøke for topp- eller bunnpunkter. Et eksempel for kode kan være:

# Importerer biblioteker
import numpy as np
import matplotlib.pyplot as plt

# Definerer funksjoner
def f(x):
    return 4 - x**2

# Antall elementer i x og y-tabellene
N = 109
# Lager x-tabell
x = np.linspace(-3,3,N)
# Lager y-tabell med tilhørende funksjonsverdier
y = f(x)
# Plotter grafen (ser om det faktisk finnes nullpkt)
plt.plot(x,y, color="red")
plt.axvline(x=0, color="black")
plt.axhline(y=0, color="black")
plt.show()
# Itererer y-tabellen for å finne topp eller bunnpunkter
for i in range(N-2):
    # sjekker for toppunkt
    if y[i] < y[i+1] and y[i+1] > y[i+2] :
        print("toppunkt x = " + str(x[i+1]) + ", y = " + str(y[i+1]))
    elif y[i] > y[i+1] and y[i+1] < y[i+2]:
        print("bunnpunkt x = " + str(x[i+1]) + ", y = " + str(y[i+1]))

toppunkt x = 0.0, y = 4.0

Løsning oppgave 2

Hvis man endrer N til 100 så er vil ikke testen:

if y[i] < y[i+1] and y[i+1] > y[i+2]

finne toppunktet, fordi det er to x-verdiene som ligger rundt symmetrilinjen ha samme verdi. Vi må legge til en ekstra test, både for topp- og bunnpunkt. Hvordan må vi endre if-testene? Løsningsforslag finnes i boksen nedenfor:

Løsning på problemet!!

I koden nedenfor er det lagt til en ekstra test for topp- og bunnpunkt. Hvordan fungerer de? Hvorfor regnes gjennomsnittet for x_verdien som vi har funnet for ekstremalpunktet? Legg også merke til at vi må bruke funksjonen f(x) for å finne toppunktet, for denne verdien finnes ikke i Arrayet y.

# Importerer biblioteker
import numpy as np
import matplotlib.pyplot as plt
# Definerer funksjoner
def f(x):
    return 4 - x**2

# Antall elementer i x og y-tabellene
N = 100 #NB
# Lager x-tabell
x = np.linspace(-3,3,N)
# Lager y-tabell med tilhørende funksjonsverdier
y = f(x)
# Plotter grafen (ser om det faktisk finnes nullpkt)
plt.plot(x,y, color="red")
plt.axvline(x=0, color="black")
plt.axhline(y=0, color="black")
plt.show()
# Itererer y-tabellen for å finne topp eller bunnpunkter
for i in range(N-2):
    # sjekker for toppunkt
    if y[i] < y[i+1] and y[i+1] > y[i+2]:
        print("toppunkt x = " + str(x[i+1]) + ", y = " + str(y[i+1]))
    elif y[i] == y[i+1] and y[i+1] > y[i+2]:
        x_snitt = (x[i]+x[i+1])/2
        print("toppunkt x = " + str(x_snitt) + ", y = " + str(f(x_snitt)))
    # sjekker for bunnpunkt
    elif y[i] > y[i+1] and y[i+1] < y[i+2]:
        print("bunnpunkt x = " + str(x[i+1]) + ", y = " + str(y[i+1]))
    elif y[i] == y[i+1] and y[i+1] < y[i+2]:
        x_snitt = (x[i]+x[i+1])/2
        print("bunnpunkt x = " + str(x_snitt) + ", y = " + str(f(x_snitt)))

Numerisk derivasjon (momentan vekstfart)¶

Med utgangspunkt i definisjonen for den deriverte, her er en Python-funksjon som finner den deriverte i ett valgt punkt på andregradsfunksjonen

Å finne vekstfarten i ett punkt¶

def f(x):
    return 4 - x**2

def deriverte(x,delta_x):
    return (f(x+delta_x)-f(x))/delta_x

x = 0
dx = 0.01

print("Den deriverte i punkt (" + str(x) + ", " + str(f(x)) + ")")
print("er: " + str(deriverte(x, dx)))

Den deriverte i punkt (0, 4)
er: -0.010000000000021103

Numeriske begrensinger¶

I matematikk bruker vi grenseverdien, når delta_x går mot 0, for å finne den deriverte analytisk. Det skulle da også være rimelig at det samme gjaldt når man regnet den deriverte numerisk med en datamaskin. Dette programmet analyserer hva forskjellen mellom den beregnede deriverte og den numeriske tilnærmingen blir for punktet (1, f(1)). Hva gjør programmet?

def f(x):
    return x**2 - 2

# analytisk derivert
def analytisk_derivert(x):
    return 2*x

def numerisk_derivert(f, x, delta_x):
    fder = (f(x + delta_x)-f(x))/delta_x
    return fder

delta_x = []

for i in range(1,16):
    dx = 10**(-i)
    delta_x.append(dx)
    numder = numerisk_derivert(f, 1, dx)
    ander = analytisk_derivert(1)
    feil = ander - numder
    print("Delta_x =", dx, " avvik:", feil)

Plotting av den deriverte¶

Dersom man har en verditabell for x-verdier, kan man finne en verditabell for den deriverte også – dette er nødvendig dersom man vil plotte den deriverte: Kjør følgende kode, og prøv å se hva den gjør

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 4 - x**2

def fderivert(x, delta_x):
    fder = (f(x + delta_x)-f(x))/delta_x
    return fder

dx = 0.00000001

x = np.linspace(-3,3,100)
y = f(x)
derivert = fderivert(x, dx)


plt.plot(x, y, color="red")
plt.plot(x, derivert, color = "blue")

plt.axvline(x=0, color="black")
plt.axhline(y=0, color="black")
plt.show()