I am coding a equation for curve fitting, but I am facing inaccuracy. I am not getting Elastic plateau in the curve. Can anyone guide me how I can address the issue.
Show older comments
clc; close all;
%% ---- Load data --
file = 'data.xlsx';
T = readtable(file); % expects columns: Strain, Stress (in GPa)
disp(T.Properties.VariableNames)
epsl = T.Strain_mm_mm_(:); % ensure column vector
sig = T.Stress_GPa_(:) * 1000; % convert GPa -> MPa
% keep only finite rows
good = isfinite(epsl) & isfinite(sig);
epsl = epsl(good);
sig = sig(good);
% (optional) trim extreme epsilon to avoid 1/(1-eps) blow-up in model
epsl = min(epsl, 0.981);
%% ---- Known strain rates -------------------------------------------------
edot = 0.000617; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
%% ---- Initial guesses & bounds ------------------------------------------
% [A, E, B, m, n, a, b]
p0 = [0.4, 2000, 0.3, 0.95, 2.0, -0.05, 0.150]; % initial params (MPa-based)
lb = [-0.1, 1.2, -.01, -1, 1, 1, 0.01 ];
ub = [2, 20, 1, 5, 5, 4.0, 0.5 ];
%% ---- Fit with lsqcurvefit (Optimization Toolbox) -----------------------
model = @(p,eps) constitutive_model(eps, edot, p(1), p(2), p(3), p(4), p(5), p(6), p(7), edot0);
try
opts = optimoptions('lsqcurvefit','Display','iter','MaxIterations',1000);
[pfit, ~] = lsqcurvefit(model, p0, epsl, sig, lb, ub, opts);
catch ME
warning('lsqcurvefit unavailable (%s). Falling back to fminsearch (unconstrained).', ME.identifier);
obj = @(p) sum((model(p,epsl) - sig).^2);
pfit = fminsearch(obj, p0);
end
sigma_fit = model(pfit, epsl);
%% ---- Goodness of fit (R^2, adj. R^2, RMSE) -----------------------------
SSE = sum((sig - sigma_fit).^2);
SST = sum((sig - mean(sig)).^2);
R2 = 1 - SSE/SST;
n = numel(sig);
k = numel(pfit); % number of fitted parameters
adjR2 = 1 - (1 - R2) * (n - 1) / max(n - k - 1, 1);
RMSE = sqrt(SSE / n);
% Print to console
fprintf('\nGoodness of fit:\n');
fprintf('R^2 = %.6f\n', R2)
fprintf('Adj R^2 = %.6f\n', adjR2);
fprintf('RMSE = %.6f MPa\n', RMSE)
%% ---- Plot and Save ------------------------------------------------------
fig = figure; hold on;
plot(epsl, sig, 'k-', 'MarkerFaceColor','k','MarkerSize',2, 'DisplayName','Experimental');
plot(epsl, sigma_fit,'r-', 'LineWidth',1, 'DisplayName','Fitted model');
xlabel('Strain \epsilon');
ylabel('Stress \sigma_d (MPa)');
title(['Fitted Constitutive Model at \epsilon-dot = ', num2str(edot), ' s^{-1}']);
legend show; grid on;
% Add a readable stats box
ax = gca;
xL = xlim(ax); yL = ylim(ax);
txt = {sprintf('R^2 = %.4f', R2), ...
sprintf('Adj R^2 = %.4f', adjR2), ...
sprintf('RMSE = %.3g MPa', RMSE)};
text(xL(1) + 0.02*range(xL), yL(2) - 0.04*range(yL), txt, ...
'VerticalAlignment','top', 'BackgroundColor','w', ...
'Margin',6, 'EdgeColor',[0.1 0.1 0.1], 'FontWeight','bold');
% Save the figure
saveas(fig, 'fitted_model_plot.png'); % PNG
% For crisper output you can also use:
% exportgraphics(fig,'fitted_model_plot.png','Resolution',300);
%% ---- Print parameters ---------------------------------------------------
fprintf('\nFitted Parameters:\n');
fprintf('A = %.4f MPa\n', pfit(1));
fprintf('E = %.4f MPa\n', pfit(2));
fprintf('B = %.4f MPa\n', pfit(3));
fprintf('m = %.4f\n', pfit(4));
fprintf('n = %.4f\n', pfit(5));
fprintf('a = %.4f\n', pfit(6));
fprintf('b = %.4f\n', pfit(7));
%% ===================== Local functions (MUST be at file end) =====================
function sigma_d = constitutive_model(epsilon, edot, A, E, B, m, n, a, b, edot0)
% Vector-safe guards
epsilon = min(max(epsilon, 0.001), 0.99); % keep in [0,1)
one_me = max(1 - epsilon, eps);
% Static part
term1 = A .* (1 - exp(-(E./A) .* epsilon .* (1 - epsilon))).^m;
term2 = B .* (epsilon ./ one_me).^n;
sigma_static = term1 + term2;
% Strain-rate factor
rate_factor = 1 + (a + b .* epsilon) .* log(edot ./ edot0);
sigma_d = sigma_static .* rate_factor;
end
Answers (1)
You need a better initial guess. Also, you should exclude the data from the linear elastic region, since clearly that is not captured by your model equation.
file = 'data.xlsx';
[epsl,sig]=readvars(file) ;
sig=sig*1000;% convert GPa -> MPa
% keep only finite and mid-range data
good = isfinite(epsl) & isfinite(sig) & epsl>0.15 & epsl<0.981;
epsl = epsl(good);
sig = sig(good);
%% ---- Known strain rates -------------------------------------------------
edot = 0.000617; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
%% ---- Initial guesses & bounds ------------------------------------------
% [A, E, B, m, n, a, b]
p0 = [0.04, 100, 0.03, 0.95, 2.0, -0.05, 0.150]; % initial params (MPa-based)
lb = [-0.1, 1.2, -.01, -1, 1, 1, 0.01 ];
ub = [2, 20, 1, 5, 5, 4.0, 0.5 ];
%% ---- Fit with lsqcurvefit (Optimization Toolbox) -----------------------
model = @(p,eps) constitutive_model(eps, edot, p(1), p(2), p(3), p(4), p(5), p(6), p(7), edot0);
try
opts = optimoptions('lsqcurvefit','Display','iter','MaxIterations',1000);
[pfit, ~] = lsqcurvefit(model, p0, epsl, sig, lb, ub, opts);
catch ME
warning('lsqcurvefit unavailable (%s). Falling back to fminsearch (unconstrained).', ME.identifier);
obj = @(p) sum((model(p,epsl) - sig).^2);
pfit = fminsearch(obj, p0);
end
sigma_fit = model(pfit, epsl);
postAnalysis(epsl,sig,pfit,sigma_fit,edot)
%% ===================== Local functions (MUST be at file end) =====================
function sigma_d = constitutive_model(epsilon, edot, A, E, B, m, n, a, b, edot0)
% Vector-safe guards
epsilon = min(max(epsilon, 0.001), 0.99); % keep in [0,1)
one_me = max(1 - epsilon, eps);
% Static part
term1 = A .* (1 - exp(-(E./A) .* epsilon .* (1 - epsilon))).^m;
term2 = B .* (epsilon ./ one_me).^n;
sigma_static = term1 + term2;
% Strain-rate factor
rate_factor = 1 + (a + b .* epsilon) .* log(edot ./ edot0);
sigma_d = sigma_static .* rate_factor;
end
function postAnalysis(epsl,sig,pfit, sigma_fit,edot)
%% ---- Goodness of fit (R^2, adj. R^2, RMSE) -----------------------------
SSE = sum((sig - sigma_fit).^2);
SST = sum((sig - mean(sig)).^2);
R2 = 1 - SSE/SST;
n = numel(sig);
k = numel(pfit); % number of fitted parameters
adjR2 = 1 - (1 - R2) * (n - 1) / max(n - k - 1, 1);
RMSE = sqrt(SSE / n);
% Print to console
fprintf('\nGoodness of fit:\n');
fprintf('R^2 = %.6f\n', R2)
fprintf('Adj R^2 = %.6f\n', adjR2);
fprintf('RMSE = %.6f MPa\n', RMSE)
%% ---- Plot and Save ------------------------------------------------------
fig = figure; hold on;
plot(epsl, sig, 'k-', 'MarkerFaceColor','k','MarkerSize',2, 'DisplayName','Experimental');
plot(epsl, sigma_fit,'r-', 'LineWidth',1, 'DisplayName','Fitted model');
xlabel('Strain \epsilon');
ylabel('Stress \sigma_d (MPa)');
title(['Fitted Constitutive Model at \epsilon-dot = ', num2str(edot), ' s^{-1}']);
legend show; grid on;
% Add a readable stats box
ax = gca;
xL = xlim(ax); yL = ylim(ax);
txt = {sprintf('R^2 = %.4f', R2), ...
sprintf('Adj R^2 = %.4f', adjR2), ...
sprintf('RMSE = %.3g MPa', RMSE)};
text(xL(1) + 0.02*range(xL), yL(2) - 0.04*range(yL), txt, ...
'VerticalAlignment','top', 'BackgroundColor','w', ...
'Margin',6, 'EdgeColor',[0.1 0.1 0.1], 'FontWeight','bold');
axis padded
end
14 Comments
JAMEEL
on 2 Sep 2025
Thank you for your time and guidance. But I need curve fitting according to provided data. Equation works well, I have taken that equation from another research paper. I just need to fit my data with model equation
But I need curve fitting according to provided data.
Is that different from what my answer proposes? Or do you mean you are not allowed to discard sections of the data, like the linear elastic part, from the fit?
If the latter, it seems futile. There doesn't appear to be any way that your equation can model a piecewise curve containing a linear elastic region at low epsl. If you think it can, show us a choice of parameters E,A,B,m,n from the paper that does so.
Why would it even be beneficial if you could? You know from Hooke's law that that part of the curve should be linear, so why try to approximate it with a nonlinear model?
JAMEEL
on 2 Sep 2025
Yes. Your it is different. I want to get elastic region too in the curve. I have to design another model.
I have to design another model.
Can't you just add a linear piece to the model? Example:
file = 'data.xlsx';
[epsl,sig]=readvars(file) ;
sig=sig*1000;% convert GPa -> MPa
% keep only data that is finite and substantially <1
good = isfinite(epsl) & isfinite(sig) & eps>0 & epsl<0.981;
epsl = epsl(good);
sig = sig(good);
%% ---- Known strain rates -------------------------------------------------
edot = 0.000617; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
%% ---- Initial guesses & bounds ------------------------------------------
% [A, E, B, m, n, a, b]
p0 = [0.04, 100, 0.03, 0.95, 2.0, -0.05, 0.150,0.15]; % initial params (MPa-based)
hooke=epsl<p0(end);
q=polyfit(epsl(hooke), sig(hooke),1);
p0=[p0,q];
lb = [-0.1, 1.2, -.01, -1, 1, 1, 0.01 ,min(epsl),0,-inf ];
ub = [2, 20, 1, 5, 5, 4.0, 0.5 ,inf ,inf,inf];
%% ---- Fit with lsqcurvefit (Optimization Toolbox) -----------------------
edotRatio=edot/edot0;
model = @(p,eps) constitutive_model(p, eps, edotRatio);
opts = optimoptions('lsqcurvefit','Display','iter','MaxIterations',1000);
[pfit, ~] = lsqcurvefit(model, p0, epsl, sig, lb, ub, opts)
sigma_fit = model(pfit, epsl);
postAnalysis(epsl,sig,pfit,sigma_fit,edot)
%% ===================== Local functions (MUST be at file end) =====================
function sigma_d = constitutive_model(p, epsilon, edotRatio)
[A, E, B, m, n, a, b, c]=deal(p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8));
elastic= (epsilon<=c);
sigmaHooke=p(9)*epsilon(elastic)+p(10);
epsilon=epsilon(~elastic);
one_me = max(1 - epsilon, eps);
% Static part
term1 = A .* (1 - exp(-(E./A) .* epsilon .* (1 - epsilon))).^m;
term2 = B .* (epsilon ./ one_me).^n;
sigma_static = term1 + term2;
% Strain-rate factor
rate_factor = 1 + (a + b .* epsilon) .* log(edotRatio);
sigmaInelastic = sigma_static .* rate_factor;
sigma_d=nan(size(elastic));
sigma_d(elastic)=sigmaHooke;
sigma_d(~elastic)=sigmaInelastic;
end
function postAnalysis(epsl,sig,pfit, sigma_fit,edot)
%% ---- Goodness of fit (R^2, adj. R^2, RMSE) -----------------------------
SSE = sum((sig - sigma_fit).^2);
SST = sum((sig - mean(sig)).^2);
R2 = 1 - SSE/SST;
n = numel(sig);
k = numel(pfit); % number of fitted parameters
adjR2 = 1 - (1 - R2) * (n - 1) / max(n - k - 1, 1);
RMSE = sqrt(SSE / n);
% Print to console
fprintf('\nGoodness of fit:\n');
fprintf('R^2 = %.6f\n', R2)
fprintf('Adj R^2 = %.6f\n', adjR2);
fprintf('RMSE = %.6f MPa\n', RMSE)
%% ---- Plot and Save ------------------------------------------------------
fig = figure; hold on;
plot(epsl, sig, 'k.','MarkerSize',2, 'DisplayName','Experimental');
plot(epsl, sigma_fit,'r-', 'LineWidth',1, 'DisplayName','Fitted model');
xlabel('Strain \epsilon');
ylabel('Stress \sigma_d (MPa)');
title(['Fitted Constitutive Model at \epsilon-dot = ', num2str(edot), ' s^{-1}']);
legend show; grid on;
% Add a readable stats box
ax = gca;
xL = xlim(ax); yL = ylim(ax);
txt = {sprintf('R^2 = %.4f', R2), ...
sprintf('Adj R^2 = %.4f', adjR2), ...
sprintf('RMSE = %.3g MPa', RMSE)};
text(xL(1) + 0.02*range(xL), yL(2) - 0.04*range(yL), txt, ...
'VerticalAlignment','top', 'BackgroundColor','w', ...
'Margin',6, 'EdgeColor',[0.1 0.1 0.1], 'FontWeight','bold');
axis padded
end
JAMEEL
on 2 Sep 2025
Matt J
on 2 Sep 2025
You're quite welcome, but since you seem satisfied with the answer, please Accept-click it.
clc
close all
clear all
file = 'PUF28_5mmmin.xlsx';
[epsl,sig]=readvars(file) ;
% sig=sig;% convert MPa -> MPa
% keep only data that is finite and substantially <1
good = isfinite(epsl) & isfinite(sig) & (epsl > 0) & (epsl < 0.981);
epsl = epsl(good);
sig = sig(good);
%% ---- Known strain rates -------------------------------------------------
edot = 0.00308; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
%% ---- Initial guesses & bounds ------------------------------------------
% [A, E, B, m, n, a, b]
p0 = [0.9,500, 0.08, 0.5, 2.0, 0.05, 50,0.07]; % initial params (MPa-based)
hooke=epsl<p0(end);
q=polyfit(epsl(hooke), sig(hooke),1);
p0=[p0,q];
lb = [-0.9, 50, -.01, -1, 1, 1, 0.01 ,min(epsl),0,-inf ];
ub = [5, 500, 2, 3, 2, 4.0, 0.1 ,inf ,inf,inf];
%% ---- Fit with lsqcurvefit (Optimization Toolbox) -----------------------
edotRatio=edot/edot0;
model = @(p,eps) constitutive_model(p, eps, edotRatio);
opts = optimoptions('lsqcurvefit','Display','iter','MaxIterations',1000);
[pfit, ~] = lsqcurvefit(model, p0, epsl, sig, lb, ub, opts)
sigma_fit = model(pfit, epsl);
postAnalysis(epsl,sig,pfit,sigma_fit,edot)
%% ===================== Local functions (MUST be at file end) =====================
function sigma_d = constitutive_model(p, epsilon, edotRatio)
% p = [A, E, B, m, n, a, b, c, slope, intercept]
[A, E, B, m, n, a, b, c] = deal(p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8));
elastic = (epsilon <= c);
sigmaHooke = p(9)*epsilon(elastic) + p(10);
e_in = epsilon(~elastic);
one_me = max(1 - e_in, eps);
% --- CORRECT term1 per paper: m is inside the exponential on (1-e)^m
term1 = A .* (1 - exp(-(E./A) .* e_in .* (one_me).^m));
term2 = B .* (e_in ./ one_me).^n;
sigma_static = term1 + term2;
rate_factor = 1 + (a + b .* e_in) .* log(edotRatio);
sigmaInelastic = sigma_static .* rate_factor;
sigma_d = nan(size(epsilon));
sigma_d(elastic) = sigmaHooke;
sigma_d(~elastic) = sigmaInelastic;
end
function postAnalysis(epsl,sig,pfit, sigma_fit,edot)
%% ---- Goodness of fit (R^2, adj. R^2, RMSE) -----------------------------
SSE = sum((sig - sigma_fit).^2);
SST = sum((sig - mean(sig)).^2);
R2 = 1 - SSE/SST;
n = numel(sig);
k = numel(pfit); % number of fitted parameters
adjR2 = 1 - (1 - R2) * (n - 1) / max(n - k - 1, 1);
RMSE = sqrt(SSE / n);
% Print to console
fprintf('\nGoodness of fit:\n');
fprintf('R^2 = %.6f\n', R2)
fprintf('Adj R^2 = %.6f\n', adjR2);
fprintf('RMSE = %.6f MPa\n', RMSE)
%
%% ---- Plot and Save ------------------------------------------------------
fig = figure; hold on;
plot(epsl, sig, 'b.','MarkerSize',3, 'LineWidth',3, 'DisplayName','Experimental');
plot(epsl, sigma_fit,'r-', 'LineWidth',3, 'DisplayName','Fitted model');
xlabel('Strain \epsilon');
ylabel('Stress \sigma_d (MPa)');
title(['Fitted Constitutive Model at \epsilon-dot = ', num2str(edot), ' s^{-1}']);
legend show; grid on;
% ---- Set axis scaling and ticks ----
axis([0 0.97 0 0.05]); % X: 0–1, Y: 0–0.05
xticks(0:0.1:1); % X increments of 0.02
yticks(0:0.01:0.05); % Y increments of 0.01
% Add a readable stats box
ax = gca;
xL = xlim(ax); yL = ylim(ax);
txt = {sprintf('R^2 = %.4f', R2), ...
sprintf('Adj R^2 = %.4f', adjR2), ...
sprintf('RMSE = %.3g MPa', RMSE)};
text(xL(1) + 0.02*range(xL), yL(2) - 0.04*range(yL), txt, ...
'VerticalAlignment','top', 'BackgroundColor','w', ...
'Margin',6, 'EdgeColor',[0.1 0.1 0.1], 'FontWeight','bold');
end
JAMEEL
on 4 Sep 2025
Matt J
on 4 Sep 2025
Can't find PUF28_5mmmin.xlsx. Please attach it to your last code presentation.
JAMEEL
on 4 Sep 2025
Both of your new .xlsx file attachments appear to contain identical data. Not sure if that was intentional...
Aside from that, the attached modification seems to work better on the new data:
file = 'data_Copy.xlsx';
%---- Known strain rates -------------------------------------------------
edot = 0.00308; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
% ---- Initial guesses & bounds ------------------------------------------
% [A, E, m, n, a, b, c]
p0 = [0.1, 500, 0.5, 2.0, 0.05, 50,0.1]; % initial params (MPa-based)
% [A, E, m, n, a, b]
lb = [-0.9, 50, -1, 1, 1, 0.01 ];
ub = [5, 500, 3, 2, 4.0, 0.1 ];
epsmax=0.85; %Maximum epsilon to fit
doFit(file, edot,edot0, p0,lb,ub,epsmax)
The old data seems to still work as well:
file = 'data.xlsx';
% ---- Known strain rates -------------------------------------------------
edot = 0.000617; % test strain rate [s^-1]
edot0 = edot; % reference strain rate [s^-1]
% ---- Initial guesses & bounds ------------------------------------------
% [A, E, m, n, a, b, c]
p0 = [0.9,500, 0.5, 2.0, 0.05, 50,0.07]; % initial params (MPa-based)
% [A, E, m, n, a, b]
lb = [-0.1, 1.2, -1, 1, 1, 0.01 ];
ub = [2, 20, 5, 5, 4.0, 0.5 ];
epsmax=0.95; %Maximum epsilon to fit
doFit(file, edot,edot0, p0,lb,ub,epsmax)
JAMEEL
on 4 Sep 2025
Hi Matt,
Thank you for your assistance. I’m a bit confused about why the same code you recently provided doesn’t fit the curve at different strain rates. Which settings should I adjust? When I use another strain rate, I don’t obtain the same curve you showed earlier
To demonstrate problems, please get into the habit of running the code here in the online environment, as I have been doing. A screen shot of the results alone doesn't tell us anything, because we cannot see what has changed in the way you are running the code.
By "another strain rate", do you mean a different .xlsx data set?
One thing that might be relevant to your last question. The known edot and edot0 parameters are not doing anything meaningful in your model, and should probably be removed. All they do is change the fitted values of a and b by a factor of log(edot/edot0). Attached is a version which removes edot,edot0 and also includes a routine to auto-guess the 'c' parameter (so, you no longer have to provide an initial guess or that parameter).
file = 'data_Copy.xlsx';
% [A, E, m, n, a, b]
p0 = [0.1, 500, 0.5, 2.0, 0, 0]; % initial params (MPa-based)
% [A, E, m, n, a, b]
lb = [-0.1, 1.2, -1, 1, -inf, -inf ];
ub = [2, 20, 5, 5, +inf, +inf ];
epsmax=0.85; %Maximum epsilon to fit
doFit(file, p0,lb,ub,epsmax)
%%
file = 'data.xlsx';
% ---- Initial guesses & bounds ------------------------------------------
% [A, E, m, n, a, b]
p0 = [0.9,500, 0.5, 2.0, 0, 0]; % initial params (MPa-based)
% [A, E, m, n, a, b]
lb = [-0.1, 1.2, -1, 1, -inf,-inf ];
ub = [2, 20, 5, 5, +inf, +inf ];
epsmax=0.85; %Maximum epsilon to fit
doFit(file, p0,lb,ub,epsmax)
Categories
Find more on Simulink Report Generator in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
