Curve fitting for jonscher power law
4 views (last 30 days)
Show older comments
Navaneeth Tejasvi Mysore Nagendra
on 10 May 2022
Commented: Navaneeth Tejasvi Mysore Nagendra
on 10 May 2022
I am trying to fit jonscher power law equation, sigma_ac = sigma_dc + A w^n in which I have values for sigma_ac and w. I had tried to fit the plot using fittype function and not able to understand what 'StartPoint ' in the fit is and what values should I give. Based on this my fit is varying a lot and I am not able to find a suitable value for this.
Any insight and help is appreciable.
my program is;
clear;close all;clc;
opts = spreadsheetImportOptions("NumVariables", 6);
% Specify sheet and range
opts.Sheet = "Sheet1";
opts.DataRange = "A2:F202";
% Specify column names and types
opts.VariableNames = ["f1", "sigmaAcZT1", "f2", "SIGMAACZT2", "F2", "SIGMAACZT3"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double"];
% Import the data
data = readtable("C:\Users\navaneeth\Downloads\IONIC CONDUCTIVITY.xlsx", opts, "UseExcel", false);
% Convert to output type
data = table2array(data);
% Clear temporary variables
clear opts
w = 2*pi*data(:,1);%w = 2*pi*f
AC = data(:,2);% AC conductivity
ft = fittype('c+a*w^b','dependent',{'Conductivity'},'independent',{'w'},'coefficients',{'a','b','c'});% AC_conductivity = Dc_conductivity +A*w^n
f = fit(w,AC,ft,'StartPoint',[0,1,1]);
plot(f,w,AC)
4 Comments
Torsten
on 10 May 2022
Edited: Torsten
on 10 May 2022
do you know how can I fix the scaling on the plot to its original powers,
The coefficients c, a and b for the original data are
c = 1e-4*c_scaled
a = a_scaled*10^(-4-6*b_scaled)
b = b_scaled
where c_scaled, a_scaled and b_scaled are the coefficients to fit AC/1e-4 against w/1e6.
You can see this if you take
AC/1e-4 = c_scaled + a_scaled*(w/1e6)^b_scaled
in the form
AC = c + a*w^b
Answers (1)
Mathieu NOE
on 10 May 2022
hello
my suggestion - and because I don't have the Curve Fitting Tbx, simply with fminsearch
i prefered to look at the data in log log scale (btw , w is log spaced) and it's also more appropriate when dealing with power models (IMHO)
a = 3.4719e-14
n = 1.4408
dc = 3.3719e-06
clear;close all;clc;
opts = spreadsheetImportOptions("NumVariables", 6);
% Specify sheet and range
opts.Sheet = "Sheet1";
opts.DataRange = "A2:F202";
% Specify column names and types
opts.VariableNames = ["f1", "sigmaAcZT1", "f2", "SIGMAACZT2", "F2", "SIGMAACZT3"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double"];
% Import the data
% data = readtable("C:\Users\navaneeth\Downloads\IONIC CONDUCTIVITY.xlsx", opts, "UseExcel", false);
data = readtable("IONIC CONDUCTIVITY.xlsx", opts, "UseExcel", false);
% Convert to output type
data = table2array(data);
% Clear temporary variables
clear opts
w = 2*pi*data(:,1);%w = 2*pi*f
AC = data(:,2);% AC conductivity
% ft = fittype('c+a*w^b','dependent',{'Conductivity'},'independent',{'w'},'coefficients',{'a','b','c'});% AC_conductivity = Dc_conductivity +A*w^n
% f = fit(w,AC,ft,'StartPoint',[0,1,1]);
% plot(f,w,AC)
%% 3 parameters fminsearch optimization
f = @(a,n,dc,x) dc + a.*(x.^n) ; % AC_conductivity = Dc_conductivity +A*w^n
obj_fun = @(params) norm(f(params(1), params(2), params(3),w)-AC);
sol = fminsearch(obj_fun, [AC(end)/w(end),1,AC(1)]);
a = sol(1);
n = sol(2)
dc = sol(3);
ACfit = f(a,n,dc, w);
Rsquared = my_Rsquared_coeff(AC,ACfit); % correlation coefficient
figure(1);
loglog(w, AC, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
loglog(w, ACfit, '-');
title(['Data fit - R squared = ' num2str(Rsquared)]);
function Rsquared = my_Rsquared_coeff(data,data_fit)
% R2 correlation coefficient computation
% The total sum of squares
sum_of_squares = sum((data-mean(data)).^2);
% The sum of squares of residuals, also called the residual sum of squares:
sum_of_squares_of_residuals = sum((data-data_fit).^2);
% definition of the coefficient of correlation is
Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
end
See Also
Categories
Find more on Linear and Nonlinear Regression 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!