Matlab: Kurvanpassning

[ElectricMan]

Vi har mätt vattenflöde i ett rör där vi har varierat diametern.
Nu vill vi ha en anpassad funktion till vårat data för att bestämma, enkelt förklarat, k-värdet (lutningen) på funktionen.

Efter test i Excel med hjälp av linjärisering av datan fick vi fram att k-värdet är 4,25.

Men nu är det MATLAB som gäller.

Jag får det att fungera skitbra för ett annat samband, men här vill det inte alls anpassa sig efter funktionen...
Känns som om jag har testat det mesta nu, men jag får verkligen inte ut värdet i MATLAB.

Vi har fått en liten guide för hur man gör detta, men den är _minst_ sagt bristfällig...

Såhär ser koden ut;

Kod: Markera allt

t = 30;                               % Tid (s)
m = [6    25   39   15   2];          % Vattenmassa per 30 sekunder. (g)
Q2 = ((m./1000)/998.2)/t              % Vattenflöde (m^3/s)
d = [1.07 1.54 1.75 1.37 0.88]./1000  % Rördiametrar (m)

plot(d,Q2,'o')
xlabel('d [m]');
ylabel('Q [m^3/s]');
title('Uppmätt flöde för olika diametrar.');


E = @(b) sum((d-b(1)*Q2.^b(2)).^2);
c = fminsearch(@(b) E(b), [5 1])

l_fin=linspace(0,0.001,200);
t_fin=c(1)*l_fin.^c(2);
plot(d,Q2,'o',l_fin,t_fin,'r');
xlabel('d [m]')
ylabel('Q [m^3/s]');
title('Funktion anpassad för \pi1 (Flöde för olika diametrar)');
legend('Mätdata','Anpassad funktion');

Ser ni var felet finns? Finns det bättre sätta att göra detta på?

[Fisksmisk Ola]

Det låter som ett typiskt jobb för polyfit:
http://www.mathworks.se/help/matlab/ref/polyfit.html

Eller så kanske jag missuppfattade vad som eftersöktes?


/Olof

Edit:
Ja, det gjorde jag.

Edit:
Fast jag tror ändå polyfit skulle funka.
skriv om:
y = a*x^b
som
log(y) = log(a)+b*log(x)

Här kommer lite kod:

Kod: Markera allt

format long

t = 30;                               % Tid (s)
m = [6    25   39   15   2];          % Vattenmassa per 30 sekunder. (g)
Q2 = ((m./1000)/998.2)/t              % Vattenflöde (m^3/s)
d = [1.07 1.54 1.75 1.37 0.88]./1000  % Rördiametrar (m)
plot(d,Q2,'o'), hold on

a = polyfit(log(d), log(Q2), 1);
c = [exp(a(2)), a(1)]

x = 0:0.00001:0.002;
fitfunc = c(1)*x.^c(2);
plot(x, fitfunc)

I 'c' borde de värden du söker finnas, det framgår ganska tydligt vilket som är vilket.
Tyvärr hittade jag inte riktigt vad som var galet i din kod... Men det c som din kod genererar är ganska olikt det jag får.

[ElectricMan]

Fast den där ger ju k-värdet som 7,63*10^5. Vilket är aningen mer än 4,25.

Nåja, tror det blir Excel.
Orkar inte rackla med MATLAB mer. Det var mest på skoj jag ville testa använda MATLAB.

[blueint]

Perl?

[SvenW]

Märkligt. För min del var det ett stort steg framåt när jag gick från Excel till Matlab-clonen Octave.
Som jag ser det är Excel mest till för ekonomikalkyler och liknande. Inte för tekniska beräkingar
och absolut inte för analys av data.

[Fisksmisk Ola]

Fast den där ger ju k-värdet som 7,63*10^5. Vilket är aningen mer än 4,25.

Nåja, tror det blir Excel.
Orkar inte rackla med MATLAB mer. Det var mest på skoj jag ville testa använda MATLAB.

Njae, då tittar du nog på fel värde. Den andra komponenten i c-vektorn (c(2)), dvs den som är variaveln 'b' i

y = a*x^b

är = 10^5*0,00004257... = 4,257...

Jag får utskriften(fast jag har förkortat den lite):

Kod: Markera allt

Q2 =

...

d =

...

c =

  1.0e+005 *

   7.639613103679824   0.000042570989994

där det sista som skrivs ut, alltså c, är en vektor med två element. Det du kallar k-värdet borde alltså vara c(2).
I figuren ser man också att funktionen 'fitfunc' ligger väldigt nära mätvärdena.

Det jag gjort är i princip transformerat funktionen till en linjär funktion och anpassat parametrarna till den, och sedan transformerat tillbaka.
(Transformerat som i 'ln' och sedan 'e^'.)


/Olof

[ElectricMan]

Doh...

Att jag inte såg det där!


SvenW;
Just därför jag ville använda MATLAB. Men när det inte fungerar så blir det svårt.