Numpy: labb 2

NumPy-laboration 2

Mål

Efter laborationen ska du

kunna hantera polynom och matriser i NumPy,
kunna beräkna integraler, samt
kunna beräkna minima, maxima och nollställen.

I laborationens temauppgift kommer du även att få lära dig mer om Mandelbrotmängden.

Kom igång

I den här labben ska vi använda NumPy paketet polynomial, matplotlib och Pythonbiblioteket SciPy. Skapa ett nytt skript lab2.py och börja med att importera paketen med lämpliga aliaser:

import numpy as np import matplotlib.pyplot as plt import numpy.polynomial.polynomial as polynomial import scipy.optimize as optimize

Polynom

Det finns ingen generell analytisk metod för att bestämma rötterna till ett femtegradspolynom. Men med numeriska verktyg, kan du bestämma dessa rötter relativt enkelt.

Betrakta följande polynom:

Bestäm polynomets rötter.

Läs om polyroots

Kontrollera att polynomets värden i dessa nollställen verkligen är noll (eller åtminstone mycket nära).

Kontrollera även med ett par andra x-värden, och bekräfta att polynomet har ett annat värde än noll där.

Använd polyval för att räkna ut polynomets värde i rötterna.

y = polynomial.polyval(roots, poly) räknar ut polynomets värde för alla x-värden i roots för polynomet med koefficienterna poly. Det går också bra att använda en skalär som argument för att räkna ut polynomets värde för ett specifikt x-värde, så här: y = polynomial.polyval(x, poly) där x är ett tal.
Du kan även bestämma ett polynom för givna rötter.

Bestäm ett polynom px som är av lägsta möjliga grad och har följande rötter:

Läs om polyfromroots

Vi kan naturligtvis kontrollera rötterna på samma sätt som i föregående uppgift. Här ska vi emellertid göra på ett annat sätt.

Plotta polynomets värden i intervallet –50 ≤ x ≤ 50. Kontrollera att nollställena verkar ligga rätt.

För att räkna ut en vektor y = p(x) för en vektor x av värden, använd polyval så här: y = polynomial.polyval(x, poly).

Saknas det nollställen?

np.linspace ger vanligtvis bara 50 värden, då kanske vi missar något nollställe helt! Öka antalet värden genom att ange ett argument num, till exempel för att skapa 1000 värden mellan -50 och 50: np.linspace(-50, 50, num = 1000).

Hur kan vi justera diagrammets y-axel för att nollställena ska kunna urskiljas?

Till exempel plt.ylim(-10, 10)

Integraler

Funktionen ft definieras som följer:

ft

Vi söker nu värdet av följande integral:

S

0

ft dt

Beräkna integralvärdet S numeriskt med trapetsmetoden i NumPy.

Skapa funktionen f med parametern t med ett lambda-uttryck: f = lambda t: ....

np.trapezoid(f(x), x).

NumPy kommer att anropa din funktion med en vektor som parameter. Funktionen måste därför använda elementvisa operationer, till exempel np.sin(...).

Plotta funktionen och bedöm om det beräknade integralvärdet verkar rimligt.

Hjälp! Bara den förra bilden syns!

plt.show() blockerar resten av exekveringen tills fönstret har stängts, stäng fönstret så kör koden vidare.
Kan jag visa många bilder samtidigt?

Du kan visa många bilder så här:
1. Lägg till plt.figure() före varje ny plot.
2. Ersätt plt.show() med plt.show(block = False) för alla förutom den sista figuren.
(Klicka inte förrän du avgjort om ditt svar är rimligt.)

S =

Minima, maxima och nollställen

Vi vill veta hur funktionen gx beter sig i intervallet 2 ≤ x ≤ 3.

gx

cos e^x

1 – x

Beräkna alla lokala minima för gx i det givna intervallet. Notera både x-värde(n) och motsvarande funktionsvärde(n). Använd scipy.optimize.minimize för att minimera funktionen.

Kom ihåg att plotta funktionen för att få en uppfattning om var lokala maxima och minima finns. Det är praktiskt med ett rutnät, och ett sådant får du med plt.grid().

Vad returnerar minimize?

Utskriften från minimize innehåller mycket information, vilket värde är intressant och hur kan vi läsa ut det?

Om vi sparar resultatet från funktionsanropet i variabeln result = optimize.minimize(...) kan vi få ut vektorn med x-värde(n) med result.x och funktionsvärdet med result.fun .

Syns det inget in din plot?

