Viktad anpassning till r¨ at linje - Mer om programmering och funktioner

6.6 Mer om programmering och funktioner

6.6.2 Viktad anpassning till r¨ at linje

Nästa uppgift blir att med utg˚angspunkt i samma schema designa och skriva en funktion som gör en viktad anpassning till en rät linje, y=a + bx.

Algoritm

Parametrarna a och b, samt os¨akerheten i dessa ges av:

a =

Pw_ix²_i ^Pw_iy_i−^Pw_ix_i ^Pw_ix_iy_i

∆ och b =

Pw_i ^Pw_ix_iy_i−^Pw_ix_i^Pw_iy_i

∆

σ_a= s

Pw_ix²_i

∆ σ_b =

sP w_i

∆ d¨ar wi = 1

σ_i² och ∆ =^Xwi

Xwix²_i −^Xwixi

62 KAPITEL 6. PROGRAMMERING Gr¨anssnitt och syntax

Indata till funktionen är lämpligen tre vektorer, en med x-värden, en med y-värden och en med felen i y-värdena. Vilka resultat vi vill f˚a ut ur funktionen kan variera en del. Vi kommer alltid att vilja se de anpassade värdena för parametrarna a och b, samt osäkerheten i dessa.

Eventuellt vill vi ocks˚a f˚a en uppfattning om kvaliteten i anpassningen, det vill säga chi-kvadratsumman per frihetsgrad och antalet frihetsgrader. Eftersom det, i motsats till lincorr ovan, inte är nödvändigt att räkna ut dessa senare tal för att funktionen skall kunna köra, s˚a lönar det sig att testa med nargout om dessa värden verkligen är efterfr˚agade eller inte innan de beräknas. Syntaxen för funktionen blir d˚a: function [parA, sigmaA, parB, sigmaB, ChiSquared, Dof] =linfit(x, y, sigma).

Felkontroll

Ett möjligt fel som m˚aste kollas är om de tre vektorerna som är argument till funktionen inte har samma längd. Vidare bör man kontrollera att ∆ inte är noll s˚a att divisionerna ger ett odefinierat resultat. Om chi-kvadrat skall beräknas m˚aste vi ocks˚a kontrollera att det finns minst tre mätta punkter, annars kommer divisionen när chi-kvadrat per frihetsgrad beräknas att bli odefinierad, samt att inget av mätfelen felaktigt är angivet som 0, i vilket fall bidraget till chi-kvadratsumman fr˚an den punkten blir odefinierat.

Testa funktionen

Du kan testa den f¨ardiga funktionen p˚a datam¨angden x =[1.097 1.700 3.39 1.399 3.70 2.190 2.5 4.00];

x=x*0.001;

y =[0.429 0.668 1.328 0.547 1.445 0.856 0.973 1.557];

dy =[0.001429 0.001668 0.0020 0.0015 0.0020 0.0010 0.0020 0.0030];

Svaret skall bli: a = (1.9347 ± 1.540)10⁻³, b = (389.99 ± 0.665), N = 6, χ²/d.o.f. = 2.38 6.6.3 Viktad anpassning och plottning av r¨at linje.

Här skall vi skriva en funktion som levererar en graf med data plottade med felstaplar och i samma graf den räta linjen som är ett resultat av anpassningen till en rät linje, med parametrarna för anpassningen utskrivna.

Algoritm

Anpassningen till en rät linje är inte n˚agot problem, vi l˚ater bara denna funktion i sin tur anropa linfit som vi skrev i förra avsnittet. Det gäller sedan ”bara” att f˚a grafen att se snygg ut utan att göra n˚agra antaganden i förväg om hur data ser ut. Det kan man

˚astadkomma genom att först plotta data, sedan flytta ut axlarna med 10% av det värde som den automatiska skalningen ger för att se till att alla punkter ligger väl inne i grafen.

D¨arefter drar vi den r¨ata linjen som definieras av anpassningen.

