Curve fitting for jonscher power law

7 views (last 30 days)
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
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

Sign in to comment.

Answers (1)

Mathieu NOE
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
  1 Comment
Navaneeth Tejasvi Mysore Nagendra
Thank you for the detailed explanation, yes I need to plot it in log scale but struggling to do so.
This makes it easier to learn what I need to do.

Sign in to comment.

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!