opgaver:Uge1

For en hjælp til hvordan man plotter data, se "MATLAB® Primer" (tidligere kaldet "Getting Started with MATLAB") guiden, der ligger på Absalon. Kapitel 4 indeholder en masse plotning generelt. Prøv at slå funktionen errorbar op i MATLAB hjælpen bagefter.

Plottet som det burde se ud efter at have kørt denne kode.

For at plotte datapunkterne med usikkerheder skal data først indlæses. Usikkerhederne er kvadratroden af tælletallet, så

t = 0:0.5:2;

nu = [214 134 101 61 54];

nu_usikk = sqrt(nu);

Herefter kan data plottes vha. MATLAB funktionen errorbar(x,y,e), der plotter $(x,y)$ data, hvor der er usikkerheden $e$ på $y$-dataen. En figur åbnes med

figure

der plottes med

errorbar(t,nu,nu_usikk,'ro')

og til sidst laves der en

hold on

for at blive ved med at plotte oven i det allerede lavede (ellers vil de næste kommandoer overskrive errorbar plottet). Jeg sætter nu også akser, titel og akselabels med

axis([-0.5 2.5 40 240])

title('Opgave 1 - Radioaktivt henfald')

xlabel('Tid i timer')

ylabel('Tælletal')

For at plotte en funktion skal man først definere sin førsteakse som en vektor, f.eks. x = 0:0.01:1. Man kan også bruge linspace(x,y,n) funktionen til at lave en liste fra $x$ til $y$ med $n$ punkter. Antallet af punkter definerer opløsningen på det plot man vil lave.

Herefter kan man udregne sine værdier til andenaksen som sin funktion af førsteakseværdierne, f.eks. y = sin(x).

Det endelige plot, der sammenligner data med usikkerheder med en teoretisk linje.

Man definerer linjen for teorien rent numerisk ved at lave en tids-akse som man selv bestemmer hvor detaljeret er, f.eks.

tt = -0.25:0.01:2.25;

og så kan man udregne den teoretiske linje ved

nu_0 = 200;

r = 0.693;

nu_exp = nu_0 * exp(-r*tt);

Hvis man huskede sin hold on tidligere burde man nu bare kunne plotte oven i sit errorbar plot med

plot(tt,nu_exp,'b:')

Som en sidste krone på værket kan man tilføje en legend med

legend('Datapunkter','Teori-linje')

Prøv at lave plottet igen med forskellige layouts af punkterne og linjen. F.eks. er der de forskellige farver r (rød), g (grøn), b (blå), y (gul), m (magenta) og k (sort); punkttyperne . (lille prik), x (krydser), o (store prikker) og * (stjerner); og linjetyperne - (solid), : (prikker) og -- (stiplet), og selvfølgelig mange flere. Prøv at finde disse options i MATLAB hjælpen under LineSpec.

Du kan også lege med egenskaber som MarkerSize, LineWidth, MarkerEdgeColor og MarkerFaceColor for at gøre linjen og punkterne mere tydelige.

Plottet som før, men nu også med en fit-linje.

Lad os sige at data allerede er indlæst i variablene t og nu (med usikkerhed nu_usikk), som ovenfor. Disse er også allerede plottet som det blev gjort i spørgsmål a.

For at kunne fitte en linje skal man selvfølgelig først definere sin fitte-kurve, der i dette tilfælde vil være $\nu = \nu_0 \exp{(-rt)}$, hvor $\nu_0$ og $r$ nu er fitte-parametre. En funktion kan defineres på mange måder i MATLAB, men når man fitter med simple funktionsudtryk i MATLAB er det nemmest simpelthen bare at fodre fitte-rutinen med den direkte.

cfun = fit(t',nu','a*exp(-b*x)')

Bemærk at t' betyder den transponerede af t. Funktionen fit kræver nemlig søjlevektorer som data, ikke rækkevektorer, som både t og nu er defineret som. Objektet cfun indeholder nu information om fittet, og man kan plotte kurven med de bedste værdier af a og b med

plot(tt,cfun(tt),'g-')

man kan også få parametrene ud ved at kalde variablene cfun.a og cfun.b. Det er også muligt at gå dybere ned i cfun objektet, men dette skal ikke uddybes her.