Om det anges i den kallande rutinen s˚a skall linplot skriva ut parametrarna fr˚an anpass-ningen. För att skriva ut parametrarna formaterar vi med formatet %0.3g för att det skall se snyggt ut oberoende av vilket resultatet av anpassningen blir. Vi vill sedan att texten skall hamna p˚a ett vettigt ställe och väljer därför att börja skriva 5% av grafens totala bredd in fr˚an vänsterkanten och översta raden 5% fr˚an toppen. Om riktningskoefficienten fr˚an anpassningen är negativ kommer datapunkter att ligga i det övre vänstra hörnet, vi väljer d˚a istället att i horisontell ledd börja att skriva i mitten av grafen. (Ett alternativ till denna algoritm kan vara att lägga till extra 20% i höjdled när vi skalar om vertikalt och använda

6.6. MER OM PROGRAMMERING OCH FUNKTIONER 63 denna area till att placera in texten.) Vi bör däremot undvika att använda kommandon som p˚averkar var grafen plottas, som figure(#) och subplot, det är bättre att det kallande programmet gör s˚adana val, eftersom önskem˚alen kan vara vitt skilda för olika tillämpningar.

Gr¨anssnitt och syntax

Gränssnittet för den här funktionen är detsamma som för linfit, de parametrar som skall skickas vidare till linfit, vektorerna x, y och sigma m˚aste ges som argument i anropet till linplot. En lyx-variant av funktionen skulle kunna ha som tillval att man även ger en sträng som anger hur grafen skall plottas av plot-kommandot, som till exempel ’-or’, men i första omg˚angen bryr vi oss inte om det. Däremot m˚aste vi skicka med tre textsträngar som ger text-strängar för titel, x-etikett och y-etikett. Om vi vill kunna välja om parametervärden skrivs ut i plotten eller inte s˚a m˚aste vi sätta in flaggor för detta i argumenten till funktionsanropen, parflag som är 1 om vi vill att parametervärdena skall skrivas in och 0 annars, chiflag styr p˚a samma sätt om chi-kvadrat och frihetsgrader skall skrivas ut (kom ih˚ag att variabler inte skall ha samma namn som funktioner, vi kan allts˚a inte använda title och xlabel som namn p˚a variablerna). Utdata fr˚an linplot är de samma som för linfit, allts˚a i första hand värdet p˚a parametrarna och osäkerheten i dessa, samt som tillval chi-kvadratsumman och antalet frihetsgrader. Vi kan dessutom tänka oss att linplot ibland kommer att kallas utan att den kallande rutinen använder utdata alls, i de fall det intressanta bara är att f˚a en graf. Syntaxen för funktionen blir d˚a : function [parA, sigmaA, parB, sigmaB, ChiSquared, Dof]

=linfit(x, y, sigma, titel, labelx, labely, parflag, chiflag)

Eftersom linplot skall kalla linfit m˚aste vi ocks˚a fundera över gränssnittet mellan dessa bägge funktioner. I princip är det ju inget som hindrar linplot fr˚an att kalla linfit med alla sex utdata parametrarna i anropet. Men om linplot inte är kallad med parametrarna ChiSquared och Dof, och om chiflag är noll s˚a behöver man ju inte beräkna chi-kvadratsumman i linfit. Att änd˚a kalla linfit med alla sex utdata parametrar betyder d˚a att programmet spenderar tid p˚a att beräkna saker som aldrig kommer att användas, vilket är slöseri. Vi m˚aste allts˚a i linplot testa om s˚a är fallet eller inte och bara anropa linfit med alla sex utdata parametrar om det är nödvändigt, annars använder vi [parA, sigmaA, parB, SigmaB]

Felkontroll

Det fel vi behöver testa för är om vektorerna x, y och sigma har samma längd, i annat fall kraschar errorbar-kommandot. Man kan tycka att vi inte behöver testa för det i linplot eftersom linfit testar för det. Men man skall försöka att skriva varje funktion s˚a komplett och sluten som möjligt. Det kan ju hända att vi i framtiden vill använda linplot s˚a att den kallas av n˚agon annan funktion än linfit, i en s˚adan situation kan det vara omöjligt att minnas att vi m˚aste kolla om den nya funktionen gör den nödvändiga testen.

64 KAPITEL 6. PROGRAMMERING Testa funktionen

Om funktionen fungerar b¨or kodsnutten nedan clear all;

x =[1.097 1.700 3.39 1.399 3.70 2.190 2.5 4.00];

x=x*0.001;

y =[0.429 0.668 1.328 0.547 1.445 0.856 0.973 1.557];

dy =[0.001429 0.001668 0.0020 0.0015 0.0020 0.0010 0.0020 0.0030];

dy=dy*10;

clf

[parA, sigmaA, parB, sigmaB, chi, N] = linplot . . . (x, y, dy,’Resistor 1’,’Str¨om (A)’,’Sp¨anning (V)’,1,1)

skapa en graf som ser ut som:

testa ocks˚a att texten dyker upp p˚a rätt ställe när riktningskoefficienten är negativ, och att vi kan sl˚a av och p˚a utskrift av parametervärden med flaggorna som det är tänkt.

6.7. OVNINGSUPPGIFTER¨ 65

6.7 Ovningsuppgifter ¨

1. Skapa en vektor med tio element som best˚ar av likformigt f¨ordelade slumptal (minns sektion 5.3.2) som antar heltalsv¨arden fr˚an 1 till 6.

2. Skriv ett program som med hj¨alp av en for-slinga genererar utskriften:

Hej - nu ar X:

1.200

Hej - nu ar X:

2.400

Hej - nu ar X:

3.600

3. Skriv en funktion som returnerar ett slumptal mellan 1 och 100. Talet skall vara ett heltal, dvs inga decimaler skall finnas med. F¨orsta raden i filen kan se ut s˚a h¨ar:

function y = slumptal()

4. Skriv en programsnutt som använder en while-slinga för att räkna hur m˚anga g˚anger du kallat slumptalsgeneratorn innan slumptalet blir större än 0.75.

5. Skriv ett program som genererar en slumpmässig tipsrad. Anta att sannolikheten för 1, x och 2 alla är 1/3, och l˚at programmet skriva ut tipsraden (en tipsrad best˚ar av tretton stycken tecken). Ett exempel p˚a en körning kan se ut s˚a här:

>> tipsrad 2

x 2 1 2 1 2 x 2 1 x 2 2

6. En metod f¨or att approximera kvadratroten av ett tal ¨ar Newtons metod som ger

√x ≈ y_k med y_k+1 = ¹₂y_k+_y^x

och y1 = 1. Skriv en funktion y=Newton(x) som använder denna metod för att beräkna kavdratroten ur ett godtyckligt positivt tal.

Kadratroten skall anges korrekt avrundad med tre decimalsiffrors noggrannhet.

Ledning: eftersom vi inte har n˚agot facit m˚aste vi i funktionen själva avgöra när serien har konvergerat tillräckligt för att vi skall ha den sökta noggrannheten. Detta kan vi t.ex. göra genom att kontrollera att |y_k− y_k−1| < 0.0001

Vi m˚aste dessutom hitta p˚a ett trick för att runda av korrekt, de rutiner vi har till hands rundar ju bara av heltal. Detta kan vi komma runt genom att transformera v˚art decimaltal till ett lämpligt heltal som vi sedan rundar av och därefter transformerar tillbaks till ett decimaltal.

66 KAPITEL 6. PROGRAMMERING 7. Skriv ett program som ber¨aknar en approximation till e, basen f¨or den naturliga

loga-ritmen, fr˚an summan.

e =

∞

k=0

1 k!

8. Skriv tv˚a funktionsfiler, medel och standavv, som beräknar medelvärde och stan-dardavvikelse för en vektor. Kontrollera att de ger samma resultat som mean och std!

Kapitel 7

N¨ ar det blir fel

Nästan inga program blir helt rätt fr˚an början. När det uppst˚ar fel i enklare program kan man ofta hitta orsaken bara genom att g˚a igenom programmet en g˚ang till med vetskap om vad som gick fel. Men i takt med att programmen blir större och mer komplicerade blir denna procedur mer och mer besvärlig och tidskrävande.

Ett fel i ett datorprogram kallas p˚a engelska ”bug”. Detta uttryck sägs härstamma fr˚an datorernas barndom. Vid ett tillfälle hade ett team som jobbade med en av de första stora beräkningsmaskinerna (rörbestyckad) länge försökt att ta reda p˚a vad sOctaveom gick fel i ett program. Till slut upptäckte man att orsaken inte l˚ag i programmet utan att en död insekt inuti maskinen hade orsakat en kortslutning. Denna historia m˚a vara sann eller inte, men ett fel i ett datorprogram heter p˚a engelska ”bug”, och processen att hitta och korrigera fel heter debugging. P˚a svenska säger man oftast bugg och debugga, men ibland kan man möta termen ”avlusa”.

I nästan alla programmeringsmiljöer finns det speciella verktyg att använda för att leta efter fel i sina program, s˚a kallade debuggers. Gemensamt för de flesta är att man med hjälp av verktyget kan köra sitt program rad för rad och vid valfri tidpunkt stanna upp och inspektera värdet p˚a de variabler som finns i programmet. Man kan ocks˚a köra programmet tills värdet p˚a en viss variabel ändras. Med hjälp av s˚adana verktyg kan man nästa alltid hitta ˚atminstone de fel i ett program som orsakas av att man tänkt fel när man skrivit programmet.

Att använda en debugger kräver en hel del träning för att det skall bli effektivt. I den här kursen ser vi detta som överkurs. I stället ger vi här n˚agra enkla tips om hur man med hjälp av de felmeddelanden vi kan se i kommandofönstret kan komma en bra bit p˚a vägen för att avhjälpa de vanligaste nybörjarfelen.

• A = Inf

A har beräknats p˚a ett sätt som “ger oändligheten”, ett fel som nästan alltid uppst˚ar för att man försöker dividera med noll.

• B = NaN

B ¨ar “not a number”, resultatet av en operation som inte ¨ar definierad, som t.ex. 0/0 eller cos(8/0),

• error: ’x’ undefined near line n column m error: called from:

error: Z:\slash...\slash funk.m at line n column m

Detta fel uppst˚ar när vi försöker använda en variabel som inte är initialiserad. En variabel m˚aste definieras antingen genom tilldelning (x = konstant) eller beräkning (x

= 3*y). Försöker vi använda en variabel innan den definierats f˚ar vi felmeddelandet 67

68 KAPITEL 7. N ÄR DET BLIR FEL ovan. Octave kan inte riktigt avgöra om vi försöker använda en icke-definierad variabel eller en icke-definierad funktion. Att det st˚ar “function” i felmeddelandet skall man allts˚a inte ta för bokstavligt.

• error: operator * : nonconformant arguments (op1 is nxm, op2 is mxn) error: called from:

error Z:\slash...\slash funk.m

Detta fel uppst˚ar när vi försöker multiplicera tv˚a matriser som inte har korrekta di-mensioner med varandra, dvs. antalet kolumner i den vänstra matrisen är inte lika med antalet rader i den högra.

• error: for A^b, A must be square

Igen ser vi i hyfsad klartext var felet ligger. H¨ar vill man g¨ora endera av tv˚a saker.

Antingen att ta skalärprodukten av x med sig själv, i vilket fall vi m˚aste komma ih˚ag att transponera en av vektorerna. Men ofta uppkommer det här felet när vi glömt att ta med punkten i den elementvisa operationen. Ofta är det istället x.∧2 vi vill beräkna.

• Syntax error parse error:

syntax error

>>> 3+/ 3

Syntax error kan vi f˚a p˚a m˚anga andra sätt, varje g˚ang vi inte följer Octaves spr˚akliga regler, dvs reglerna för hur vi skall formera v˚ara rader s˚a kommer vi att f˚a ett syntaxfel.

• error: invalid call to cos. Correct usage is:

-- Mapping Function: cos(X)

Detta felmeddelande är ocks˚a lätt att först˚a. Vi har förmodligen glömt att Octave använder decimalpunkt och inte decimalkomma, s˚a när vi skriver till exempel cos(2,4) tolkas det som tv˚a tal, och eftersom cosinusfunktionen bara accepterar ett argument s˚a f˚ar vi ett fel.

Kapitel 8

2D-Grafik

I det här passet kommer vi att koncentrera oss p˚a att skapa grafik: histogram, plottar och annat som kan användas för att ˚ask˚adliggöra experimentella resultat. Vi skall försöka avst˚a fr˚an frestelsen att komma med pekpinnar om hur grafik skall utformas för att vara bra, utan i stället koncentrera oss p˚a att g˚a igenom tekniken, s˚a att vi kan ˚astadkomma grafik som ser ut som vi vill.

Efter att ha g˚att igenom det här passet skall du själv kunna producera en plot som den ovan och ocks˚a kunna b˚ade skriva ut den och kopiera in den i ett dokument för att bifogas i en lab-rapport.

70 KAPITEL 8. 2D-GRAFIK

8.1 Kurvor i tv˚ a dimensioner

Att skapa enklare tv˚a-dimensionella plottar i Octave är relativt enkelt. En skillnad mellan Octave och en del andra grafikprogram är att vi inte kan instruera Octave att till exempel rita kurvan y=sin(x). Istället m˚aste vi ge Octave en uppsättning punkter med en x-koordinat och en y-koordinat. Octave ritar sedan upp en kurva som förbinder punkterna. Vill vi rita upp y = sin(x) m˚aste vi därför först själva skapa en x-vektor och sen en y-vektor¹ där y=sin(x) och instruera Octave att rita en kurva mellan de punkter som beskrivs av dessa vektorerna. L˚at oss börja med n˚agot enkelt som att rita grafen för y = 2x + 3.

>x = [1 8];

> y = 2*x + 3;

> plot (x,y)

8.1.1 Att växla mellan fönster, och städa i fönster.

När vi trycker p˚a return efter att ha givit plot kommandot öppnas ett nytt fönster, Octaves grafikfönster, där plotten visas. När vi jobbar med grafik kommer vi att vilja pendla mellan m˚anga fönster. Ofta kommer vi att utföra samma lunk mellan editorns fönster där vi just

andrat n˚agon detalj i en m-fil till kommandofönstret för att köra filen till grafikfönstret för att se hur det blev. Det finns flera metoder för att snabbt växla mellan fönster, vilken som passar bäst just för stunden beror litet p˚a exakt vad man h˚aller p˚a med och hur man föredrar att jobba.

• Det kanske enklaste sättet är att helt enkelt lägga fönstren bredvid varandra, eller i alla fall se till inget fönster ligger helt bakom de andra. Genom att klicka och dra p˚a kanten av ett fönster kan du flytta det vart du vill p˚a datorns skrivbord, och genom att klicka och dra i ett aktivt fönsters nedre högra hörn kan du ändra dess storlek. Om du lägger upp de fönster du jobbar i s˚a att det alltid sticker fram n˚agon bit av fönstret

aven när det är de andra fönstren som är aktiva s˚a kan du enkelt byta fönster bara genom att klicka med musen i det fönster du vill skall bli det aktiva.

• I windows finns det en funktion för att bläddra mellan alla fönster som är öppna p˚a datorns skrivbord. H˚all ”alt”-tangenten nedtryckt och tryck p˚a TAB tangenten. Varje g˚ang du trycker ned TAB bläddras ett nytt fönster fram som aktivt, vilket fönster som görs aktivt visas i en liten ruta.

Vill du städa upp i grafikfönstret skriver du (i kommandofönstret) clf (CLear Figure window), vill du städa upp i kommandofönstret ger du kommandot clc (CLear Command window).

Du kan st¨anga grafikf¨onstret genom att ge kommandot close.

8.1.2 plot-kommandot

L˚at oss prova att plotta en litet mer komplicerad kurva, y = x². Pr¨ova

> clear all;

> x = [1 4];

> y = x.∧2;

> plot(x,y)

Resultatet blir inte direkt vad vi väntat oss, men egentligen är det inte s˚a märkligt: plot kommandot ser ju bara de punkter vi vill rita en graf mellan och har ingen aning om huru-vida de tv˚a punkterna (1,1) och (4,16) är tv˚a punkter p˚a kurvan y=x², eller om de ligger p˚a

1Eftersom det Octave anv¨ander i plottandet ¨ar punkter, vars x- och y-koordinater kommer fr˚an vektorerna x och y s˚a m˚aste vektorerna ha samma antal element.

8.1. KURVOR I TV˚A DIMENSIONER 71 kurvan y = 5x-4. Mer information m˚aste allts˚a till, pr¨ova att l¨agga till punkten ( x,y)

= (2,4). Inte heller särskilt imponerande. För att kurvan skall bli jämn och fin och följa det funktionssamband vi önskar m˚aste vi ge Octaves plot-kommando vektorer där punkterna ligger ganska tätt. Det görs enklast om vi använder oss av Octaves vektorfunktioner. Antag till exempel att vi vill plotta y=x²mellan 0 och 4.5. Vi kan d˚a antingen använda kolonoper-atorn för att skaffa oss en x-vektor där elementen har värden fr˚an 0 till 4.5 med vilket avst˚and vi vill, till exempel 0.1: x = 0:0.1:4.5; (när vektorerna f˚ar s˚a här m˚anga element vill man gärna stänga av ekot genom att ge semikolon efter kommandot). Alternativt kan vi använda linspace för att skapa 50 punkter mellan 0 och 4.5: x = linspace(0,4.5,50); ger oss ocks˚a en användbar x-vektor. När x-vektorn är klar skapar vi en y-vektor med y = x.∧2 och plottar resultatet.

Prova nu att plotta n˚agra andra funktioner som polynom och sin(x) för att se att resultatet blir som du förväntar dig.

8.1.3 Styra utseendet p˚a grafiken

Det mesta av den tid man spenderar p˚a att göra grafer tenderar att läggas p˚a att f˚a lay-outen att bli vad man har tänkt sig. Eftersom korrekt och tydlig lay-out är ett viktigt verktyg för att ˚astadkomma grafik av hög kvalitet kan det vara värt besväret att redan fr˚an början försöka lära sig de kommandon som krävs för att man skall kunna styra programmet dit man vill.

Linjef¨arg, linjetyp och mark¨orer

Kommandot plot kan ges med tre argument: plot(x,y,S). Argumentet S är en textsträng som anger hur linjen som ritas skall formateras. Det finns tre typer av information att styra, färgen p˚a linjen, hur linjen dras (heldragen, prickad mm), och slutligen markör, eller vilken symbol som skall ritas ut i de punkter som anges av vektorerna. Skriver du help plot kan du bland annat hitta följande lista p˚a vilka möjligheterna är:

Linjef¨arg Mark¨or

1/r red * star

2/g green s square 3/b blue v triangle 4/m magenta o circle

5/c cyan + plus

y yellow . dot

k black p pentagram

Man kan i plot kommandot använda en godtycklig kombination av dessa (men bara ett val för varje kategori!). Prova med litet olika kombinationer för att se hur det fungerar, en gul kurva med heldragen linje mellan sm˚a cirklar ges t.ex av plot(x,y,’yo-’). Märk att om man ger formatkommandot för linjefärg och markör men inte n˚agon kod för linjetyp som

In document Octave i kursen experimentella metoder (Page 67-0)