vilse i octave, data till ekvation

[jsiei97]

Hej

Satt och pillade lite med NTC motstånd, en underhållande detalj är att dessa är olinjära i sin natur.
Hur som helst, tillverkarna av ntc motstånd ger ofta en datamängd som i praktiken berättar att vid den här temperaturen har vi detta motstånd.

Men eftersom tabeller är lite småtråkiga så har jag tidigare stoppat in utvalda data in i min miniräknare (en TI-89),
och sedan låtit den göra en approximation på vilken ekvation som passa bäst mot denna datamängd.
Och sedan kontrollerat detta genom att både rita upp datapunkterna och ekvationen, så man enkelt kan se om det stämmer eller inte.
(Det beror mycket på vilken data, om den är linjär eller olinjär osv osv)

Denna gången fick jag dock en liten undran om man inte kan använda octave (som verkar bra häftigt) att göra detta åt mig istället?

Så jag har börjat med att installera QtOctave (och en massa paket som heter octave-XXX),
och sedan klipp och klistrat in data in i en tabell vilket kan återskapas med.

Kod: Markera allt

ntcC=[0; 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; ]
ntcR=[28380; 27130; 25940; 24810; 23740; 22720; 21750; 20830; 19950; 19120; 18320; 17570; 16840; 16160; 15500; 14880; 14280; 13710; 13170; 12650; 12160; 11690; 11240; 10810; 10390; 10000; 9623; 9263; 8918; 8588; 8272; 7970; 7680; 7402; 7136; 6881; 6636; 6402; 6177; 5961; 5754; 5555; 5365; 5182; 5006; 4837; 4674; 4518; 4368; 4224; 4085; ]

I ntcC lagrar jag grader Celcius och i ntcR den resistans som ntc motståndet har vi denna temperatur (dessa värden kommer direkt från databladet)

och sedan kan jag rita upp datamängden med t.ex.

Kod: Markera allt

plot(ntcC,ntcR)
bar(ntcC,ntcR)
print -dpng '~/ntc_raw.png'

Vilket ger mig detta fina diagram

ntc_raw.png (3.01 KiB) Visad 2128 gånger

Men sedan börjar det roliga, hur gör jag för att octave ska titta på datan och göra en "exponential regression", "linear regression" m.fl. ?
(och skapar ett nytt diagram baserat på den skapade ekvationen?)

/Johan

[jsiei97]

En notis, om jag gör som jag brukar med TI-89:an och använder ExpReg så får jag

Kod: Markera allt

y=27247*0.962^x

dvs

Kod: Markera allt

ntcR=27247*0.962^ntcC

Vilket är hyfsat nära (trots att jag avrundade lite för att det skulle gå att läsa)...

Frågan är dock fortfarande hur jag får fram ovan med octave?

[Konservburk]

polyfit och polyval är användbara...

Kod: Markera allt

% linjär
plot(ntcC,[polyval(polyfit(ntcC, ntcR, 1), ntcC) ntcR])

% exponentiell
plot(ntcC,[exp(polyval(polyfit(ntcC, log(ntcR), 1), ntcC)) ntcR])

Öka på ordningen om du vill ha en bättre kurvanpassning:

Kod: Markera allt

plot(ntcC,[exp(polyval(polyfit(ntcC, log(ntcR), 3), ntcC)) ntcR])

Nu frågar du dig förmodligen hur du får ut resultatet som ekvationer:

Kod: Markera allt

> polyfit(ntcC, ntcR, 1)
ans =

     -446.217194570136      23198.9004524887

ntcR = 23199 - 446 * ntcC

Kod: Markera allt

> exp(polyfit(ntcC, log(ntcR), 1))
ans =

     0.962066577426195      26974.7761431943

ntcR = 26975 * 0.962 ^ ntcC

[jsiei97]

polyfit och polyval är användbara...

Det där var trevliga saker

Argument listan känns inte helt intuitiv, men det är nog en vanesak... (som med så mycket annat)

Det jag tänker på är exp och log, att du lägger en exp före för att berätta att vi vill göra en exponentiell kan verka logiskt,
men log() sedan i mitten verkar inte helt uppenbar.

[Konservburk]

Det där med log() och exp() är bara en mappning fram och tillbaka före och efter en vanlig linjär kurvanpassning.

