Curve fitting for jonscher power law

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)

[Torsten]

Change the units. Scale w by 1e6, scale conductivity by 1e-4 before the fitting process.

[Navaneeth Tejasvi Mysore Nagendra]

Thank you it works, though the scalling changes. I don't know this is what you suggested?

change:

clear opts

f = data(:,1);

w = 2*pi*data(:,1)./1e6;%w = 2*pi*f

AC = data(:,2)./1e-4;% 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)

do you know how can I fix the scaling on the plot to its original powers,

and in meanwhile I used power2 fittype with Levenberg - Marquardt algorithm to fit the data and based on that I had got the values for my coefficients as,

program,

function [fitresult, gof] = createFit(w, AC1)

%CREATEFIT(W,AC)

% Create a fit.

%

% Data for 'untitled fit 1' fit:

% X Input : w

% Y Output: AC

% Output:

% fitresult : a fit object representing the fit.

% gof : structure with goodness-of fit info.

%

% See also FIT, CFIT, SFIT.

% Auto-generated by MATLAB on 10-May-2022 11:38:11

%% Fit: 'untitled fit 1'.

clc;clear,close all;

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

f = data(:,1)

w = 2*pi*data(:,1);%w = 2*pi*f

AC1 = data(:,2);%AC conductivity sample 1

AC2 = data(:,4);%AC conductivity sample 2

AC3 = data(:,6);%AC conductivity sample 3

[xData1, yData1] = prepareCurveData( w, AC1 );

[xData2, yData2] = prepareCurveData( w, AC2 );

[xData3, yData3] = prepareCurveData( w, AC3 );

% Set up fittype and options.

ft = fittype( 'power2' );

opts = fitoptions( 'Method', 'NonlinearLeastSquares' );

opts.Algorithm = 'Levenberg-Marquardt';

opts.Display = 'Off';

opts.StartPoint = [9.32271357188677e-06 0.0429375187609614 2.50935751249309e-06];

% Fit model to data.

[fitresult1, gof1] = fit( xData1, yData1, ft, opts );

[fitresult2, gof2] = fit( xData2, yData2, ft, opts );

[fitresult3, gof3] = fit( xData3, yData3, ft, opts );

% Plot fit with data.

figure( 'Name', 'untitled fit 1' );

h = plot(fitresult1,'k',xData1,yData1,'ok');

hold on

g = plot(fitresult2,'k',xData2,yData2,'or');

hold on

f = plot(fitresult3,'k',xData3,yData3,'ob');

legend('ZT-3','fit', 'Location', 'NorthWest')

set(h, 'LineWidth',2)

set(g, 'LineWidth',2)

set(f, 'LineWidth',2)

% Label axes

set(gca,'fontweight','bold','fontsize',16);

xlabel( '\omega');

ylabel('\sigma_{AC} (S/cm)');

fitresult1

fitresult2

fitresult3

grid on

for which,

fitresult1 =

General model Power2:

fitresult1(x) = a*x^b+c

Coefficients (with 95% confidence bounds):

a = 2.326e-14 (1.972e-14, 2.68e-14)

b = 1.467 (1.457, 1.477)

c = 3.475e-06 (3.335e-06, 3.615e-06)

but these values changes from the first fit I did (apart from the power),

val =

General model:

val(w) = c+a*w^b

Coefficients (with 95% confidence bounds):

a = 0.1476 (0.1452, 0.1499)

b = 1.467 (1.457, 1.477)

c = 0.03475 (0.03335, 0.03615)

I dont know why

[Torsten]

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

[Navaneeth Tejasvi Mysore Nagendra]

Understood thank you

Answer 1

[Mathieu NOE]

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

[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.