Bemærk også at ovenstående ikke fodrer MATLAB med et start-gæt på de to fitte-parametre. Dette bør man normalt gøre, da programmet ellers vælger tilfældige startpunkter, og her ikke altid vil konvergere. Man kan gøre det ved at bruge StartPoint option'en, og dette vil vi komme ind på senere.

Se på hvad antagelserne for at formel (4.14) gælder faktisk er - passer dette med hvad der arbejdes med i denne opgave?

For at se på det i MATLAB kan man generere en række tilfældige tal vha. randn funktionen. Slå den op i MATLAB hjælpen og se hvad den kan.

Konklusionen i den fejlagtige udledning er forkert fordi man har brugt udtrykket i formel (4.14) uden at overholde dennes antagelser - nemlig at de variable man ganger sammen er uafhængige. Og nu er $x$ jo fuldkommen afhængig af sig selv - hvorfor konklusionen altså kun kan være forkert.

For at se på dette i MATLAB kan man starte med at lave en liste med normalfordelte tal

x = 10 + randn(1,500);

Bemærk at der lægges 10 til, da vi gerne vil dividere med gennemsnittet senere (og at dividere med nul er jo ikke så godt!).

Man kan nu helt manuelt se på den relative usikkerhed på x og på x i anden potens, ved at bruge de indbyggede MATLAB funktioner for standardafvigelse og gennemsnit, std og mean:

std(x)/mean(x)

std(x.^2)/mean(x.^2)

Man finder at nummer to usikkerhed ovenfor ligner at man ganger den første med 2, ikke med $\sqrt{2}$, ved at sammenligne

2*std(x)/mean(x)

sqrt(2)*std(x)/mean(x)

Og dermed er det altså også rent numerisk vist at den naive fremgangsmåde ikke er sand.

Først defineres datasættet som en enkel liste med

nu = [10 13 8 15 8 13 14 13 19 8 13 13 7 8 6 8 11 12 8 7];

Derefter kan gennemsnittet og SD udregnes enten manuelt eller ved brug af mean funktionen. F.eks.

gn = sum(nu)/length(nu)

std(nu)

Dermed har man gennemsnittet og datasættets faktiske standardafvigelse.

Husk ikke at forvirre standardafvigelsen SD med standardafvigelsen på gennemsnittet (også kaldt SDOM, for Standard Deviation Of the Mean). Hvis man her vil finde SDOM, skal man dividere med kvadratroden af antallet af punkter,

std(nu)/sqrt(length(nu))

Kvadratroden af gennemsnittet udregnes simpelthen med (ved brug af løsningen fra spørgsmål a)

sqrt(gn)

Denne standardafvigelse er faktisk et over-estimat, men den passer ok med den udregnede std(nu) fra ovenfor.

Start med at opskrive data i to 1D arrays - et array for l og et array for T. Nu kan du udregne g vha. den givne funktion, og lave statistik på denne.

Dataen læses ind med

l = [51.2 59.7 68.2 79.7 88.3];

T = [1.448 1.566 1.669 1.804 1.896];

og tyngdeaccelerationskoefficienten udregnes ved

g = 4*pi^2*l./T.^2

Herefter kan man finde gennemsnittet, SD og SDOM af værdierne vha. de indbyggede MATLAB funktioner

g_mean = mean(g)

sigma = std(g)

sigma_mean = std(g)/sqrt(length(g))

Diskrebansen findes bare som forskellen,

g_right = 979.6;

discrepancy = g_mean - g_right

Prøv at lave en udregning hvor du ændrer din værdi for $l$ en lille smule, og se hvor langt du så skal ændre den for at få en $g$ der passer med den korrekte $g$. Det kan evt. være en god idé at se lidt på for-løkker i MATLAB, så programmet kan iterere for dig.

Plot der viser afstanden mellem den accepterede $g$ og den yderste SDOM grænse på den fundne $g$ når der er lagt $x$ til $l$.

For at finde hvor stor fejlen skulle være for at den rigtige værdi er inden for SDOM, kan man lave et loop der lægger lidt ekstra til $l$, og så udregner forskellen mellem den fundne $g$ og den rigtige,

x = 0.1:0.1:1.2;

for i=1:length(x)

g_new = 4*pi^2*(l+x(i))./T.^2;