polyfit(X, Y, 1)
ger oss koeffiecenterna A och B i ekvationen
Y = A + B * X

Men du vill ha en ekvation på formen:
y = a * b ^ x

Den går att skriva om som:
log(y) = log(a) + log(b) * x

polyfit(x, log(y), 1)
ger oss då koefficenterna log(a) och log(b)

Vi kör därför exp() på resultatet för att istället få ut a och b.

[jsiei97]

Nu börjar det klarna,
alltid lättare att förstå när man ser ett par mellanled

[Konservburk]

Jag kom just på att man faktiskt kan göra detta utan polyfit och polyval (eller andra specialfunktioner) och dessutom flera gånger snabbare. Det spelar iofs inte någon roll med så här små datamängder, men det kan vara intressant ändå.

Kod: Markera allt

% Alternativ A
FIT = polyfit(ntcC, ntcR, 1);
VAL = polyval(FIT, ntcC);

% Alternativ B
x = [ntcC, ntcC.^0];
fit = x \ ntcR;
val = x * fit;

Tidstest:

Kod: Markera allt

> tic; for i = 1:10000; FIT = polyfit(ntcC, ntcR, 1); VAL = polyval(FIT, ntcC); end; toc
Elapsed time is 5.59800102957524 seconds.

> tic; for i = 1:10000; x = [ntcC, ntcC.^0]; fit = x \ ntcR; val = x * fit; end; toc
Elapsed time is 0.679437005426735 seconds.

[jsiei97]

Bra mångsidigt program

Satt och roade mig med att titta på delmängder (jag svarar mig själv så jag kan titta här nästa gång ifall jag glömmer hur det gick till ),
säg att man vill använda så enkel ekvation som möjligt och vet inom vilket område man ska jobba i.
Då är det ofta bättre att "zoma" in så mycket som möjligt och sedan där göra sin uppskattning, och i den mindre datamängden är kanske t.o.m. en linjär uppskattning bra nog (men då får man zoma in rätt mycket)

Låt oss använda all data som exempel.
Här har vi temperatur i första kolumnen, och sedan har vi 3st resistanser som verkar vara övre och undre gräns samt normal fallet.

Kod: Markera allt

