Matlab-laboration 1
EDAA55, LTH

Mål

Efter laborationen ska du

Laborationen avslutas med en temauppgift, där du får träna på att konstruera script i Matlab.

    Grunder: komplexa tal och enkla script

  1. Ubuntu: ikonen 'Search your computer' Logga in på datorn och starta Matlab. På LTH:s Linux-datorer klickar du på Ubuntu-ikonen, som finns överst till vänster på skärmen (och ser ut som bilden här till höger). Skriv matlab och tryck Return.

    På samma Linux-datorer kan du enkelt få dina två fönster (Matlab + denna webbläsare) sida vid sida. Klicka och dra i Matlab-fönstrets titelrad, och dra fönstret till skärmens kant (vänster eller höger). När du ser en "skugga" över halva skärmen släpper du fönstret där. Gör därefter likadant för webbläsaren, fast dra den istället till den andra kanten.

    Gör några enkla beräkningar (de vanliga räknesätten) i kommandofönstret (Command Window) och prova att använda variabler.

  2. Med Matlab kan man enkelt hantera såväl reella som komplexa tal.

    Skriv ut följande värden:

    1. π
    2. √3
    3. √–3
    4. √i
    5. eπi

  3. Komplexa tal presenteras alltså i rektangulär form (det vill säga med real- och imaginärdel). Ibland är man istället intresserad av talets polära form (absolutbelopp och argument).

    Tilldela variabeln z något av de komplexa värdena i uppgift D2(c)-(e) ovan. Skriv ut absolutbelopp och argument för z.

    Varför funkar inte arg(z)?

    Funktionen angle(z) ger argumentet (vinkeln) för ett komplext tal z.

    Tilldela z ett annat av de tre komplexa värdena och skriv ut real- och imaginärdel för z.

  4. Ofta är det praktiskt att spara en sekvens av kommandon i en fil. Då har man kommandona sparade till senare, och kan man köra dem tillsammans igen. Sådana script är alltså ett slags små program.

    Bild: knappen New/Script i Matlab

    Klicka på New→Script (se bilden ovan). Ett nytt delfönster dyker upp med ett script som än så länge heter "untitled". Skriv in följande text, klicka på Save-knappen och spara scriptet som lab1.m. Låt scriptet förbli öppet – du kommer att bygga vidare på det.

    %% Komplexa tal   a = 3 + 2i; b = 5 - i; x = a * b

    Den första raden är en rubrik. I nästa uppgift kommer du att få se hur sådana rubriker används.

    Notera att scriptet lab1.m nu återfinns bland dina Matlab-filer i delfönstret Current Folder. Därifrån kan du öppna det senare om du vill. Det finns kvar även efter du stängt Matlab.

    Klicka på Run. Alla tre satserna i ditt script körs, men du kan bara se resultatet från en av dem. Vilken? Varför?

    Du kan även köra samma script genom att ange dess namn i kommandofönstret. Skriv lab1 i kommandofönstret och kontrollera att du får samma svar.

    Alla tre variablerna är definierade i Matlab, även om bara en är utskriven. Pröva gärna följande i kommandofönstret:

    b 5.0000 - 1.0000i

    Här har det du ska skriva in markerats i fetstil (b). Den övriga texten är det förväntade resultatet från Matlab.

    5. Funktioner

  5. Inför de följande uppgifterna ska du skapa ett nytt avsnitt i ditt script.

    Lägg till en rad med två procenttecken och en rubrik sist i ditt script, exempelvis så här:

    %% Funktioner

    En sådan rad som inleds med två procenttecken (%%) markerar början på ett nytt avsnitt i scriptet. Detta avsnitt kan köras separat, utan att man behöver köra hela scriptet på en gång. Då använder man knappen Run Section, som du finner strax till höger om Run-knappen (se bilden nedan).

    Bild: knappen Run Section i Matlab

    Placera textmarkören i scriptets första avsnitt (Komplexa tal) och tryck på Run Section. Kontrollera att samma sak som i förra uppgiften sker.

    Placera därefter markören i scriptets andra avsnitt (Funktioner) och tryck på Run Section. Nu händer ingenting, eftersom detta avsnitt än så länge är tomt.

  6. Låt oss definiera funktionen gx som följer:

    gxe–0.1x cos x

    Definiera gx i ditt script. Gör definitionen på en enda rad med en s.k. anonym funktion, så som beskrivs i Matlab Primer.

    Glöm inte köra script-avsnittet med funktionsdefinitionen.

    Kontrollera (i kommandofönstret) att funktionen ger rätt resultat:

    g(0) 1 g(2) -0.3407

    Om något inte stämmer, justera ditt script och tryck Run Section igen.

  7. Vi är nu intresserade av funktionens värden i intervallet 0 ≤ x ≤ 10.

    Lägg till en rad i ditt script, där en vektor av 100 sådana x-värden skapas, jämnt fördelade.

    Kontrollera (i kommandofönstret) att värdena verkar rimliga:

    x(1) 0 x(27) 2.6263 x(100) 10

    Undersök Matlabs Workspace-delfönster. Om du inte ser något sådant delfönster kan du öppna det genom att klicka på knappen Layout och välja Workspace (under rubriken Show).

    I Workspace-fönstret ser du vilka variabler Matlab känner till just nu, och vilka värden de har.

    Vilken information om x kan du få ur Workspace-fönstret? Är x en rad- eller kolonnvektor? Hur många element har den? Kan du se funktionen du definierade ovan?

    Nu vill vi beräkna värdet av funktionen g för alla dessa x-värden.

    Lägg till en rad i ditt script, så att en motsvarande vektor y av funktionsvärden skapas.

    I Matlab Primer finns ett exempel där en serie funktionsvärden beräknas. Kan du komma på hur du kan använda samma teknik här?

    y = g(x);
    Vad är 'Incorrect dimensions for matrix multiplication'?

    Du kan behöva justera din funktion så att den använder elementvisa operationer.

    Kontrollera att även dessa värden verkar stämma:

    y(27) -0.6692 g(x(27)) -0.6692

    I det senare kommandot är x(27) värdet av element 27 i vektorn x, och det värdet används i sin tur som parameter till funktionen g. Resultatet ska överensstämma med element 27 i vektorn y.

    Låt ditt script plotta funktionen (gx enligt ovan) för värdena 0 ≤ x ≤ 10.

    Ditt Funktioner-avsnitt ska nu se ut i stil med följande (fast du ersatt ... med något annat):

    %% Funktioner   g = ...; x = ...; y = ...; pl...

    I fortsättningen kan du själv välja när du vill skapa nya avsnitt i ditt script. Du kan också själv välja vilka kommandon du skriver i ditt script, och vilka du skriver direkt i kommandofönstret. Fråga gärna handledaren om du är osäker.

    8. Något om syntax i Matlab och Java

    Matlab och Java är två olika programspråk för lite olika ändamål: Matlab är inriktat på små program av numerisk/matematisk karaktär, och Java är inriktat på större program och generellare tillämpningar. De två språken har både likheter och skillnader.

    Matlabs funktioner motsvaras i Java närmast av metoder, även om det finns viktiga skillnader. I Java finns även vektorer, som liknar Matlabs vektorer på många sätt. (Om du läser F, I eller Pi har du nog ännu inte kommit i kontakt med Java-vektorer, men det kommer du snart att göra.)

    Matlab använder samma syntax till både funktioner och vektorer, trots vi här ser dem som två ganska olika saker. I förra uppgiften skrev vi g(0), och då är g en funktion och 0 ett parametervärde. Här skriver vi x(1), där x är en vektor och 1 ett index. Vad ett uttryck som a(3) betyder beror alltså på om a är en funktion eller en vektor.

    I Java skrivs dessa två saker olika: ett metodanrop inom klassen kan exempelvis skrivas g(0), medan en vektor indexeras med klamrar, exempelvis x[1]. I Java indexeras vektorer dessutom från 0, medan Matlab indexerar dem från 1.

  8. Ibland vill man (approximativt) beräkna derivatan för en given funktion. Då kan man utgå från derivatans definition:

    f'x≈
    fx + h – fx
    h
    (om h är litet)

    I analysen låter man ju h gå mot noll. Här måste vi bestämma oss för ett värde, och låter därför istället h vara ett "lagom" litet tal. Vi kan exempelvis sätta h till 10–6:

    h = 1e-6;

    Definiera en funktion deriv som tar två parametrar: en funktion f, och ett x-värde. Funktionen deriv(f,x) ska beräkna (approximera) f'x enligt ovan.

    Du definierar funktionen med något i stil med

    deriv = @(f,x) ...

    där den saknade delen (...) hämtas från derivatans definition ovan. Parametern f är en funktion, och x är ett x-värde vi vill beräkna derivatan för.

    Derivatans definition kan skrivas i Matlab som

    deriv = @(f,x) (f(x+h) - f(x)) / h;

    Använd en funktion med känd derivata för att testa din deriv-funktion för något x-värde. Hur nära kommer approximationen?

    Vi vet ju exempelvis att derivatan av sin x är cos x. När du anger sin som parameter i Matlab måste du skriva @sin (se Matlab Primer, "Function Handles").

    Vad ger följande kommandon för resultat? Är resultaten rimliga?

    deriv(@sin,0.3) cos(0.3) deriv(@sin,0.3) - cos(0.3)

  9. Ett annat sätt att övertyga sig om att deriv ger rimliga resultat är att plotta en funktion tillsammans med sin (approximerade) derivata. Här kan vi göra detta för funktionen g(x) i uppgift D6.

    gxe–0.1x cos x

    Plotta gx och g'x i samma diagram, för 0 ≤ x ≤ 10. Använd deriv för att approximera g'x, och gör grafen för g'x röd.

    Tips om hur du väljer färg hittar du i Matlab Primer, "Specifying Line Styles and Colors".

    Hur får man två grafer i samma bild?

    För att få flera grafer i samma bild kan man skriva hold on efter man gjort den första grafen. När man sedan gör nästa graf ritas den tillsammans med den första, istället för att ersätta den.

    Verkar den approximerade derivatan rimlig?

    Derivatan ska exempelvis vara noll vid funktionens lokala maxima/minima. Derivatan ska i sin tur ha ett lokalt maximum då funktionen växer som brantast.

    För att se axlarna tydligare kan du använda axis on.

    Behåll fönstret med graferna öppet. Vi ska nu se hur man kan utrusta det med rubriker och spara bilden som en fil.

    Lägg till följande rader i ditt script, och notera hur fönstret med graferna förändras:

    title('Funktionen och dess derivata') legend('funktion', 'derivata') print('funktionsbild','-dpng')

    Det sista kommandot ovan skapade en fil som heter funktionsbild.png och hamnat bland dina Matlab-filer (se i delfönstret Current Folder). Sådana png-bilder kan importeras i ett ordbehandlingsprogram (LaTeX, Word eller något annat) och användas i exempelvis en labbrapport.

  10. Bent, Alva och Kit försöker derivera funktionen gx ovan analytiskt. Det är länge sedan de läste analys, och endast ett vagt minne av deriveringsreglerna återstår. De kommer fram till olika svar.

    h1x
    e–0.1x cos x – 0.1e–0.1x sin x
    (Bent)
    h2x
    –0.1e–0.1x cos x – e–0.1x sin x
    (Alva)
    h3x
    –0.1e–0.1x cos x – 0.1e–0.1x sin x
    (Kit)

    Vi ska nu använda Matlab för att försöka jämföra dessa tre kurvor h1x, h2x och h3x med din beräknade derivata deriv(g,x). Du ska alltså inte använda dina kunskaper om deriveringsreglerna för den här uppgiften.

    Skapa fyra vektorer av derivatavärden: en från deriv(g,x), en från h1x, en från h2x och en från h3x.

    y0 = deriv(g,x); y1 = h1(x); y2 = ...; % på motsvarande sätt y3 = ...; % på motsvarande sätt

    Skapa tre nya vektorer, en för var och en av de tre funktionerna. Varje vektor ska innehålla differenserna mellan funktionens värden och den beräknade derivatans.

    Varje vektor anger alltså hur mycket respektive funktion avviker från den "riktiga" (beräknade) derivatan.

    Med beteckningarna från föregående ledtråd blir

    d1 = y1 - y0; d2 = ...; % på motsvarande sätt d3 = ...; % på motsvarande sätt

    Vi har nu tre differensvektorer, en för varje funktion. De är rätt jobbigt att läsa igenom och jämföra tre vektorer om hundra element. Istället vill vi ha ett enda mått på respektive differensvektors storlek.

    Här är funktionen norm, som beräknar en vektors magnitud (euklidisk norm), användbar. Den är ett enkelt och bra storleksmått på en vektor. Läs gärna om den euklidiska normen på Wikipedia, eller pröva help norm.

    Vilken av de tre funktionerna kan anses ligga närmst den beräknade derivatan?

    Med beteckningarna från föregående ledtråd kan man göra

    norm(d1) norm(d2) norm(d3)
    Ett litet resultat from norm betyder att vektorn innehåller små värden. En av normerna ligger nära noll, och motsvarar alltså den korrekta derivatan.

    11. Matriser och ekvationssystem

  11. Vi söker en matris med följande utseende:

    M

    Bilda matrisen M ovan med hjälp av Matlab-funktionerna eye och ones.

    Ovan till höger finns en knapp []. Om du trycker på den så genereras en ny uppgift. Pröva!

  12. Betrakta följande ekvationssystem (med tre obekanta , och ):

    Lös ekvationssystemet genom att skriva om det på matrisform och använda operatorn "\".

    Matrisform? Vad är det?

    Ekvationssystem på matrisform beskrivs i Matlab Primer, avsnitt 3.

    Ekvationssystemet kan skrivas A·xb, där

    A
               x
               b

    Med dessa beteckningar söker vi kolonnvektorn x.

    Kontrollera att din lösning uppfyller ekvationerna. Även det kan göras enkelt i Matlab — hur då?

    På matrisform skriver vi ju A·xb, där x är lösningen du fann ovan. Med dessa beteckningar ska både vänsterledet A*x och högerledet b ska ha samma värde.

    Kontrollera detta genom att skriva ut dessa båda värden i Matlab.

    x
      =  

  13. Betrakta matrisen A nedan:

    A
    11–1
    211
    43–1

    Bestäm determinanten för A. Använd en av Matlabs inbyggda funktioner.

    Determinant? Vad är det?

    Determinanter behandlas i kursen i linjär algebra, som du kanske ännu inte läst. Det är inget bekymmer här.

    Determinanten är ett värde som man kan beräkna för en given matris, och som säger oss något om matrisens egenskaper. Exempelvis är det intressant, visar det sig, att veta huruvida matrisens determinant är noll eller ej.

    Vilken inbyggd funktion ska jag använda?

    En bra gissning är att funktionens namn påminner om "determinant". Pröva att söka efter "determinant" i Matlab Primer.

    Om matrisen heter exempelvis A kan man skriva

    det(A)

    Använd Matlab för att lösa följande ekvationssystem:

    x1 + x2 – x3
    2
    2x1 + x2 + x3
    3
    4x1 + 3x2 – x3
    4

    Du får sannolikt fram ett märkligt felmeddelande, och värdena i resultatet ser nog också konstiga ut. Hur ska man tolka detta? Jämför med matrisen A ovan.

    I kursboken i linjär algebra kan man läsa följande:

    Huvudsatsen. För kvadratiska matriser A är följande villkor ekvivalenta:

     
    [...]
    (c)
    ekvationssystemet AX=Y är lösbart för alla Y,
    (d)
    A är inverterbar,
     
    [...]
    (f)
    det A ≠ 0 .

    Ur: Gunnar Sparr, Linjär algebra, andra upplagan, Studentlitteratur 1995, s. 212.

    Ekvationssystemet saknar lösning (vilket också indikeras av att determinanten är noll).

    14. Tema: från kilskrift till Matlab-script

  14. Du har redan sett att Matlab kan användas för att beräkna kvadratrötter, till och med komplexa sådana. Och i Java kan vi naturligtvis använda Math.sqrt(x) för att få fram roten ur reella x. Här ska vi titta närmare på hur en sådan beräkning kan fungera.

    Kvadratroten av ett positivt, reellt tal S kan beräknas iterativt så här (mer om algoritmen nedan):

    xk+1
    xk + ( S / xk )
    2

    Här beräknas alltså successivt allt bättre approximationer xk+1 till √S, hela tiden utifrån den föregående approximationen xk. Med "tillräckligt" många upprepningar får vi en approximation av √S. Som startvärde kan vi exempelvis gissa

    x0S /  2
    Låt oss beräkna √50 med sex decimaler. Då blir

    x0 S / 225
    x1
    x0 + ( S / x0 )
    2
    25 + (50 / 25)
    2
    13.5
    x2
    x1 + ( S / x1 )
    2
    13.5 + (50 / 13.5)
    2
    8.601851
    x3
    7.207277

    x4
    7.072355

    x5
    7.071068

    x6
    7.071068


    I det här fallet får vi alltså sex korrekta decimaler efter fem iterationer. Den sjätte iterationen bekräftar att vi är klara.

    Skriv ett Matlab-script som läser in ett tal S från tangentbordet, och därefter skriver ut de tio första approximationerna av √S med algoritmen ovan.

    Använd gärna knappen New→Script för att komma igång.

    Ditt script ska inte använda sig av Matlabs sqrt-funktion.

    Kontrollera att ditt script fungerar, exempelvis genom att jämföra med exemplet ovan.

  15. Som du kanske ser konvergerar beräkningen ganska fort. Det är sällan nödvändigt med tio iterationer, men det är svårt att på förhand avgöra exakt hur många som behövs.

    Ändra ditt script så att det inte gör ett fixt antal iterationer (tio ovan), utan beräknar nya värden tills serien konvergerar. Låt scriptet skriva ut antalet iterationer som behövdes.

    Använd följande kriterium för konvergens:

    xk – xk–1 < 10–6

    Tips: Om ditt program fastnar i en oändlig loop kan du bryta det genom att klicka i kommandofönstret och trycka Ctrl-C.

    Du behöver två x-värden för iterationen: en variabel för nuvarande x-värde, och en för det föregående. Välj deras startvärden på ett sådant sätt att while-satsen inte avbryts från början.

    Du använde Matlabs funktion för absolutbelopp i en av laborationens första uppgifter.

    Om ditt script inte verkar konvergera, utan fastnar i en oändlig loop, fundera då på ditt while-villkor.

    • Lägg in en ny rad i din loop, där båda variablerna skrivs ut (om variablerna exempelvis heter x1 och x2):

      fprintf('%.6f %.6f\n', x1, x2)

      Här ska de första raderna vara olika, men efter några iterationer ska de två värdena bli lika. Om värdena är helt lika redan från början så behöver du titta närmare på i vilken ordning dina variabler får sina värden. Hur ska du göra för att komma ihåg både det nuvarande och det förra värdet?

    • Använder du < eller > i ditt while-villkor? Vilket ska villkoret vara för att loopen ska fortsätta?

    Tips: När man skriver Matlab-script kan det vara praktiskt att inleda med en serie kommandon som återställer Matlab till ett känt ursprungsläge. Exempelvis kan man inleda sitt script med följande kommandon:

    clear        % glöm alla tidigare variabelvärden close all    % stäng alla figurer clc          % rensa kommandofönstret

    Pröva gärna att stoppa in kommandona ovan i ditt kvadratrotsscript och se att de fungerar som du förväntar dig.

  16. När du kommit hit är du klar med den första Matlab-laborationen. Läs igenom uppgifterna igen och se till att du förstår dina lösningar inför redovisningen.

    17. Bakgrund: mer om algoritmen

    YBC 7289
    Lerskivan YBC 7289 (Yale Babylonian Collection; bild: Bill Casselman; licens: CC BY 2.5).
    På diagonalen står 1·600 + 24·60–1 + 51·60–2 + 10·60–3 ≈ 1.414213.
    Tecknen är jämförelsevis stora (ca 8mm), vilket ofta tolkats som att de skrivits av en student.

    Approximationstekniken ovan användes av babylonierna (i nuvarande Irak) kring 1600–1800 f.Kr., även om de naturligtvis inte använde samma symboler som oss. De använde ett positionsbaserat talsystem baserat på basen 60, och alltsedan dess delar vi upp timmar i 60 minuter, och minuter i 60 sekunder. Och de hade alltså även en metod för att beräkna kvadratrötter.

    Det är inte helt klart om babylonierna gjorde mer än en iteration i sina kvadrotsberäkningar, men vi vet att de kände till approximationen ovan. De kände också till √2 med en noggrannhet motsvarande ungefär sex decimaler. I alla händelser beskrev Heron (ca 10–70 e.Kr.), verksam i Alexandria i nuvarande Egypten, den iterativa algoritmen ovan. Idag kan vi även komma fram till samma algoritm om vi utgår från Newton-Raphsons metod. Då ser vi problemet som att söka nollställen till funktionen

    fxx2 – S

    Persondatorer har idag ofta en särskild instruktion för kvadratrotsberäkning, men även sådan hårdvara använder en iterativ algoritm av ovanstående typ. I många tekniska tillämpningar använder man dessutom enklare processorer, som saknar sådana instruktioner. Då kan man beräkna kvadratrötter på det sätt som du prövat här, precis som babylonierna gjorde för närmare fyratusen år sedan.

    Mer läsning (för den intresserade)

    David Fowler och Eleanor Robson, "Square Root Approximations in Old Babylonian Mathematics: YBC 7289 in Context": Historia Mathematica 25 (4), 1998, s. 366–378.
    Bo Göran Johansson, Matematikens historia, andra upplagan, Studentlitteratur 2013, avsnitt 1.6 och 3.6. (Finns på UB.)
    Jan Thompson, Matematiken i historien, Studentlitteratur 1996, kap. 3. (Finns på UB.)
2021-11-17 08:59:06