Difference in computational time
Hi all, I am working on a curve fitting project and I have a couple of curves to fit. For 2 of the curves, the fitting function involves an integration and I was told that I will need to make each of the data loop through the function (using the for loop). The below is my first code, which takes less than a second to complete.
tic
% Preparation
clear; clc
load Steady_State_Data.mat
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Clean up of data to select range of values
Absorbance = log10(T_substrate./T_sample);
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero
n = length(E) + 1;
% Data for fitting
E_p = E(n-N:n-1);
Abs = Absorbance(n-N:n-1) – min(Absorbance);
% Fitting function: Elliott’s Model for Steady State with non-parabolic factor
function F = EM_SS_wR(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*…
(integral(@(E)sech(((E_p – E)./p(6)))*(1 + 10*p(5)*(E – p(3)) + …
126*(p(5))^2*(E – p(3))^2)/(1 – exp(-2*pi*sqrt(p(4)/(E – p(3))))), p(3), Inf, ‘ArrayValued’, 1)) + …
p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(…
(1/1^3)*sech((E_p – p(3) + p(4)/1^2)./p(6)) + …
(1/2^3)*sech((E_p – p(3) + p(4)/2^2)./p(6)) + …
(1/3^3)*sech((E_p – p(3) + p(4)/3^2)./p(6)) + …
(1/4^3)*sech((E_p – p(3) + p(4)/4^2)./p(6)) + …
(1/5^3)*sech((E_p – p(3) + p(4)/5^2)./p(6)) + …
(1/6^3)*sech((E_p – p(3) + p(4)/6^2)./p(6)) + …
(1/7^3)*sech((E_p – p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
% Carrier Contribution
function F = EM_SS_wR_CC(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*…
(integral(@(E)sech(((E_p – E)./p(6)))*(1 + 10*p(5)*(E – p(3)) + …
126*(p(5))^2*(E – p(3))^2)/(1 – exp(-2*pi*sqrt(p(4)/(E – p(3))))), p(3), Inf, ‘ArrayValued’, 1));
end
F = F(:);
end
% Excitoninc Contribution
EM_SS_wR_EC = @(p, E_p) p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(…
(1/1^3)*sech((E_p – p(3) + p(4)/1^2)./p(6)) + …
(1/2^3)*sech((E_p – p(3) + p(4)/2^2)./p(6)) + …
(1/3^3)*sech((E_p – p(3) + p(4)/3^2)./p(6)) + …
(1/4^3)*sech((E_p – p(3) + p(4)/4^2)./p(6)) + …
(1/5^3)*sech((E_p – p(3) + p(4)/5^2)./p(6)) + …
(1/6^3)*sech((E_p – p(3) + p(4)/6^2)./p(6)) + …
(1/7^3)*sech((E_p – p(3) + p(4)/7^2)./p(6)));
%% Curve fitting options
% Initial parameter guess and bounds
lb = [-Inf, -Inf, 0, 0, -Inf, -Inf]; ub = [Inf, Inf, 55, 0.030, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘levenberg-marquardt’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, 1000, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘trust-region-reflective’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘interior-point’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, 1000, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% Solver for lsqcurvefit
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS_wR, p0, E_p, Abs, lb, ub);
%% Error bars calculation
ci = nlparci(p, residual, ‘Jacobian’, jacobian);
PCI = table(ci(:,1), p(:), ci(:,2),’VariableNames’,{‘CI 0.025′,’p’,’CI 0.975′});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) – Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) – Parameter_CI(:,1))./2;
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
e1 = lower_bar(1,1);
e2 = lower_bar(2,1);
e3 = lower_bar(3,1);
e4 = lower_bar(4,1);
e5 = lower_bar(5,1);
e6 = lower_bar(6,1);
X1 = [‘ A1 = ‘, num2str(p1), char(177), num2str(e1)];
X2 = [‘ A2 = ‘, num2str(p2), char(177), num2str(e2)];
X3 = [‘ Eg = ‘, num2str(p3), char(177), num2str(e3)];
X4 = [‘ Eb = ‘, num2str(p4), char(177), num2str(e4)];
X5 = [‘ R = ‘, num2str(p5), char(177), num2str(e5)];
X6 = [‘ g = ‘, num2str(p6), char(177), num2str(e6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
% Table of parameter values
parameter_name = {‘A1’; ‘A2’; ‘Eg’; ‘Eb’; ‘R’; ‘g’};
parameter_values = [p1; p2; p3; p4; p5; p6];
error = [e1; e2; e3; e4; e5; e6];
T = table(parameter_values, error, ‘RowNames’,parameter_name);
disp(T)
%% Plot command
plot(E_p, Abs, ‘o’, ‘Color’,’b’) % scatter plot
hold on
plot(E_p, EM_SS_wR(p, E_p), ‘LineWidth’, 1.0, ‘Color’,’red’) % curve fitting
plot(E_p, EM_SS_wR_CC(p, E_p), ‘LineWidth’, 1.0,’Color’,’black’) % carrier contribution
plot(E_p, EM_SS_wR_EC(p, E_p), ‘LineWidth’, 1.0, ‘Color’,’green’) % excitonic contribution
% Table of parameter values
TString = evalc(‘disp(T)’);
TString = strrep(TString,'<strong>’,’bf’);
TString = strrep(TString,'</strong>’,’rm’);
TString = strrep(TString,’_’,’_’);
FixedWidth = get(0,’FixedWidthFontName’);
annotation(gcf,’Textbox’,’String’,TString,’Interpreter’,’tex’,…
‘FontName’,FixedWidth,’Units’,’Normalized’,’Position’,[0.15, 0.8, 0.1, 0.1]);
xlabel(‘Probe energy (eV)’)
ylabel(‘Absorbance.O.D’)
legend(‘Experimental Data’, ‘Fitted Curve’, ‘Carrier Contribution’,’Excitonic Contribution’, ‘Location’,’southeast’)
hold off
%% Miscellaneous
% For a rectangular pic:
% drag horizontally until it just covers the letter ‘c’ in the word ‘col’
% For a square pic:
% drag horizontall until it just covers the letter ‘o’ in the word ‘Contribution’ in the next line)
% legend(‘Experimental Data’, ‘Fitted Curve’, ‘Carrier Contribution’,’Excitonic Contribution’)
toc
However, when using my second code (shown below) it takes around 5 minutes to complete running. May I please ask for the difference in the computational time?
tic
%% Preparation
clc; clear
data = importdata("FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis
load Steady_State_Parameter_Values.mat
%% Preamble
% Fundamental constants
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
% Data
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
delay_t = data(1, :);
% Parameter values from Steady state fiting
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
t1 = 0.5342;
t2 = 1.0257;
t3 = 1.9673;
t4 = 3.9554;
carrier_T = [1198.8, 816.7, 446.8, 328.7];
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
% Data for fitting
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
delta_Abs1 = data(Range_W,col1);
delta_Abs2 = data(Range_W,col2);
delta_Abs3 = data(Range_W,col3);
delta_Abs4 = data(Range_W,col4);
% Fitting function: Elliott’s Model for Transient Absorption
function F = EM_TA_wR1(x, e_p)
load Steady_State_Parameter_Values.mat
kB = 8.617333268*10^-5; % units: eV/ K
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 1;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)./(E_p).*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)./E_p.*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR2(x, e_p)
data = importdata("Experimental dataTransient AbsorptionFCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 2;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR3(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 3;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR4(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 4;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
% Curve fitting
lb = [Eg, Eb, g, 0.3, 0.5, 0.2]; ub = [55, 0.05, 0.05, 20, 2, 1];
x0 = [1.65, 0.03, 0.03, 1.3, 1, 0.3];
optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘levenberg-marquardt’, ‘MaxFunctionEvaluations’,10^5, ‘MaxIterations’, 10^5, ‘FunctionTolerance’,10^-20, ‘StepTolerance’, 10^-20);
delta_Abs(:,1)= data_new2(:,1);
x1 = lsqcurvefit(@EM_TA_wR1, x0, E_p, delta_Abs1, lb, ub, optim_lsq);
x2 = lsqcurvefit(@EM_TA_wR2, x1, E_p, delta_Abs2, lb, ub, optim_lsq);
x3 = lsqcurvefit(@EM_TA_wR3, x2, E_p, delta_Abs3, lb, ub, optim_lsq);
x4 = lsqcurvefit(@EM_TA_wR4, x3, E_p, delta_Abs4, lb, ub, optim_lsq);
plot(E_p, delta_Abs1, ‘o’,’color’,’blue’)
hold on
plot(E_p, delta_Abs2, ‘o’, ‘color’, ‘yellow’)
plot(E_p, delta_Abs3, ‘o’, ‘color’, ‘green’)
plot(E_p, delta_Abs4, ‘o’, ‘color’, ‘magenta’)
plot(E_p, EM_TA_wR1(x1, E_p), ‘color’, ‘blue’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR2(x2, E_p), ‘color’, ‘yellow’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR3(x3, E_p), ‘color’, ‘green’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR4(x4, E_p), ‘color’, ‘magenta’, ‘LineWidth’, 4.0)
xlabel(‘Probe Photon Energy (eV)’)
ylabel(‘Delta A (O.D.)’)
legend(‘0.5 ps’, ‘1.0 ps’, ‘2.0 ps’, ‘4.0 ps’, ‘Location’, ‘southeast’)
hold off
disp(x1)
disp(x2)
disp(x3)
disp(x4)
tocHi all, I am working on a curve fitting project and I have a couple of curves to fit. For 2 of the curves, the fitting function involves an integration and I was told that I will need to make each of the data loop through the function (using the for loop). The below is my first code, which takes less than a second to complete.
tic
% Preparation
clear; clc
load Steady_State_Data.mat
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Clean up of data to select range of values
Absorbance = log10(T_substrate./T_sample);
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero
n = length(E) + 1;
% Data for fitting
E_p = E(n-N:n-1);
Abs = Absorbance(n-N:n-1) – min(Absorbance);
% Fitting function: Elliott’s Model for Steady State with non-parabolic factor
function F = EM_SS_wR(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*…
(integral(@(E)sech(((E_p – E)./p(6)))*(1 + 10*p(5)*(E – p(3)) + …
126*(p(5))^2*(E – p(3))^2)/(1 – exp(-2*pi*sqrt(p(4)/(E – p(3))))), p(3), Inf, ‘ArrayValued’, 1)) + …
p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(…
(1/1^3)*sech((E_p – p(3) + p(4)/1^2)./p(6)) + …
(1/2^3)*sech((E_p – p(3) + p(4)/2^2)./p(6)) + …
(1/3^3)*sech((E_p – p(3) + p(4)/3^2)./p(6)) + …
(1/4^3)*sech((E_p – p(3) + p(4)/4^2)./p(6)) + …
(1/5^3)*sech((E_p – p(3) + p(4)/5^2)./p(6)) + …
(1/6^3)*sech((E_p – p(3) + p(4)/6^2)./p(6)) + …
(1/7^3)*sech((E_p – p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
% Carrier Contribution
function F = EM_SS_wR_CC(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*…
(integral(@(E)sech(((E_p – E)./p(6)))*(1 + 10*p(5)*(E – p(3)) + …
126*(p(5))^2*(E – p(3))^2)/(1 – exp(-2*pi*sqrt(p(4)/(E – p(3))))), p(3), Inf, ‘ArrayValued’, 1));
end
F = F(:);
end
% Excitoninc Contribution
EM_SS_wR_EC = @(p, E_p) p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(…
(1/1^3)*sech((E_p – p(3) + p(4)/1^2)./p(6)) + …
(1/2^3)*sech((E_p – p(3) + p(4)/2^2)./p(6)) + …
(1/3^3)*sech((E_p – p(3) + p(4)/3^2)./p(6)) + …
(1/4^3)*sech((E_p – p(3) + p(4)/4^2)./p(6)) + …
(1/5^3)*sech((E_p – p(3) + p(4)/5^2)./p(6)) + …
(1/6^3)*sech((E_p – p(3) + p(4)/6^2)./p(6)) + …
(1/7^3)*sech((E_p – p(3) + p(4)/7^2)./p(6)));
%% Curve fitting options
% Initial parameter guess and bounds
lb = [-Inf, -Inf, 0, 0, -Inf, -Inf]; ub = [Inf, Inf, 55, 0.030, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘levenberg-marquardt’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, 1000, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘trust-region-reflective’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘interior-point’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, 1000, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% Solver for lsqcurvefit
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS_wR, p0, E_p, Abs, lb, ub);
%% Error bars calculation
ci = nlparci(p, residual, ‘Jacobian’, jacobian);
PCI = table(ci(:,1), p(:), ci(:,2),’VariableNames’,{‘CI 0.025′,’p’,’CI 0.975′});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) – Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) – Parameter_CI(:,1))./2;
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
e1 = lower_bar(1,1);
e2 = lower_bar(2,1);
e3 = lower_bar(3,1);
e4 = lower_bar(4,1);
e5 = lower_bar(5,1);
e6 = lower_bar(6,1);
X1 = [‘ A1 = ‘, num2str(p1), char(177), num2str(e1)];
X2 = [‘ A2 = ‘, num2str(p2), char(177), num2str(e2)];
X3 = [‘ Eg = ‘, num2str(p3), char(177), num2str(e3)];
X4 = [‘ Eb = ‘, num2str(p4), char(177), num2str(e4)];
X5 = [‘ R = ‘, num2str(p5), char(177), num2str(e5)];
X6 = [‘ g = ‘, num2str(p6), char(177), num2str(e6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
% Table of parameter values
parameter_name = {‘A1’; ‘A2’; ‘Eg’; ‘Eb’; ‘R’; ‘g’};
parameter_values = [p1; p2; p3; p4; p5; p6];
error = [e1; e2; e3; e4; e5; e6];
T = table(parameter_values, error, ‘RowNames’,parameter_name);
disp(T)
%% Plot command
plot(E_p, Abs, ‘o’, ‘Color’,’b’) % scatter plot
hold on
plot(E_p, EM_SS_wR(p, E_p), ‘LineWidth’, 1.0, ‘Color’,’red’) % curve fitting
plot(E_p, EM_SS_wR_CC(p, E_p), ‘LineWidth’, 1.0,’Color’,’black’) % carrier contribution
plot(E_p, EM_SS_wR_EC(p, E_p), ‘LineWidth’, 1.0, ‘Color’,’green’) % excitonic contribution
% Table of parameter values
TString = evalc(‘disp(T)’);
TString = strrep(TString,'<strong>’,’bf’);
TString = strrep(TString,'</strong>’,’rm’);
TString = strrep(TString,’_’,’_’);
FixedWidth = get(0,’FixedWidthFontName’);
annotation(gcf,’Textbox’,’String’,TString,’Interpreter’,’tex’,…
‘FontName’,FixedWidth,’Units’,’Normalized’,’Position’,[0.15, 0.8, 0.1, 0.1]);
xlabel(‘Probe energy (eV)’)
ylabel(‘Absorbance.O.D’)
legend(‘Experimental Data’, ‘Fitted Curve’, ‘Carrier Contribution’,’Excitonic Contribution’, ‘Location’,’southeast’)
hold off
%% Miscellaneous
% For a rectangular pic:
% drag horizontally until it just covers the letter ‘c’ in the word ‘col’
% For a square pic:
% drag horizontall until it just covers the letter ‘o’ in the word ‘Contribution’ in the next line)
% legend(‘Experimental Data’, ‘Fitted Curve’, ‘Carrier Contribution’,’Excitonic Contribution’)
toc
However, when using my second code (shown below) it takes around 5 minutes to complete running. May I please ask for the difference in the computational time?
tic
%% Preparation
clc; clear
data = importdata("FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis
load Steady_State_Parameter_Values.mat
%% Preamble
% Fundamental constants
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
% Data
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
delay_t = data(1, :);
% Parameter values from Steady state fiting
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
t1 = 0.5342;
t2 = 1.0257;
t3 = 1.9673;
t4 = 3.9554;
carrier_T = [1198.8, 816.7, 446.8, 328.7];
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
% Data for fitting
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
delta_Abs1 = data(Range_W,col1);
delta_Abs2 = data(Range_W,col2);
delta_Abs3 = data(Range_W,col3);
delta_Abs4 = data(Range_W,col4);
% Fitting function: Elliott’s Model for Transient Absorption
function F = EM_TA_wR1(x, e_p)
load Steady_State_Parameter_Values.mat
kB = 8.617333268*10^-5; % units: eV/ K
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 1;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)./(E_p).*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)./E_p.*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR2(x, e_p)
data = importdata("Experimental dataTransient AbsorptionFCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 2;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR3(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 3;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR4(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 4;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
% Curve fitting
lb = [Eg, Eb, g, 0.3, 0.5, 0.2]; ub = [55, 0.05, 0.05, 20, 2, 1];
x0 = [1.65, 0.03, 0.03, 1.3, 1, 0.3];
optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘levenberg-marquardt’, ‘MaxFunctionEvaluations’,10^5, ‘MaxIterations’, 10^5, ‘FunctionTolerance’,10^-20, ‘StepTolerance’, 10^-20);
delta_Abs(:,1)= data_new2(:,1);
x1 = lsqcurvefit(@EM_TA_wR1, x0, E_p, delta_Abs1, lb, ub, optim_lsq);
x2 = lsqcurvefit(@EM_TA_wR2, x1, E_p, delta_Abs2, lb, ub, optim_lsq);
x3 = lsqcurvefit(@EM_TA_wR3, x2, E_p, delta_Abs3, lb, ub, optim_lsq);
x4 = lsqcurvefit(@EM_TA_wR4, x3, E_p, delta_Abs4, lb, ub, optim_lsq);
plot(E_p, delta_Abs1, ‘o’,’color’,’blue’)
hold on
plot(E_p, delta_Abs2, ‘o’, ‘color’, ‘yellow’)
plot(E_p, delta_Abs3, ‘o’, ‘color’, ‘green’)
plot(E_p, delta_Abs4, ‘o’, ‘color’, ‘magenta’)
plot(E_p, EM_TA_wR1(x1, E_p), ‘color’, ‘blue’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR2(x2, E_p), ‘color’, ‘yellow’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR3(x3, E_p), ‘color’, ‘green’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR4(x4, E_p), ‘color’, ‘magenta’, ‘LineWidth’, 4.0)
xlabel(‘Probe Photon Energy (eV)’)
ylabel(‘Delta A (O.D.)’)
legend(‘0.5 ps’, ‘1.0 ps’, ‘2.0 ps’, ‘4.0 ps’, ‘Location’, ‘southeast’)
hold off
disp(x1)
disp(x2)
disp(x3)
disp(x4)
toc Hi all, I am working on a curve fitting project and I have a couple of curves to fit. For 2 of the curves, the fitting function involves an integration and I was told that I will need to make each of the data loop through the function (using the for loop). The below is my first code, which takes less than a second to complete.
tic
% Preparation
clear; clc
load Steady_State_Data.mat
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Clean up of data to select range of values
Absorbance = log10(T_substrate./T_sample);
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero
n = length(E) + 1;
% Data for fitting
E_p = E(n-N:n-1);
Abs = Absorbance(n-N:n-1) – min(Absorbance);
% Fitting function: Elliott’s Model for Steady State with non-parabolic factor
function F = EM_SS_wR(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*…
(integral(@(E)sech(((E_p – E)./p(6)))*(1 + 10*p(5)*(E – p(3)) + …
126*(p(5))^2*(E – p(3))^2)/(1 – exp(-2*pi*sqrt(p(4)/(E – p(3))))), p(3), Inf, ‘ArrayValued’, 1)) + …
p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(…
(1/1^3)*sech((E_p – p(3) + p(4)/1^2)./p(6)) + …
(1/2^3)*sech((E_p – p(3) + p(4)/2^2)./p(6)) + …
(1/3^3)*sech((E_p – p(3) + p(4)/3^2)./p(6)) + …
(1/4^3)*sech((E_p – p(3) + p(4)/4^2)./p(6)) + …
(1/5^3)*sech((E_p – p(3) + p(4)/5^2)./p(6)) + …
(1/6^3)*sech((E_p – p(3) + p(4)/6^2)./p(6)) + …
(1/7^3)*sech((E_p – p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
% Carrier Contribution
function F = EM_SS_wR_CC(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*…
(integral(@(E)sech(((E_p – E)./p(6)))*(1 + 10*p(5)*(E – p(3)) + …
126*(p(5))^2*(E – p(3))^2)/(1 – exp(-2*pi*sqrt(p(4)/(E – p(3))))), p(3), Inf, ‘ArrayValued’, 1));
end
F = F(:);
end
% Excitoninc Contribution
EM_SS_wR_EC = @(p, E_p) p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(…
(1/1^3)*sech((E_p – p(3) + p(4)/1^2)./p(6)) + …
(1/2^3)*sech((E_p – p(3) + p(4)/2^2)./p(6)) + …
(1/3^3)*sech((E_p – p(3) + p(4)/3^2)./p(6)) + …
(1/4^3)*sech((E_p – p(3) + p(4)/4^2)./p(6)) + …
(1/5^3)*sech((E_p – p(3) + p(4)/5^2)./p(6)) + …
(1/6^3)*sech((E_p – p(3) + p(4)/6^2)./p(6)) + …
(1/7^3)*sech((E_p – p(3) + p(4)/7^2)./p(6)));
%% Curve fitting options
% Initial parameter guess and bounds
lb = [-Inf, -Inf, 0, 0, -Inf, -Inf]; ub = [Inf, Inf, 55, 0.030, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘levenberg-marquardt’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, 1000, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘trust-region-reflective’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘interior-point’, ‘MaxFunctionEvaluations’,1000, ‘MaxIterations’, 1000, ‘FunctionTolerance’,10^-9, ‘StepTolerance’, 10^-9);
% Solver for lsqcurvefit
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS_wR, p0, E_p, Abs, lb, ub);
%% Error bars calculation
ci = nlparci(p, residual, ‘Jacobian’, jacobian);
PCI = table(ci(:,1), p(:), ci(:,2),’VariableNames’,{‘CI 0.025′,’p’,’CI 0.975′});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) – Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) – Parameter_CI(:,1))./2;
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
e1 = lower_bar(1,1);
e2 = lower_bar(2,1);
e3 = lower_bar(3,1);
e4 = lower_bar(4,1);
e5 = lower_bar(5,1);
e6 = lower_bar(6,1);
X1 = [‘ A1 = ‘, num2str(p1), char(177), num2str(e1)];
X2 = [‘ A2 = ‘, num2str(p2), char(177), num2str(e2)];
X3 = [‘ Eg = ‘, num2str(p3), char(177), num2str(e3)];
X4 = [‘ Eb = ‘, num2str(p4), char(177), num2str(e4)];
X5 = [‘ R = ‘, num2str(p5), char(177), num2str(e5)];
X6 = [‘ g = ‘, num2str(p6), char(177), num2str(e6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
% Table of parameter values
parameter_name = {‘A1’; ‘A2’; ‘Eg’; ‘Eb’; ‘R’; ‘g’};
parameter_values = [p1; p2; p3; p4; p5; p6];
error = [e1; e2; e3; e4; e5; e6];
T = table(parameter_values, error, ‘RowNames’,parameter_name);
disp(T)
%% Plot command
plot(E_p, Abs, ‘o’, ‘Color’,’b’) % scatter plot
hold on
plot(E_p, EM_SS_wR(p, E_p), ‘LineWidth’, 1.0, ‘Color’,’red’) % curve fitting
plot(E_p, EM_SS_wR_CC(p, E_p), ‘LineWidth’, 1.0,’Color’,’black’) % carrier contribution
plot(E_p, EM_SS_wR_EC(p, E_p), ‘LineWidth’, 1.0, ‘Color’,’green’) % excitonic contribution
% Table of parameter values
TString = evalc(‘disp(T)’);
TString = strrep(TString,'<strong>’,’bf’);
TString = strrep(TString,'</strong>’,’rm’);
TString = strrep(TString,’_’,’_’);
FixedWidth = get(0,’FixedWidthFontName’);
annotation(gcf,’Textbox’,’String’,TString,’Interpreter’,’tex’,…
‘FontName’,FixedWidth,’Units’,’Normalized’,’Position’,[0.15, 0.8, 0.1, 0.1]);
xlabel(‘Probe energy (eV)’)
ylabel(‘Absorbance.O.D’)
legend(‘Experimental Data’, ‘Fitted Curve’, ‘Carrier Contribution’,’Excitonic Contribution’, ‘Location’,’southeast’)
hold off
%% Miscellaneous
% For a rectangular pic:
% drag horizontally until it just covers the letter ‘c’ in the word ‘col’
% For a square pic:
% drag horizontall until it just covers the letter ‘o’ in the word ‘Contribution’ in the next line)
% legend(‘Experimental Data’, ‘Fitted Curve’, ‘Carrier Contribution’,’Excitonic Contribution’)
toc
However, when using my second code (shown below) it takes around 5 minutes to complete running. May I please ask for the difference in the computational time?
tic
%% Preparation
clc; clear
data = importdata("FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis
load Steady_State_Parameter_Values.mat
%% Preamble
% Fundamental constants
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
% Data
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
delay_t = data(1, :);
% Parameter values from Steady state fiting
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
t1 = 0.5342;
t2 = 1.0257;
t3 = 1.9673;
t4 = 3.9554;
carrier_T = [1198.8, 816.7, 446.8, 328.7];
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
% Data for fitting
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
delta_Abs1 = data(Range_W,col1);
delta_Abs2 = data(Range_W,col2);
delta_Abs3 = data(Range_W,col3);
delta_Abs4 = data(Range_W,col4);
% Fitting function: Elliott’s Model for Transient Absorption
function F = EM_TA_wR1(x, e_p)
load Steady_State_Parameter_Values.mat
kB = 8.617333268*10^-5; % units: eV/ K
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 1;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)./(E_p).*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)./E_p.*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR2(x, e_p)
data = importdata("Experimental dataTransient AbsorptionFCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 2;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR3(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 3;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
function F = EM_TA_wR4(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 4;
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p – E)./g))*(1 + 10*R*(E – x(1)) + …
126*(R)^2*(E – x(1))^2)/(1 – exp(-2*pi*sqrt(x(2)/(E – x(1))))), x(1), Inf, ‘ArrayValued’, 1))*(1 – 1/(1 + exp((E_p – x(4))/(kB*carrier_T(1, n)))))^2 – …
A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*…
(integral(@(e)sech(((E_p – e)./g))*(1 + 10*R*(e – Eg) + …
126*(R)^2*(e – Eg)^2)/(1 – exp(-2*pi*sqrt(Eb/(e – Eg)))), Eg, Inf, ‘ArrayValued’, 1)) + …
x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(…
(1/1^3)*sech((E_p – x(1) + x(2)/1^2)./x(3)) + …
(1/2^3)*sech((E_p – x(1) + x(2)/2^2)./x(3)) + …
(1/3^3)*sech((E_p – x(1) + x(2)/3^2)./x(3)) + …
(1/4^3)*sech((E_p – x(1) + x(2)/4^2)./x(3)) + …
(1/5^3)*sech((E_p – x(1) + x(2)/5^2)./x(3)) + …
(1/6^3)*sech((E_p – x(1) + x(2)/6^2)./x(3)) + …
(1/7^3)*sech((E_p – x(1) + x(2)/7^2)./x(3))) – …
A2*(4*pi*Eb^3/2)*1/g*(…
(1/1^3)*sech((E_p – Eg + Eb/1^2)./g) + …
(1/2^3)*sech((E_p – Eg + Eb/2^2)./g) + …
(1/3^3)*sech((E_p – Eg + Eb/3^2)./g) + …
(1/4^3)*sech((E_p – Eg + Eb/4^2)./g) + …
(1/5^3)*sech((E_p – Eg + Eb/5^2)./g) + …
(1/6^3)*sech((E_p – Eg + Eb/6^2)./g) + …
(1/7^3)*sech((E_p – Eg + Eb/7^2)./g));
end
F = F(:);
end
% Curve fitting
lb = [Eg, Eb, g, 0.3, 0.5, 0.2]; ub = [55, 0.05, 0.05, 20, 2, 1];
x0 = [1.65, 0.03, 0.03, 1.3, 1, 0.3];
optim_lsq = optimoptions(‘lsqcurvefit’, ‘Algorithm’, ‘levenberg-marquardt’, ‘MaxFunctionEvaluations’,10^5, ‘MaxIterations’, 10^5, ‘FunctionTolerance’,10^-20, ‘StepTolerance’, 10^-20);
delta_Abs(:,1)= data_new2(:,1);
x1 = lsqcurvefit(@EM_TA_wR1, x0, E_p, delta_Abs1, lb, ub, optim_lsq);
x2 = lsqcurvefit(@EM_TA_wR2, x1, E_p, delta_Abs2, lb, ub, optim_lsq);
x3 = lsqcurvefit(@EM_TA_wR3, x2, E_p, delta_Abs3, lb, ub, optim_lsq);
x4 = lsqcurvefit(@EM_TA_wR4, x3, E_p, delta_Abs4, lb, ub, optim_lsq);
plot(E_p, delta_Abs1, ‘o’,’color’,’blue’)
hold on
plot(E_p, delta_Abs2, ‘o’, ‘color’, ‘yellow’)
plot(E_p, delta_Abs3, ‘o’, ‘color’, ‘green’)
plot(E_p, delta_Abs4, ‘o’, ‘color’, ‘magenta’)
plot(E_p, EM_TA_wR1(x1, E_p), ‘color’, ‘blue’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR2(x2, E_p), ‘color’, ‘yellow’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR3(x3, E_p), ‘color’, ‘green’, ‘LineWidth’, 4.0)
plot(E_p, EM_TA_wR4(x4, E_p), ‘color’, ‘magenta’, ‘LineWidth’, 4.0)
xlabel(‘Probe Photon Energy (eV)’)
ylabel(‘Delta A (O.D.)’)
legend(‘0.5 ps’, ‘1.0 ps’, ‘2.0 ps’, ‘4.0 ps’, ‘Location’, ‘southeast’)
hold off
disp(x1)
disp(x2)
disp(x3)
disp(x4)
toc curve fitting, computational time MATLAB Answers — New Questions