Syns det inget? Tänk på att din funktion anropas med en vektor som parameter.
Nu söker vi både minima och nollställen.

Använd scipy.optimize paketet för att beräkna alla nollställen och minima (x-värde(n) och motsvarande funktionsvärde(n)) för funktionen hx nedan.

hx

Läs om optimize.fsolve.

optimize.fsolve anropas med funktionen h och en gissning g som argument. För att hitta bra gissningar kan det vara bra att plotta polynomet som i tidigare uppgifter.

Minimum:
h =

Nollställen:
h =
h = 0

Mer om matriser

Bandmatrisen M nedan har storleken . Matrisen har värdet på diagonalen och på subdiagonalerna.

M

Skapa M med hjälp av NumPy-funktion eye, notera att den valfria parametrarn k anger diagonalen.

Skriv M som summan av tre termer, där varje term beräknas med eye(...) + eye(..., k = -1) + ... etc.

Du kan skapa en 4⨯4- matris med 7:or på den undre diagonalen så här: np.eye(4, k = -1) * 7.

Bestäm determinanten för M.

Determinanten är .

Inför programmeringsuppgiften Mandelbrotmängden (se nedan) behöver vi en matris med komplexa tal inom ett visst talområde. Nedan visas en figur över ett sådant (rektangulärt) talområde och en motsvarande matris. Vi kan se det som att matrisen spänner upp ett "nät" i det komplexa planet.

För enkelhets skull är matrisen kvadratisk (lika många rader som kolonner).

–2+i	–1.25+i	–0.5+i	0.25+i	1+i
–2+0.5i	–1.25+0.5i	–0.5+0.5i	0.25+0.5i	1+0.5i
–2	–1.25	–0.5	0.25	1
–2–0.5i	–1.25–0.5i	–0.5–0.5i	0.25–0.5i	1–0.5i
–2–i	–1.25–i	–0.5–i	0.25–i	1–i

Vi behöver nu en funktion, complexmat, som skapar en sådan matris. Funktionen ska ha tre parametrar:

N: Matrisens storlek (matrisen får storleken N ⨯ N och är alltså kvadratisk)
z0: Områdets övre vänstra punkt i det komplexa talplanet
z1: Områdets nedre högra punkt i det komplexa talplanet

Så här kan funktionen användas för att skapa en 5⨯5-matris för området från –2+i till 1–i, precis som på bilden ovan (kom ihåg att Python använder j för imaginära tal):

m = complexmat(5, -2+1j, 1 - 1j) print(m) [[-2.00+1.00j -1.25+1.00j -0.50+1.00j 0.25+1.00j 1.00+1.00j] [-2.00+0.50j -1.25+0.50j -0.5+0.50j 0.25+0.50j 1.00+0.50j] [-2.00+0.00j -1.25+0.00j -0.5+0.00j 0.25+0.00j 1.00+0.00j] [-2.00-0.50j -1.25-0.50j -0.5-0.50j 0.25-0.50j 1.00-0.50j] [-2.00-1.00j -1.25-1.00j -0.5-1.00j 0.25-1.00j 1.00-1.00j]]

Kommentar: för att betämma hur många decimaler som ska visas i utskriften av en np.array kan vi använda set_printoptions, till exempel för att visa exakt två decimaler som ovan: np.set_printoptions(precision = 2, floatmode='fixed').

En nästan färdig version av complexmat finns nedan. Skapa en ny fil och spara den som mandelbrot.py. Du kan till exempel använda tangentbordskommandona (Windows, Linux) ctrl-N eller (Mac) cmd--N och (Windows, Linux) ctrl-S eller (Mac) cmd--S.

def complexmat(N, z0, z1): # skapar en N x N-matris med komplexa tal a + bi # där re(z0) <= a <= re(z1) och im(z1) <= b <= im(z0) # (jämnt fördelade i matrisen) xs = np.linspace(z0.real, z1.real, num = N) # realdelar ys = np.linspace(z0.imag, z1.imag, num = N) # imaginärdelar # skapa två matriser med real- respektive imaginärdelar [X, Y] = np.meshgrid(xs, ys) # matrisen X innehåller resultatets realdelar # matrisen Y innehåller resultatets imaginärdelar return X # <<< denna rad behöver utökas end

Du kan enkelt kopiera koden ovan genom att placera muspilen över den och trycka (Windows, Linux) ctrl-C eller (Mac) cmd--C.