discr_new(i) = mean(g_new) + std(g_new)/sqrt(length(g_new)) - g_right;

end

Den nemmeste måde at se resultatet på er ved at plotte det

plot(x,discr_new,'.b')

Man kan se at den skærer 0 ca. omkring 1 cm (se figuren), hvilket netop var hvad der skulle vises.

Man kunne selvfølgelig også sammenligne den teoretiske værdi for $g$ med gennemsnits-værdien for sin udregnede $g$, hvilket selvfølgelig bare skal gøres ved at erstatte discr_new udregningen ovenfor med

discr_new(i) = mean(g_new) - g_right;

For at rette den systematiske fejl kan man bare lægge den ene cm til sine $l$-værdier,

l = l + 1

og så udregner man sin $g$-værdi forfra, ligesom det blev gjort før:

g = 4*pi^2*l./T.^2;

g_mean = mean(g)

sigma_mean = std(g)/sqrt(length(g))

Plot af to forskellige Gauss-funktioner.

Først definerer man sin $x$-akse, f.eks. ved

x = -2:0.05:6;

Derefter giver man sine parametre

X1 = 2;

sigma1 = 1;

X2 = 3;

sigma2 = 0.3;

til sidst kan man så udregne sine funktionsværdier ved

y1 = 1/( sigma1*sqrt(2*pi) ) * exp(-(x-X1).^2/(2*sigma1^2));

y2 = 1/( sigma2*sqrt(2*pi) ) * exp(-(x-X2).^2/(2*sigma2^2));

Man kunne også definere en funktion for sin Gauss-fordeling ved definitionen

MinGauss = @(x,a,b) 1/( b*sqrt(2*pi) ) * exp(-(x-a).^2/(2*b^2));

og med denne kan man bare skrive

y1 = MinGauss(x,X1,sigma1);

y2 = MinGauss(x,X2,sigma2);

Til sidst plotter man de to funktioner sammen. Her kan man bruge at plot godt kan finde ud af at plotte flere funktioner af gangen:

plot(x,y1,'r-',x,y2,'b-')

Først og fremmest skal datasættet lige indlæses,

t = [63 58 74 78 70 74 75 82 68 69 76 62 72 88 65 81 79 77 66 76 86 72 79 ...

77 60 70 65 69 73 77 72 79 65 66 70 74 84 76 80 69];

Standardafvigelsen findes så simpelt med den indbyggede MATLAB funktion,

std(t)

For at binne med fire målinger af gangen er det nemmest at arbejde med et loop der gør det automatisk for dig. For at bruge et generelt antal punkter at binne med defineres tallet først:

ncol = 4;

Herefter startes der to tomme arrays op - et til gennemsnittene, og et til standardafvigelserne (zeros(n,m) laver et $n\times m$ array der indeholder nuller i alle elementer):

gn = zeros(1,10);

sd = zeros(1,10);

Et loop laver udregningerne:

for n = 0:9

gn(n+1) = sum(t(n*ncol+1:n*ncol + ncol))/ncol;

sd(n+1) = std(t(n*ncol+1:n*ncol + ncol));

end

Man kan så endelig se gn og sd ved simpelthen at printe dem ud via kommandolinjen.

Ud fra del (a) ville man forvente at SD er fordelt omkring den SD man fandt tidligere, og det er også det der sker.

De to histogrammer til sammenligning.

For at plotte histogrammer der kan sammenlignes skal man bruge histc, hvor man kan vælge bin-grænser manuelt. Histogrammerne laves ved først at definere grænserne for bins'ene ved

nbins = 5;  %antal bins

min = 58;  %position af første bin grænse

dist = 6;  %bredde af bins

max = min+nbins*dist;

edges = min:dist:max;

For at kunne plotte bins'ene med en korrekt førsteakse skal man også finde centrum af hvert bin. Dette gøres med

centers = min+dist*0.5:dist:max+dist*0.5;

Herefter udregnes histogrammerne med

hist40 = histc(t,edges);

hist10 = histc(gn,edges);

og de kan nu plottes med bar funktionen,

figure

bar(centers,hist40,0.5,'y')

hold on

bar(centers,hist10,0.3,'r')

Man kan tydeligt se at hist10 er mere centreret, som forventet.