Curve fitting a power law function

9 views (last 30 days)
Yen Tien Yap
Yen Tien Yap on 24 Sep 2021
Answered: Alan Stevens on 24 Sep 2021
%% 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.

Answers (1)

Alan Stevens
Alan Stevens on 24 Sep 2021
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

Categories

Find more on MATLAB in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!