Kör funktionen för exemplet ovan. Resultatet stämmer inte riktigt med det önskade. Gör klart funktionen. (En liten förändring behövs på den utpekade raden ovan, där resultatet skapas.)

Du kan läsa mer om funktionen meshgrid här (länk). Sidans första exempel liknar vårt.

Om du kör funktionen i omodifierat skick, för samma indata som ovan, får du följande resultat:

m = complexmat(5, -2+1j, 1-1j) print(m) [[-2.00 -1.25 -0.50 0.25 1.00] [-2.00 -1.25 -0.50 0.25 1.00] [-2.00 -1.25 -0.50 0.25 1.00] [-2.00 -1.25 -0.50 0.25 1.00] [-2.00 -1.25 -0.50 0.25 1.00]]

Vad är det som saknas? Hur kan du förändra den utpekade raden i complexmat ovan så att resultatet blir rätt?

Resultatet saknar imaginärdel.

return X + 1j * Y

Tema: Mandelbrotmängden

Bakgrund

Mandelbrotmängden är den mängd av komplexa tal c för vilka följande serie konvergerar:

z₀

z_k+1

z2
k + c

Vi bildar alltså ett nytt värde i serien genom att kvadrera det föregående samt lägga till c. Frågan är, för varje sådant c, huruvida serien konvergerar eller divergerar.

I de följande deluppgifterna ska vi försöka förstå något om Mandelbrotmängden. I ovanstående enkla samband döljer sig ett fantastiskt, oändligt detaljerat mönster, som vi ska titta närmare på.

Att testa om ett tal tillhör mängden

Vi kan undersöka om ett visst tal ingår i Mandelbrotmängden med hjälp av definitionen ovan. Två exempel:

Tillhör det komplexa talet 0.1i mängden? För att ta reda på det sätter vi c0.1i, och får då en serie som börjar så här:

z₀

c0.0+0.1i

z₁

z2
0 + c–0.01+0.1i

z₂

z2
1 + c–0.0099+0.098i

z₃

z2
2 + c–0.009506+0.09806i

z₄

z2
3 + c–0.009525+0.098136i

z₅

z2
4 + c–0.009540+0.098130i

z₆

z2
5 + c–0.009539+0.098128i

z₇

z2
6 + c–0.009538+0.098128i

z₈

z2
7 + c–0.009538+0.098128i

Denna serie konvergerar mot –0.009538+0.098128i. Vi kan därför dra slutsatsen att 0.1i ingår i Mandelbrot-mängden.
På samma sätt kan vi pröva om 126 tillhör mängden genom att sätta c126. Då får vi z som följer:

z₀

c126

z₁

z2
0 + c16002

z₂

z2
1 + c256064130

z₃

z2
2 + c65568838672657024

z₄

z2
3 + c4299272604880923324903695936536702

Serien divergerar uppenbarligen mot +∞. Talet 126 tillhör alltså inte Mandelbrot-mängden.

Att visualisera Mandelbrot-mängden

Det är inte särskilt praktiskt att själv undersöka varje möjligt tal på det här viset. Istället skulle vi vilja visualisera mängden på något sätt. Vi söker en bild av det komplexa talplanet, där vi med färg visar vilka komplexa värden som tillhör Mandelbrot-mängden, och vilka som inte gör det.

Du kommer nu att använda NumPy och matplotlib för att skapa en sådan bild.

Innan vi ger oss på den riktiga Mandelbrotmängden, låt oss först begrunda följande (något enklare) talserie:

z₀

c

z_k+1

z2
k

Den enda skillnaden är att vi utelämnar termen c när vi beräknar nästa värde i talserien.

Undersök för vilka värden på c serien konvergerar. Kanske inser du svaret direkt, men låt oss ändå pröva några värden på c. Låt exempelvis c = 0.5 och använd Python-repl:en för att testa följande (skriv kommandot python i terminalen för att start en repl):

z = 0.5 z = z * z z 0.25 z = z * z

Upprepa beräkningen z = z*z tre–fyra gånger så du ser värdet konvergera. Tips: du kan använda pil upp-tangenten för att hämta kommandon du skrivit tidigare och se nuvarande värde på z genom att skriva z på en egen rad.

Gör nu samma undersökning för några startvärden till, som

c

1.1

c

0.5 + 0.5i

Kanske har du nu en idé om för vilka värden serien konvergerar. (Det har med enhetscirkeln att göra.) Vi kan emellertid inte beräkna

lim
k → ∞ z_k