ntc = [-40, 220800, 230400, 240300; -39, 208000, 216900, 226100; -38, 196000, 204200, 212800; -37, 184800, 192500, 200400; -36, 174300, 181400, 188800; -35, 164500, 171100, 177900; -34, 155300, 161400, 167800; -33, 146700, 152400, 158300; -32, 138600, 143900, 149400; -31, 131100, 136000, 141100; -30, 124000, 128500, 133300; -29, 117300, 121600, 126000; -28, 111000, 115000, 119100; -27, 105200, 108800, 112700; -26, 99630, 103100, 106600; -25, 94420, 97630, 100900; -24, 89530, 92510, 95590; -23, 84920, 87700, 90570; -22, 80570, 83170, 85840; -21, 76480, 78910, 81390; -20, 72630, 74890, 77210; -19, 68990, 71100, 73260; -18, 65560, 67530, 69540; -17, 62320, 64160, 66040; -16, 59270, 60980, 62740; -15, 56380, 57980, 59620; -14, 53660, 55150, 56680; -13, 51080, 52480, 53910; -12, 48650, 49950, 51290; -11, 46350, 47570, 48810; -10, 44170, 45310, 46470; -9, 42110, 43180, 44260; -8, 40160, 41160, 42170; -7, 38320, 39240, 40190; -6, 36560, 37430, 38320; -5, 34910, 35720, 36540; -4, 33330, 34090, 34860; -3, 31840, 32550, 33270; -2, 30420, 31090, 31760; -1, 29080, 29700, 30330; 0, 27800, 28380, 28970; 1, 26590, 27130, 27680; 2, 25430, 25940, 26450; 3, 24340, 24810, 25290; 4, 23300, 23740, 24190; 5, 22310, 22720, 23140; 6, 21360, 21750, 22140; 7, 20470, 20830, 21190; 8, 19610, 19950, 20290; 9, 18800, 19120, 19440; 10, 18030, 18320, 18620; 11, 17290, 17570, 17840; 12, 16590, 16840, 17100; 13, 15920, 16160, 16400; 14, 15280, 15500, 15730; 15, 14670, 14880, 15090; 16, 14090, 14280, 14480; 17, 13530, 13710, 13900; 18, 13000, 13170, 13340; 19, 12500, 12650, 12810; 20, 12010, 12160, 12310; 21, 11550, 11690, 11820; 22, 11110, 11240, 11360; 23, 10690, 10810, 10920; 24, 10290, 10390, 10500; 25, 9900, 10000, 10100; 26, 9523, 9623, 9723; 27, 9163, 9263, 9363; 28, 8819, 8918, 9018; 29, 8489, 8588, 8687; 30, 8174, 8272, 8371; 31, 7872, 7970, 8068; 32, 7583, 7680, 7777; 33, 7306, 7402, 7499; 34, 7041, 7136, 7232; 35, 6787, 6881, 6976; 36, 6543, 6636, 6730; 37, 6310, 6402, 6495; 38, 6086, 6177, 6269; 39, 5871, 5961, 6052; 40, 5665, 5754, 5844; 41, 5468, 5555, 5644; 42, 5278, 5365, 5452; 43, 5096, 5182, 5268; 44, 4921, 5006, 5091; 45, 4754, 4837, 4921; 46, 4592, 4674, 4757; 47, 4438, 4518, 4600; 48, 4289, 4368, 4449; 49, 4146, 4224, 4303; 50, 4008, 4085, 4163; 51, 3876, 3952, 4029; 52, 3749, 3823, 3899; 53, 3627, 3700, 3774; 54, 3509, 3581, 3654; 55, 3396, 3466, 3538; 56, 3287, 3356, 3427; 57, 3182, 3250, 3320; 58, 3081, 3148, 3216; 59, 2984, 3050, 3117; 60, 2890, 2955, 3021; 61, 2800, 2863, 2928; 62, 2713, 2775, 2839; 63, 2629, 2691, 2753; 64, 2548, 2609, 2670; 65, 2470, 2530, 2590; 66, 2395, 2454, 2513; 67, 2323, 2380, 2439; 68, 2253, 2309, 2367; 69, 2185, 2241, 2297; 70, 2120, 2174, 2230; 71, 2057, 2111, 2165; 72, 1997, 2049, 2102; 73, 1938, 1989, 2041; 74, 1881, 1931, 1983; 75, 1826, 1876, 1926; 76, 1773, 1822, 1871; 77, 1722, 1770, 1819; 78, 1673, 1720, 1767; 79, 1625, 1671, 1718; 80, 1579, 1624, 1670; 81, 1534, 1578, 1623; 82, 1491, 1534, 1578; 83, 1449, 1491, 1535; 84, 1408, 1450, 1493; 85, 1369, 1410, 1452; 86, 1331, 1372, 1413; 87, 1295, 1334, 1375; 88, 1259, 1298, 1338; 89, 1225, 1263, 1302; 90, 1192, 1229, 1267; ] 

Dvs här finns data för -40 till +90 grader, men om jag bara vill titta på -10 till +30 grader så verkar jag kunna göra så här.

Kod: Markera allt

start =  31
stopp =  71
plot(ntc([start:1:stopp],1), [ntc([start:1:stopp],2),ntc([start:1:stopp],3),ntc([start:1:stopp],4)])

Eller kanske kanske ännu snävare mellan +15 och +30 (dvs den temperatur som normalt finns inomhus),
då behöver jag bara uppdatera "start" och sedan räkna igen..

Kod: Markera allt

start =  56
plot(ntc([start:1:stopp],1), [ntc([start:1:stopp],2),ntc([start:1:stopp],3),ntc([start:1:stopp],4)])

Jisses, jag tror att jag börjar gilla det här programmet
(Både enkelt och hopplöst svårt på samma gång )

[Konservburk]

Det går såklart att använda indexeringen när du väljer start och stopp precis som du har gjort. Men om du nu är intresserad av intervallet -10 till 30 grader så blir det lite smidigare att använda sig av de faktiska värdena istället. Ett annat tips är att använda 2:end så slipper du lista alla kolumnerna separat när du plottar:

Kod: Markera allt

intervall = ntc(:, 1) >= -10 & ntc(:, 1) <= 30;
plot(ntc(intervall, 1), ntc(intervall, 2:end))

[jsiei97]

Jag vet inte varför jag får en "perl" känsla,
det finns alltid mer än ett sätt att göra det på

(fast syntaxen är snyggare och mer fantasifull)