Matlab-laboration 2
EDAA55, LTH

Mål

Efter laborationen ska du

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

    Polynom

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

    Betrakta följande polynom:

    Bestäm polynomets rötter.

    Du såg ett användbart exempel i följande film.

    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.

  2. 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:

    Man 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.

    Kanske behöver du justera diagrammets y-axel för att nollställena ska kunna urskiljas.

    help axis

    3. Integraler

  3. 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 Matlab.

    help integral

    Matlab kommer att anropa din funktion med en vektor som parameter. Funktionen måste därför använda elementvisa operationer.

    Du såg ett användbart exempel med integral i följande film.

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

    (Klicka inte förrän du avgjort om ditt svar är rimligt.)

    S =

    4. Minima, maxima och nollställen

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

    gx
    cos ex
    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).

    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 grid on.

    Syns det inget in din plot?

    Syns det inget? Tänk på att din funktion anropas med en vektor som parameter.

    Vi vill beräkna ett funktionsvärde för varje vektorelement. Då är det elementvisa operationer vi behöver använda.

  5. Nu söker vi både minima och nollställen.

    Använd Matlab 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

    Minimum:
    h =

    Nollställen:
    h =
    h = 0

    6. Mer om matriser

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

    M

    Skapa M med hjälp av Matlab-funktionerna eye och diag.

    Du kan skriva help diag för mer information. Se även laboration 1.

    Skriv M som summan av tre termer, där en term beräknas med eye och de andra två med hjälp av diag.

    Du kan skapa en radvektor med, säg, sju ettor med ones(1,7).

    Bestäm determinanten för M.

    Determinanten är .

  7. Inför programmeringsuppgiften 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 Matlab-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:

    m = complexmat(5, -2+i, 1-i) -2.00+1.00i  -1.25+1.00i  -0.50+1.00i   0.25+1.00i   1.00+1.00i -2.00+0.50i  -1.25+0.50i  -0.50+0.50i   0.25+0.50i   1.00+0.50i -2.00+0.00i  -1.25+0.00i  -0.50+0.00i   0.25+0.00i   1.00+0.00i -2.00-0.50i  -1.25-0.50i  -0.50-0.50i   0.25-0.50i   1.00-0.50i -2.00-1.00i  -1.25-1.00i  -0.50-1.00i   0.25-1.00i   1.00-1.00i

    En nästan färdig version av complexmat finns nedan. Skapa en funktionsfil (m-fil) i Matlab, exempelvis genom att klicka på knappen New→Function. Spara filen som complexmat.m.

    function result = 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 = linspace(real(z0), real(z1), N);  % realdelar     ys = linspace(imag(z0), imag(z1), N);  % imaginärdelar       % skapa två matriser med real- respektive imaginärdelar     [X, Y] = meshgrid(xs, ys);       % matrisen X innehåller resultatets realdelar     % matrisen Y innehåller resultatets imaginärdelar       result = 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. Notera dock att Matlab på LTH:s Linux-datorer har ett lite ovanligt tangentbordskommando för att klistra in (ctrl-Y).

    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+i, 1-i) -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.

    result = X + i * Y;

    8. Tema: Mandelbrotmängden

    Bakgrund

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

    z0
    c
    zk+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å med hjälp av Matlab.

    Att testa om ett tal tillhör mängden

    Vi kan undersöka om ett visst tal ingår i Mandelbrot-mä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:

      z0
      c0.0+0.1i
      z1
      z2
      0       + c      –0.01+0.1i
      z2
      z2
      1       + c      –0.0099+0.098i
      z3
      z2
      2       + c      –0.009506+0.09806i
      z4
      z2
      3       + c      –0.009525+0.098136i
      z5
      z2
      4       + c      –0.009540+0.098130i
      z6
      z2
      5       + c      –0.009539+0.098128i
      z7
      z2
      6       + c      –0.009538+0.098128i
      z8
      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:

      z0
      c126
      z1
      z2
      0       + c      16002
      z2
      z2
      1       + c      256064130
      z3
      z2
      2       + c      65568838672657024
      z4
      z2
      3       + c      4299272604880923324903695936536702

      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 Matlab för att skapa en sådan bild.

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

    z0
    c
    zk+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 Matlab för att testa följande:

    z = 0.5 z = z * z 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.

    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 → ∞     zk

    numeriskt på det här viset (inom ändlig tid). Istället behöver vi känna igen konvergens på något annat sätt.

  9. Låt oss göra det lite enklare för oss, och anta att följden divergerar om följande gäller:

    zk > 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.

    Skriv en funktionsfil (m-fil) i Matlab (New→Script). 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 (till skillnad från i den förra funktionsfilen du skrev).

    Testa din nya funktion med några värden, och se att den verkar bete sig riktigt:

    converge(0.5+0.5i) 100 converge(1.1) 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).

  10. 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+i, 1-i) -2.00+1.00i  -1.25+1.00i  -0.50+1.00i   0.25+1.00i   1.00+1.00i -2.00+0.50i  -1.25+0.50i  -0.50+0.50i   0.25+0.50i   1.00+0.50i -2.00+0.00i  -1.25+0.00i  -0.50+0.00i   0.25+0.00i   1.00+0.00i -2.00-0.50i  -1.25-0.50i  -0.50-0.50i   0.25-0.50i   1.00-0.50i -2.00-1.00i  -1.25-1.00i  -0.50-1.00i   0.25-1.00i   1.00-1.00i

    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 Matlab-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:

    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.

    Hur använder jag en funktion på alla matrisvärdena på en gång?

    Använd funktionen arrayfun för att skapa värdena. Pröva gärna help arrayfun.

    När man anger en funktionsfil som parameter till en annan funktion använder man @, exempelvis @converge. (Se Matlab Primer, "Function Handles".)

    Du utgår här från den matris m du skapat med din complexmat. Med hjälp av denna kan du få resultatet genom att göra

    v = arrayfun(@converge, m);

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

    Med Matlab-funktionen imagesc kan man 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:

    imagesc(v)
    Blir allt bara blått?

    Om Matlab inte verkar använda hela färgskalan, och istället verkar visa bilden i lite olika blåa nyanser, kan du istället pröva

    image(v,"CDataMapping","scaled")

  11. 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 1000⨯1000, 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. (Glöm inte skriva semikolon sist på raden, annars kommer 1000000 matrisvärden att skrivas ut.)

    Får du väldigt mycket utskrifter?

    Om ditt program fastnar i en oändlig loop, eller om du råkat få alla värdena utskrivna, kan du bryta funktionen genom att klicka i Matlabs kommandofönster ("Command Window") och trycka Ctrl-C.

    Får du fortfarande väldigt mycket utskrifter?

    Får du mängder av utskrifter trots att du avslutat raden med semikolon? Då kan det vara så att någon av raderna i din converge-funktion saknar semikolon. Öppna funktionen och se om ett väl placerat semikolon kan stoppa utskrifterna.

    Visa matrisen v som en bild, precis som i uppgift D10. Runt cirkeln finns en "aura" av lite avvikande färg. Vad beror den på?

  12. Nu ska vi äntligen undersöka Mandelbrotmängden. Denna mängd definieras av de komplexa tal för vilka följande serie konvergerar:

    z0
    c
    zk+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.

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

    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 Matlabs 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.7i, -0.5+0.6i);
    m = complexmat(1000, -1.4+0.48i, -1.1+0.24i);
    m = complexmat(1000, -1.8+0.03i, -1.5-0.03i);
    m = complexmat(1000, -0.18-1.02i, -0.14-1.05i);

  13. När du kommit hit är du klar med kursens andra Matlab-laboration. Läs igenom uppgifterna igen och se till att du förstår dina lösningar inför redovisningen.

    14. 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.

2021-04-20 10:29:53