numeriskt på det här viset (inom ändlig tid). Istället behöver vi känna igen konvergens på något annat sätt.
Låt oss göra det lite enklare för oss, och anta att följden divergerar om följande gäller:

z_k > 2

för något k, 0 ≤ k ≤ 100

Annorlunda formulerat: vi antar att serien konvergerar om seriens första 100 värden är tillräckligt små. (Det är naturligtvis ingen ersättning för konvergens i allmänhet, men fungerar bra här.) Med z avses absolutbeloppet av z.

Skapa en ny funktion i filen mandelbrot.py. Funktionen kan heta converge och ta ett startvärde (c) som parameter. Funktionen ska returnera antalet iterationer som krävdes för att avgöra konvergens, dock max 100 (k enligt ovan).
Funktionen ska alltså upprepa beräkningen så länge följande villkor är uppfyllda:
- Det senaste z-värdets absolutvärde är mindre än eller lika med 2 och
- Antalet gjorda iterationer är mindre än 100.
Du behöver alltså en loop.
Testa din nya funktion med några värden, och se att den verkar bete sig riktigt:

converge(0.5+0.5j) #ska ge 100 converge(1.1) # ska ge 3

För c0.5 + 0.5i hamnar alltså seriens samtliga första 100 värden innanför gränsen, och för detta c anser vi därmed att serien konvergerar. För c1.1 är det däremot bara de tre första värdena som (till absolutbeloppet) är mindre än 2. Serien divergerar alltså för det startvärdet.

Om resultaten från din converge-funktion avviker lite grand från det angivna (exempelvis 99 istället för 100, eller 4 istället för 3), så är det helt okej. Det viktiga här är att vi kan se skillnad på om serien konvergerar eller divergerar för ett visst värde, och något om hur snabbt den divergerar (om den gör det).
Nu ska vi ta fram en matris med sådana funktionsvärden. Vi utgår från en matris m, som ges av funktionen complexmat (som du gjorde färdigt tidigare under laborationen) och omfattar området från –2+i till 1–i:

m = complexmat(5, -2+1j, 1 - 1j) print(m) [[-2.00+1.00j -1.25+1.00j -0.50+1.00j 0.25+1.00j 1.00+1.00j] [-2.00+0.50j -1.25+0.50j -0.5+0.50j 0.25+0.50j 1.00+0.50j] [-2.00+0.00j -1.25+0.00j -0.5+0.00j 0.25+0.00j 1.00+0.00j] [-2.00-0.50j -1.25-0.50j -0.5-0.50j 0.25-0.50j 1.00-0.50j] [-2.00-1.00j -1.25-1.00j -0.5-1.00j 0.25-1.00j 1.00-1.00j]]

Vi ska nu beräkna en annan matris, v. Varje matriselement i v ges av funktionen converges värde för motsvarande komplexa tal i matrisen m.

Skriv python-kod som beräknar en matris v som har samma dimensioner som m, men istället består av motsvarande funktionsvärden från converge.

För 5⨯5-matrisen m ovan ska matrisen v ska se ut så här:

print(v) [[0 1 3 5 2] [0 2 100 100 3] [1 2 100 100 100] [0 2 100 100 3] [0 1 3 5 2]]

Här ser vi alltså, för ett antal olika komplexa värden, huruvida serien konvergerar eller inte. Värdet 100 betyder att serien konvergerat, värdet 0 att den divergerat omedelbart, och värden däremellan att den divergerat efter ett antal iterationer.

Beroende på exakt hur du definierat din converge-funktion kan värdena skilja sig något (±1 i några få positioner). Sådana små avvikelser är inget problem här; det viktiga är att vi kan avgöra om serien konvergerar för de givna värdena.

Kan jag använda en for-loop som itererar över alla element?

Om vi skriver vår egen for-loop kan vi inte använda NumPys optimering, så koden blir mycket långsammare, det går 2-3 gånger snabbare om vi vektoriserar funktionen med np.vectorize så här: vConverge = np.vectorize(converge).

Hur gör jag för att anropa en funktion elementvis i en matris?

Anropa en funktion elementvis på en matris

Vi vill ju iterera olika många steg beroende på elementet i matrisen, så vi kan inte beräkna alla värden på en gång som vi kunde med matematiska funktioner i förra labben (np.cos(), np.exp() osv). Funktionen np.vectorize skapar en vektoriserad version av en funktion, till exempel: vfunc = np.vectorize(func). Sen kan vfunc anropas med matrisen m som argument, v = vfunc(m), vilket returnerar en matris med funktionen applicerad elementvis.

