9.5 Mer om programmering och funktioner
9.5.2 Viktad anpassning till r¨ at linje
N¨asta uppgift blir att med utg˚angspunkt i samma schema designa och skriva en funktion som g¨or en viktad anpassning till en r¨at linje, y=a + bx.
Algoritm
Parametrarna a och b, samt os¨akerheten i dessa ges av:
a =
Pwix2i Pwiyi−Pwixi Pwixiyi
∆ och b =
Pwi Pwixiyi−PwixiPwiyi
∆
σa= s
Pwix2i
∆ σb =
sP wi
∆ d¨ar wi = 1
σi2 och ∆ =Xwi
Xwix2i −Xwixi
2
9.5. MER OM PROGRAMMERING OCH FUNKTIONER 99 Gr¨anssnitt och syntax
Indata till funktionen ¨ar l¨ampligen tre vektorer, en med x-v¨arden, en med y-v¨arden och en med felen i y-v¨ardena. Vilka resultat vi vill f˚a ut ur funktionen kan variera en del. Vi kommer alltid att vilja se de anpassade v¨ardena f¨or parametrarna a och b, samt os¨akerheten i dessa.
Eventuellt vill vi ocks˚a f˚a en uppfattning om kvaliteten i anpassningen, det vill s¨aga chi-kvadratsumman per frihetsgrad och antalet frihetsgrader. Eftersom det, i motsats till lincorr ovan, inte ¨ar n¨odv¨andigt att r¨akna ut dessa senare tal f¨or att funktionen skall kunna k¨ora, s˚a l¨onar det sig att testa med nargout om dessa v¨arden verkligen ¨ar efterfr˚agade eller inte innan de ber¨aknas. Syntaxen f¨or funktionen blir d˚a: function [parA, sigmaA, parB, sigmaB, ChiSquared, Dof] =linfit(x, y, sigma).
Felkontroll
Ett m¨ojligt fel som m˚aste kollas ¨ar om de tre vektorerna som ¨ar argument till funktionen inte har samma l¨angd. Vidare b¨or man kontrollera att ∆ inte ¨ar noll s˚a att divisionerna ger ett odefinierat resultat. Om chi-kvadrat skall ber¨aknas m˚aste vi ocks˚a kontrollera att det finns minst tre m¨atta punkter, annars kommer divisionen n¨ar chi-kvadrat per frihetsgrad ber¨aknas att bli odefinierad, samt att inget av m¨atfelen felaktigt ¨ar 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−3, b = (389.99 ± 0.665), N = 6, χ2/d.o.f. = 2.38 9.5.3 Viktad anpassning och plottning av r¨at linje.
H¨ar skall vi skriva en funktion som levererar en graf med data plottade med felstaplar och i samma graf den r¨ata linjen som ¨ar ett resultat av anpassningen till en r¨at linje, med parametrarna f¨or anpassningen utskrivna.
Algoritm
Anpassningen till en r¨at linje ¨ar inte n˚agot problem, vi l˚ater bara denna funktion i sin tur anropa linfit som vi skrev i f¨orra avsnittet. Det g¨aller sedan ”bara” att f˚a grafen att se snygg ut utan att g¨ora n˚agra antaganden i f¨orv¨ag om hur data ser ut. Det kan man
˚astadkomma genom att f¨orst plotta data, sedan flytta ut axlarna med 10% av det v¨arde som den automatiska skalningen ger f¨or att se till att alla punkter ligger v¨al 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¨or att skriva ut parametrarna formaterar vi med formatet %0.3g f¨or 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¨alle och v¨aljer d¨arf¨or att b¨orja skriva 5% av grafens totala bredd in fr˚an v¨ansterkanten och ¨oversta raden 5% fr˚an toppen. Om riktningskoefficienten fr˚an anpassningen ¨ar negativ kommer datapunkter att ligga i det ¨ovre v¨anstra h¨ornet, vi v¨aljer d˚a ist¨allet att i horisontell ledd b¨orja att skriva i mitten av grafen. (Ett alternativ till denna algoritm kan vara att l¨agga till extra 20% i h¨ojdled n¨ar vi skalar om vertikalt och anv¨anda
denna area till att placera in texten.) Vi b¨or d¨aremot undvika att anv¨anda kommandon som p˚averkar var grafen plottas, som figure(#) och subplot, det ¨ar b¨attre att det kallande programmet g¨or s˚adana val, eftersom ¨onskem˚alen kan vara vitt skilda f¨or olika till¨ampningar.
Gr¨anssnitt och syntax
Gr¨anssnittet f¨or den h¨ar funktionen ¨ar detsamma som f¨or 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 ¨aven ger en str¨ang som anger hur grafen skall plottas av plot-kommandot, som till exempel ’-or’, men i f¨orsta omg˚angen bryr vi oss inte om det. D¨aremot m˚aste vi skicka med tre textstr¨angar som ger text-str¨angar f¨or titel, x-etikett och y-etikett. Om vi vill kunna v¨alja om parameterv¨arden skrivs ut i plotten eller inte s˚a m˚aste vi s¨atta in flaggor f¨or detta i argumenten till funktionsanropen, parflag som ¨ar 1 om vi vill att parameterv¨ardena skall skrivas in och 0 annars, chiflag styr p˚a samma s¨att 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¨anda title och xlabel som namn p˚a variablerna). Utdata fr˚an linplot ¨ar de samma som f¨or linfit, allts˚a i f¨orsta hand v¨ardet p˚a parametrarna och os¨akerheten i dessa, samt som tillval chi-kvadratsumman och antalet frihetsgrader. Vi kan dessutom t¨anka oss att linplot ibland kommer att kallas utan att den kallande rutinen anv¨ander utdata alls, i de fall det intressanta bara ¨ar att f˚a en graf. Syntaxen f¨or 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 ¨over gr¨anssnittet mellan dessa b¨agge funktioner. I princip ¨ar det ju inget som hindrar linplot fr˚an att kalla linfit med alla sex utdata parametrarna i anropet. Men om linplot inte ¨ar kallad med parametrarna ChiSquared och Dof, och om chiflag ¨ar noll s˚a beh¨over man ju inte ber¨akna chi-kvadratsumman i linfit. Att ¨and˚a kalla linfit med alla sex utdata parametrar betyder d˚a att programmet spenderar tid p˚a att ber¨akna saker som aldrig kommer att anv¨andas, vilket ¨ar sl¨oseri. Vi m˚aste allts˚a i linplot testa om s˚a ¨ar fallet eller inte och bara anropa linfit med alla sex utdata parametrar om det ¨ar n¨odv¨andigt, annars anv¨ander vi [parA, sigmaA, parB, SigmaB]
Felkontroll
Det fel vi beh¨over testa f¨or ¨ar om vektorerna x, y och sigma har samma l¨angd, i annat fall kraschar errorbar-kommandot. Man kan tycka att vi inte beh¨over testa f¨or det i linplot eftersom linfit testar f¨or det. Men man skall f¨ors¨oka att skriva varje funktion s˚a komplett och sluten som m¨ojligt. Det kan ju h¨anda att vi i framtiden vill anv¨anda linplot s˚a att den kallas av n˚agon annan funktion ¨an linfit, i en s˚adan situation kan det vara om¨ojligt att minnas att vi m˚aste kolla om den nya funktionen g¨or den n¨odv¨andiga testen.
9.5. MER OM PROGRAMMERING OCH FUNKTIONER 101 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
disp(’YYYY’)
[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¨att st¨alle n¨ar riktningskoefficienten ¨ar negativ, och att vi kan sl˚a av och p˚a utskrift av parameterv¨arden med flaggorna som det ¨ar t¨ankt.
9.6 Ovningsuppgifter ¨
1. Rita f¨oljande funktioner, prova med olika kombinationer av mesh, surf, surfl och shadint interp, samt testa olika f¨argkartor.
z = q
x2+ y2 − 10 ≤ x ≤ 10, −10 ≤ y ≤ 10 z =
q
x2+ y2
− 2π ≤ x ≤ 2π, −2π ≤ y ≤ 2π
2. Kronan p˚a verket blir att skriva en funktion som med matrismetoden g¨or en minsta kvadratanpassning till ett polynom med godtyckligt gradtal. Indata till funktionen ¨ar en vektor med y-v¨arden, en med x-v¨arden och en vektor med os¨akerheten i y-v¨ardena.
Dessutom m˚aste vi vid anropet ange till vilket gradtal anpassningen skall g¨oras. Utdata
¨
ar vektorn med parameterv¨ardena och covariansmatrisen. Den felkontroll som beh¨over g¨oras ¨ar att antalet y-v¨arden ¨ar minst lika stort som antalet parametrar som skall anpassas och att inget av felen ¨ar noll.
Kapitel 10
Svar och l¨ osningar
10.2 Svar till kapitel 2
1. a) y=3*x.^2+4*x (y = 7 39 95 175 )
b) yPrim = yPrim = 6*x + 4 (yPrim = 10 22 34 46) 2. y=2*x.^3 + 4*x - 3 (y = 3.0000 49.5938 131.2380 210.0720) 3. a1 = [2;-5;2] ; a2 = [3 ; 2 ; -4]; a3 = [1 ; 0 ; 2]
A = [a1 a2 a3]
A2 = [a2’ ; a3’ ; a1’]
4. A + B =
4 -7 9 2
-1 -5 -1 7
1 -2 6 12
A-B =
-2 1 5 -2
5 -3 -1 5
-5 2 0 -2
A’*B =
-9 -2 -4 -10
3 16 -6 -10
33 -33 23 34
-3 -16 15 41
5. x = alfa*pi/180 (x = 0 0.5236 1.0472 1.5708)
s = sin(x) ( s = 0 0.5000 0.8660 1.0000 )
c=(1-s.^2).^0.5 ( c = 1.0000 0.8660 0.5000 0 )
A=[alfa’ s’ c’]
( A =
0 0 1.0000
30.0000 0.5000 0.8660 60.0000 0.8660 0.5000
90.0000 1.0000 0 )
103
6. massa = [2.0 1.5 3.0 2.0 2.5];
hastighet = [ 2.0 4.0 -3.0 1.0 6.0];
T = 0.5*massa.*hastighet.^2
(T = 4.0000 12.0000 13.5000 1.0000 45.0000 totT = sum(T) (totT = 75.5000)
7. r = 6*rand(1,10)
10.3 Svar till kapitel 3
1. Ett f¨orslag till l¨osning kan vara:
function y = slumptal()
%
% producerar ett slumpat heltal mellan 1 och 100
%
y = 1 + floor(100*rand(1,1));
2. Ett f¨orslag till l¨osning ¨ar:
function tipsrad()
%
% genererar en slumpm\"assig tipsrad
%
r = 3*rand(1,13) for i = 1:13
if (r(i) < 1) disp(’1’)
elseif (r(i) < 2) disp(’x’)
else disp(’2’)
end end
3. Ett exempel ¨ar:
function y=Newton(x)
% function to compute square root of x, using Newtons formula
% y(k+1) = 0.5*(y(k) + x/y(k)) with y(1) = 1
%
% Sten Hellman 050301 yold = 1;
delta = 1;
index = 0;
while(delta > 0.00001) index = index + 1;
10.5. SVAR TILL KAPITEL 5 105 ynew = 0.5* ( yold + x/yold);
delta = abs(yold-ynew) yold=ynew
end
y = ( round(10000*ynew) / 10000)
10.5 Svar till kapitel 5
1. Felet med programmet ¨ar att variabeln x inte ¨andras n¨ar vi har hittat en primtalsfaktor.
Varje g˚ang vi till exempel testar om 9 ¨ar j¨amnt delbart med 3 s˚a stegas AntalFaktorer upp med en enhet, men x ¨ar fortfarande 9 n¨asta g˚ang vi testar. N¨ar vi hittar en primfaktor till x m˚aste vi dividera bort den ur x innan vi forts¨atter att s¨oka fler faktorer, till exempel s˚a h¨ar:
AntalFaktorer=0;
2. Det vi ser ¨ar hur serieutvecklingen av sinusfunktionen successivt n¨armar sig funktionen vartefter man tar med fler och fler termer i summan. Grafen b¨or se ut ungef¨ar s˚a h¨ar:
0 1 2 3 4 5 6 7
Vilket f˚as genom att koda:
>> I= [3 4 ; -2 1];
10.8. SVAR TILL KAPITEL 8 107 vilket f˚as ur: >> clear all
>> B= [22 ; -21 ; 17];
>> A = [4 -5 3 ; 1 5 -6 ; 3 1 4];
>> inv(A)*B ans =
2.0000 -1.0000 3.0000 3. Det r¨acker med att l¨agga till raderna
Deltaa = sqrt(sumwx2/delta);
Deltab = sqrt(sumw/delta);
10.8 Svar till kapitel 8
1. Ett exempel ¨ar:
% Primtal-final.m
% to compute factors to arbitrary number x
% Sten Hellman 2005-03-03
%
%
clear all;
namn=input(’Hej och valkommen till PRIMFAKT - \n Vad heter du? ’,’s’);
out=[’Hej ’ namn ’\n Mata nu in det tal du vill soka faktorer for \n’];
x=input(out);
xinitial =x;
% startvrdet p x mste sparas eftersom x kan frndras senare i
% programmet AntalFaktorer=0;
i=2;
while(i<= x) kvot = x/i;
rest = i*floor(kvot)-x;
if(rest ==0)
AntalFaktorer = AntalFaktorer + 1;
Faktor(AntalFaktorer) = i;
x = kvot;
else
if(i==2) i=i+1;
else i=i+2;
end end end
out1 = sprintf(’%g har %g primfaktorer:\n’,xinitial,AntalFaktorer);
out2 = sprintf(’%g, ’,Faktor);
out3 = sprintf(’\b\b’);
out = [out1 out2 out3];
disp(out)
2. En m¨ojlig l¨osning ¨ar:
% exponent.m
% computes approximation to exp(x) from serial
% expansion
% Sten Hellman 2005-03-03 clear all;
disp(’En approximation till e berknas genom serieutveckling’);
Termer=input(’Upp till vilken term skall summeringen g? \n’);
Faktor=1;
Approx = 1;
% this establishes the zero-term for index=1:Termer
Faktor = Faktor * index;
Approx = Approx + 1/Faktor;
end
Diff = exp(1) - Approx;
out = sprintf(’Efter %g termer r approximationen= %g, skillnaden r %g’,index,Approx,Diff);
disp(out)
3. Eftersom vi inte p˚a f¨orhand vet hur m˚anga termer som kommer att beh¨ovas s˚a m˚aste vi ers¨atta for-loopen med en while-loop. Vi justerar dessutom in/output satserna:
% exponent2.m
% computes approximation to exp(x) from serial
% expansion
% Sten Hellman 2005-03-03 clear all;
disp(’En approximation till e berknas genom serieutveckling’);
Delta=input(’Upp till vilken noggrannhet summeringen g? \n’);
10.8. SVAR TILL KAPITEL 8 109 Faktor=1;
Approx = 1;
Diff = exp(1) - Approx;
index = 0;
% this establishes the zero-term while (Diff >= Delta)
index = index + 1;
Faktor = Faktor * index;
Approx = Approx + 1/Faktor;
Diff = exp(1) - Approx;
end
out = sprintf(’Efter %g termer r approximationen= %g, skillnaden r %g’,index,Approx,Diff);
disp(out)
\, 89
st¨ang av eko av kommandon, 12 ekvationer
i Word, 52 110
INDEX 111
fler grafer i samma f¨onster, 70 fler kurvor i samma graf, 69 flytta till andra program, 68
klippa ut text nollst¨allen f¨or polynom, 79 num2str, 88
INDEX 113
st¨ang av utskrift av resultat av kom-mandon, 19
fet stil, 45 figurer, 51 figurtext, 52 formatmallar, 47 fotnoter, 49
h¨ogerjusterad text, 45 infoga symboler, 52 justera text, 45 klipp ut text, 44 klistra in text, 44 kopiera text, 44 kursiv stil, 45 l¨agga in bilder, 69 mallhanterare, 47
marginaljusterad text, 45 markera text, 43
numrerade listor, 49 punktlistor, 49 rubriker, 46 sidbrytning, 47 sidnumrering, 48 stavningskontroll, 48 tabeller, 50
tabelltext, 51 typsnitt, 45
understruken text, 45 v¨ansterjusterad text, 45 workspace, 13
l¨as in, 35 rensa, 76 spara, 75
spara workspace till diskspace, 35 workspace browser, 75
xlabel, 65 ylabel, 65 zeros, 24