Curve fitting a power law function

9 visualizaciones (últimos 30 días)
Yen Tien Yap
Yen Tien Yap el 24 de Sept. de 2021
Respondida: Alan Stevens el 24 de Sept. de 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.

Respuestas (1)

Alan Stevens
Alan Stevens el 24 de Sept. de 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

Categorías

Más información sobre MATLAB en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by