För att anropa converge på hela matrisen kan vi först vektorisera den med vConverge = np.vectorize(converge) och sen beräkna v med v = vConverge(m) där m är den complexa matrisen ovan.

I exemplet ovan kan vi se att serien konvergerar för sju av värdena i m. Vilka då?

Med funktionen plt.imshow() kan vi visa en matris som en bild, där varje ruta i bilden motsvarar matrisens värde. I vårt fall får vi alltså en bild av för vilka värden serien konvergerar. Pröva:

plt.imshow(v, aspect = 2/3) plt.show()

Eftersom matrisen är kvadratisk men täcker ett område som har längd 3 i x-led (realdelen) och 2 i y-led (imaginärdelen), kan vi dra ut den horisontellt med argumentet aspect = 2/3 som ändrar bildens aspect ratio.
Bilden i föregående uppgift är minst sagt grovkornig, och det krävs en smula fantasi för att urskilja enhetscirkeln. Kan du ändå ana den?

Skapa en ny matris m, denna gång av storleken 200⨯200, med hjälp av complexmat. Utgå från samma område i det komplexa planet som tidigare.

Skapa också en ny matris v av värden, som du gjorde i uppgift D10.

Tar det väldigt lång tid?

Om ditt program fastnar i en oändlig loop, kan du bryta funktionen genom att klicka i terminalfönstret och trycka Ctrl-C.

Visa matrisen v som en bild, precis som i uppgift D10. Nu syns enhetscirkeln tydligare. Runt cirkeln finns en "aura" av lite avvikande färg. Vad beror den på?
Nu ska vi äntligen undersöka Mandelbrotmängden. Denna mängd definieras av de komplexa tal för vilka följande serie konvergerar:

z₀

c

z_k+1

z2
k + c

Som du ser är den mycket lik den talserie vi undersökt, frånsett termen c.

Modifiera funktionen converge så att den undersöker konvergens för uttrycket ovan. Beräkna om matrisen v, visa den som en bild och förundras.

Öka nu storleken på m till 1000⨯1000 med hjälp av complexmat. Utgå från samma område i det komplexa planet som tidigare. Nu tar det längre tid att rita upp figuren men den blir mer detaljrik.

Den ursprungliga talserien konvergerar för z < 1. Mandelbrotmängden är inte lika lätt att ringa in. Visst är den märklig?

Vi kan byta färgskala på figuren genom att lägga till ett argument cmap till imshow(): plt.imshow(v, aspect = 2/3, cmap = plt.get_cmap('ocean')) där 'ocean' är namnet på en av många colormaps. Välj ut en snygg colormap till din figur.

Vilka colormaps är snyggast?

Prova till exempel hot, inferno, RdPu, Spectral, cubehelix ...

Detta är en fraktal med oändliga, fantastiska mönster. Vi kan alltså "zooma" in i figuren och finna nya detaljer. (I praktiken sätter Pythons beräkningsprecision, även om den är god, gränsen för hur långt vi kan zooma.)

När man tittar närmare verkar faktiskt mönstret innehålla små kopior av sig självt! Pröva med något av följande komplexa områden (som alla tillhör den ursprungliga figuren):

Från –0.7+0.7i till –0.5+0.6i

Från –1.4+0.48i till –1.1+0.24i

Från –1.8+0.03i till –1.5–0.03i

Från –0.18–1.02i till –0.14–1.05i

Skapa en ny m-matris (gärna enligt något av förslagen ovan), beräkna om v och rita om bilden.

m = complexmat(1000, -0.7+0.7j, -0.5+0.6j)

m = complexmat(1000, -1.4+0.48j, -1.1+0.24j)

m = complexmat(1000, -1.8+0.03j, -1.5-0.03j)

m = complexmat(1000, -0.18-1.02j, -0.14-1.05j)
När du kommit hit är du klar med kursens andra NumPy-laboration. Läs igenom uppgifterna igen och se till att du förstår dina lösningar inför redovisningen.

Mer om Mandelbrot-mängden

Intresset för Mandelbrotmängden tog fart under 1970- och 1980-talen, då tillgången till datorer med färggrafik gjorde det möjligt för fler att undersöka mängden.

Många fascineras av Mandelbrotmängden, och på nätet finns många förslag på områden att undersöka. Om du är intresserad, sök exempelvis på mandelbrot interesting coordinates. Man kan även hitta fantastiska animationer, som exempelvis här till höger.