Curve fitting a power law function

%% Polyacrylamide solution curve fitting analysis 
rho=1000; %[kg/m^3]
alpha=15; %[deg]
R=1.2*10^-3; %[m]
L=1.1; %[m]
D=10*10^-3; %[m]
g=9.81; %[ms-2]

h0=0.654; %[m]
h_t=[0.654;0.628;0.604;0.582;0.56;0.54;0.52;0.501;0.482;0.465;0.447;0.43;0.415;...
    0.399;0.384;0.369;0.356;0.343;0.329;0.316;0.304;0.293;0.282;0.271;0.261;0.251;...
    0.242;0.233;0.225;0.216;0.209]; %[m]
t=transpose([0:60:1800]); %[s]

syms n k 
C1=-((pi*R^3)/(((1/n)+3)*tand(alpha)*D))*(rho*g*R/(2*L*k))^(1/n)
C2=0.654^(2*n-1/n)
h2_t=((2*n-1/n)*(C1.*t+C2)).^(n/2*n-1)

figure(4)
plot(t,h_t,'r','LineWidth',1.5)
hold on
scatter(t,h2_t)

Hi, I want to have a curve fitting (regression) plot and find the values of n and k for the power law function. How can I do this? Thank you.

Answer 1

[Alan Stevens]

Like this?

h0=0.654; %[m] This seems to be unused
h_t=[0.654;0.628;0.604;0.582;0.56;0.54;0.52;0.501;0.482;0.465;0.447;0.43;0.415;...
    0.399;0.384;0.369;0.356;0.343;0.329;0.316;0.304;0.293;0.282;0.271;0.261;0.251;...
    0.242;0.233;0.225;0.216;0.209]; %[m]
t=transpose(0:60:1800); %[s]

X0 = [1, 1];  % Initial guesses [n0, k0]
X = fminsearch(@(X) fitfn(X, t, h_t), X0);
n = X(1); k = X(2);
disp(['n = ' num2str(n)])

n = 2.2345

disp(['k = ' num2str(k)])

k = 0.00020214

figure(4)
plot(t,h_t,':r','LineWidth',1.5)
hold on
scatter(t,h2_t(X,t))

function F = fitfn(X, t, h_t)
         F = norm(h2_t(X,t) - h_t);
end

function h2t = h2_t(X,t)
        %% Polyacrylamide solution curve fitting analysis 
        rho=1000; %[kg/m^3]
        alpha=15; %[deg]
        R=1.2*10^-3; %[m]
        L=1.1; %[m]
        D=10*10^-3; %[m]
        g=9.81; %[ms-2]
        n = X(1); k = X(2);
        C1=-((pi*R^3)/(((1/n)+3)*tand(alpha)*D))*(rho*g*R/(2*L*k))^(1/n);
        C2=0.654^(2*n-1/n);
        h2t=((2*n-1/n)*(C1.*t+C2)).^(n/2*n-1);
end