How do I determine the linear part of a force-displacement curve?

I have a force-displacement graph (the attached file) and I want to find out exactly at what point the linear portion of the graph ends. Please see the picture for more clarification.

I've already smoothed out the curve, so it is noiseless.

I did a quick search but couldn't find anything useful.

Any help would be highly appreciated.

Answer 1

[Star Strider]

Try this:

T1 = readtable('Moein M force-displacement.xlsx');
T1.Properties.VariableNames = {'Force','Displacement'};

ThrVal = 2.5E+11;
[TF,Sl,Ic] = ischange(T1.Displacement, 'linear','Threshold',ThrVal);
Idx = find(TF);

figure
plot(T1.Force, T1.Displacement)
hold on
plot(T1.Force(Idx(1)), T1.Displacement(Idx(1)), '+r')
hold off
grid
text(T1.Force(Idx(1)), T1.Displacement(Idx(1)), sprintf('\\leftarrowLinear Region End:\n    Force              = %10.8f\n    Displacement = %10.2f', T1.Force(Idx(1)), T1.Displacement(Idx(1))), 'HorizontalAlignment','left', 'VerticalAlignment','top')

producing:

Experiment with ‘ThrVal’ to get the result you want. The value I chose appears to do what you want, however it may need some minor adjustments.

The ischange function was introduced in R2017b. Another option is the Signal Processing Toolbox findchangepts function (with different argumnents and behaviour), introduced in R2016a.

[Moein M]

Wonderful. Thank you!

[Star Strider]

As always, my pleasure!

[schilakwad]

Hi, how did you select this particular threshold value? I'm trying to do something very similar and am unsure of the appropriate ThrVal to go for.

[Star Strider]

schilakwad —

‘... how did you select this particular threshold value?’

I doubt that there is any specific way of actually approximating it through calculation. Note that I got all 3 outputs from ischange, so one of the criteria I use is the second one (slope, with the 'linear' option), plot it separately, and see where it changes, usually with respect to a tolerance. (I do not remember exactly how I arrived at the specific value of ‘ThrVal’ here.)

[schilakwad]

Thank you!

[Marzieh Nodeh]

Hi

I`m tring to use these codes for my Load-displacment curve to determine the linear part of the curve but I receive the following error.

Could you please help me.

I attache the excel file. Thank you in advance.

plot(T1.Force(Idx(1)), T1.Displacement(Idx(1)), '+r')
hold off
grid
Index exceeds the number of array elements (0).

[Zoe Man]

Hi @Star Strider

I'm also getting this issue where it says this so I checked what idx was and:

idx = 0x1 empty double column vector

index exceeds the number of array elements. Index must not exceed 0

Would you be able to help please?

Thanks

[Star Strider]

@Marzieh Nodeh (and probably @Zoe Man as well) — Change ‘ThrVal’ to a value that makes sense in the context of the data, then choose the appropriate element of ‘idx’. (I doubt that there is a way to make my code automatically adapt to the data.)

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/993930/Marzieh.xlsx');
T1.Properties.VariableNames = {'Force','Displacement'}

T1 = 8361×2 table
      Force      Displacement
    _________    ____________

            0              0 
    0.0015162     0.00070096 
    0.0030368     0.00054254 
    0.0043883     0.00036289 
    0.0046311     0.00092839 
    0.0076054      0.0010386 
    0.0094601     0.00066402 
     0.011589     0.00086091 
     0.013755      0.0015124 
     0.015513      0.0011822 
     0.017806      0.0012508 
     0.018115      0.0015146 
     0.018571       0.001629 
     0.020005      0.0012286 
     0.024438      0.0018267 
      0.02327      0.0019344 

ThrVal = 1.0E-1;                                                        % <— Changed
[TF,Sl,Ic] = ischange(T1.Displacement, 'linear','Threshold',ThrVal);
Idx = find(TF)

Idx = 3×1

852 3440 7971

figure
plot(T1.Force, T1.Displacement)
hold on
plot(T1.Force(Idx(1)), T1.Displacement(Idx(1)), '+r')
hold off
grid
text(T1.Force(Idx(1)), T1.Displacement(Idx(1)), sprintf('\\leftarrowLinear Region End:\n    Force              = %10.4f\n    Displacement = %10.4f', T1.Force(Idx(1)), T1.Displacement(Idx(1))), 'HorizontalAlignment','left', 'VerticalAlignment','top')

.

[Florian]

Hi all

I would be interested in this (and the most state of the art solution with R2025B) as well! (:

Currntly investigating my force / displacement data and want to derive insights on the behaviour to validate my FEM.

What is the best way to calculate the turning points of the curves in a force / displacement data set? Still manual selection or any tool-based solution? Tried to use the derivative and ended up with zig-zags, which become worse for any additional derivative.

Is there any way to find a simplified interpolation first and then calculate the derivative?

Thanks

FW

[Florian]

this is my try with copilot and some assumption (smoothening between each derivation, constraint for starting with 0, ending with ~0 for d2 and d3 and so on

any more reliable idea?

script:

%% ROW 1: Top — EXACT Fig.5 curves (no x labels/ticks)

ax1 = subplot(4,1,1); hold(ax1,'on'); grid(ax1,'on');

title(ax1, 'Mean functions (Fig.5 reused)');

ylabel(ax1,'Force [kN]');

xlim(ax1, xlimAll); ylim(ax1, ylimForce);

xlabel(ax1,''); set(ax1,'XTickLabel',[]);

for idx = 1:4

if ~exist('Fig5Curves','var') || numel(Fig5Curves) < idx, continue; end

% ----- (O) -----

if isfield(Fig5Curves(idx),'O') && isfield(Fig5Curves(idx).O,'x')

xO = Fig5Curves(idx).O.x(:).'; yO = Fig5Curves(idx).O.y(:).';

% Visual marker in top plot (centered position for readability)

[~,iM_topO] = min(abs(xO - mean([min(xO) max(xO)])));

hO = plot(ax1, xO, yO, 'Color',[1 0.4 0],'LineWidth',2.4, ...

'Marker',markers{idx}, 'MarkerIndices',iM_topO, 'MarkerFaceColor',[1 0.4 0], ...

'DisplayName', sprintf('O: %s', Fig5Curves(idx).name));

legendHandles(end+1) = hO; legendLabels{end+1} = get(hO,'DisplayName');

% Derivative chain (no left padding). Operate on uniform Fig.5 grid.

N = numel(xO);

w1 = oddwin(N, SG_W_D1); w2 = oddwin(N, SG_W_D2); w3 = oddwin(N, SG_W_D3);

% dF/dx: numerical gradient -> SG smoothing -> right-end taper to 0

d1_raw = gradient(yO, xO);

d1_sm = sgolayfilt(d1_raw, SG_P, w1);

d1_sm = f6_taper_right_to_zero(xO, d1_sm, TAPER_RIGHT_PTS);

d1_sm = sgolayfilt(d1_sm, SG_P, oddwin(numel(d1_sm), SG_W_D1));

% d2: from smoothed d1 -> SG smoothing

% leftmost points set exactly to 0; right-end taper to 0

d2_raw = gradient(d1_sm, xO);

d2_sm = sgolayfilt(d2_raw, SG_P, w2);

d2_sm = f6_set_left_zero(d2_sm, LEFT_ZERO_PTS);

d2_sm = f6_taper_right_to_zero(xO, d2_sm, TAPER_RIGHT_PTS);

d2_sm = sgolayfilt(d2_sm, SG_P, oddwin(numel(d2_sm), SG_W_D2));

% d3: from smoothed d2 -> SG smoothing

% leftmost points set exactly to 0; right-end taper to 0

d3_raw = gradient(d2_sm, xO);

d3_sm = sgolayfilt(d3_raw, SG_P, w3);

d3_sm = f6_set_left_zero(d3_sm, LEFT_ZERO_PTS);

d3_sm = f6_taper_right_to_zero(xO, d3_sm, TAPER_RIGHT_PTS);

d3_sm = sgolayfilt(d3_sm, SG_P, oddwin(numel(d3_sm), SG_W_D3));

Xvals(idx).O = xO;

Deriv1(idx).O = d1_sm;

Deriv2(idx).O = d2_sm;

Deriv3(idx).O = d3_sm;

end

% ----- (B) -----

if isfield(Fig5Curves(idx),'B') && isfield(Fig5Curves(idx).B,'x')

xB = Fig5Curves(idx).B.x(:).'; yB = Fig5Curves(idx).B.y(:).';

[~,iM_topB] = min(abs(xB - mean([min(xB) max(xB)])));

hB = plot(ax1, xB, yB, 'Color',[0 0.2 1],'LineWidth',2.4, ...

'Marker',markers{idx}, 'MarkerIndices',iM_topB, 'MarkerFaceColor',[0 0.2 1], ...

'DisplayName', sprintf('B: %s', Fig5Curves(idx).name));

legendHandles(end+1) = hB; legendLabels{end+1} = get(hB,'DisplayName');

N = numel(xB);

w1 = oddwin(N, SG_W_D1); w2 = oddwin(N, SG_W_D2); w3 = oddwin(N, SG_W_D3);

d1_raw = gradient(yB, xB);

d1_sm = sgolayfilt(d1_raw, SG_P, w1);

d1_sm = f6_taper_right_to_zero(xB, d1_sm, TAPER_RIGHT_PTS);

d1_sm = sgolayfilt(d1_sm, SG_P, oddwin(numel(d1_sm), SG_W_D1));

d2_raw = gradient(d1_sm, xB);

d2_sm = sgolayfilt(d2_raw, SG_P, w2);

d2_sm = f6_set_left_zero(d2_sm, LEFT_ZERO_PTS);

d2_sm = f6_taper_right_to_zero(xB, d2_sm, TAPER_RIGHT_PTS);

d2_sm = sgolayfilt(d2_sm, SG_P, oddwin(numel(d2_sm), SG_W_D2));

d3_raw = gradient(d2_sm, xB);

d3_sm = sgolayfilt(d3_raw, SG_P, w3);

d3_sm = f6_set_left_zero(d3_sm, LEFT_ZERO_PTS);

d3_sm = f6_taper_right_to_zero(xB, d3_sm, TAPER_RIGHT_PTS);

d3_sm = sgolayfilt(d3_sm, SG_P, oddwin(numel(d3_sm), SG_W_D3));

Xvals(idx).B = xB;

Deriv1(idx).B = d1_sm;

Deriv2(idx).B = d2_sm;

Deriv3(idx).B = d3_sm;

end

end

%% ROW 2: dF/dx — markers at global extremum (max |.|); manual y-limits

ax2 = subplot(4,1,2); hold(ax2,'on'); grid(ax2,'on');

title(ax2,'1st derivative dF/dx (smoothed; right end → 0 near x≈0)');

ylabel(ax2,'dF/dx [kN/mm]'); xlim(ax2, xlimAll);

xlabel(ax2,''); set(ax2,'XTickLabel',[]);

ylim(ax2, USER_YLIM_D1);

for idx = 1:4

if isfield(Xvals(idx),'O')

[~,iM] = max(abs(Deriv1(idx).O));

plot(ax2, Xvals(idx).O, Deriv1(idx).O, 'Color',[1 0.4 0], 'LineWidth',2, ...

'Marker',markers{idx}, 'MarkerIndices',iM, 'MarkerFaceColor',[1 0.4 0]);

end

if isfield(Xvals(idx),'B')

[~,iM] = max(abs(Deriv1(idx).B));

plot(ax2, Xvals(idx).B, Deriv1(idx).B, 'Color',[0 0.2 1], 'LineWidth',2, ...

'Marker',markers{idx}, 'MarkerIndices',iM, 'MarkerFaceColor',[0 0.2 1]);

end

end

%% ROW 3: d^2F/dx^2 — markers at global extremum; left=0, right→0

ax3 = subplot(4,1,3); hold(ax3,'on'); grid(ax3,'on');

title(ax3,'2nd derivative d^2F/dx^2 (from smoothed dF/dx; left=0, right→0)');

ylabel(ax3,'d^2F/dx^2 [kN/mm^2]'); xlim(ax3, xlimAll);

xlabel(ax3,''); set(ax3,'XTickLabel',[]);

ylim(ax3, USER_YLIM_D2);

for idx = 1:4

if isfield(Xvals(idx),'O')

[~,iM] = max(abs(Deriv2(idx).O));

plot(ax3, Xvals(idx).O, Deriv2(idx).O, 'Color',[1 0.4 0], 'LineWidth',2, ...

'Marker',markers{idx}, 'MarkerIndices',iM, 'MarkerFaceColor',[1 0.4 0]);

end

if isfield(Xvals(idx),'B')

[~,iM] = max(abs(Deriv2(idx).B));

plot(ax3, Xvals(idx).B, Deriv2(idx).B, 'Color',[0 0.2 1], 'LineWidth',2, ...

'Marker',markers{idx}, 'MarkerIndices',iM, 'MarkerFaceColor',[0 0.2 1]);

end

end

%% ROW 4: d^3F/dx^3 — markers at global extremum; left=0, right→0 (x label shown)

ax4 = subplot(4,1,4); hold(ax4,'on'); grid(ax4,'on');

title(ax4,'3rd derivative d^3F/dx^3 (from smoothed curvature; left=0, right→0)');

ylabel(ax4,'d^3F/dx^3 [kN/mm^3]'); xlabel(ax4,'Displacement [mm]');

xlim(ax4, xlimAll);

ylim(ax4, USER_YLIM_D3);

for idx = 1:4

if isfield(Xvals(idx),'O')

[~,iM] = max(abs(Deriv3(idx).O));

plot(ax4, Xvals(idx).O, Deriv3(idx).O, 'Color',[1 0.4 0], 'LineWidth',2, ...

'Marker',markers{idx}, 'MarkerIndices',iM, 'MarkerFaceColor',[1 0.4 0]);

end

if isfield(Xvals(idx),'B')

[~,iM] = max(abs(Deriv3(idx).B));

plot(ax4, Xvals(idx).B, Deriv3(idx).B, 'Color',[0 0.2 1], 'LineWidth',2, ...

'Marker',markers{idx}, 'MarkerIndices',iM, 'MarkerFaceColor',[0 0.2 1]);

end

end

% Shared legend across the four subplots (placed at bottom)

lgd = legend(legendHandles, legendLabels, 'Orientation','horizontal', 'NumColumns', 4, 'Box','off');

set(lgd,'Units','normalized'); set(lgd,'Position',[0.1 0.01 0.8 0.04]);

saveas(figure6, fullfile(outputFolder, 'Figure6_chainSmooth_globalExtrema_manualYLim.svg'));

disp('All figures (1–6) generated successfully.');

%% ========================================================================

% LOCAL FUNCTIONS

% ========================================================================

function [disp, force] = g_readDatFile(fp)

% Parse MTS ASCII files; capture displacement and force (kN) columns.

fid = fopen(fp,'r');

if fid<0, error('Cannot open file: %s', fp); end

disp = []; force = []; start=false;

while true

L = fgetl(fid);

if ~ischar(L), break; end

if contains(L,'Axial Weg')

start=true; fgetl(fid); continue;

end

if ~start, continue; end

L = strrep(L, ',', '.'); % decimal comma -> dot

vals = sscanf(L,'%f %f %f %f'); % parse four columns

if numel(vals)==4

disp(end+1,1)=vals(1); % displacement

force(end+1,1)=vals(3); % force (kN)

end

end

fclose(fid);

end

function [typ,id] = g_parseCellIdentifier(fn)

% Extract pack identifier: e.g., 'O35' or 'B12' -> type 'O'/'B', numeric id.

t = regexp(fn,'(O|B)([0-9]+)','tokens','once');

typ = t{1};

id = str2double(t{2});

end

function g_drawConnector(xcurve, ycurve, labelX, labelY)

% Draw a thin connector from curve to label staging area near x = -2 mm.

targetX = -2;

[~, idx] = min(abs(xcurve - targetX));

plot([xcurve(idx) labelX], [ycurve(idx) labelY], ...

'Color',[0.6 0.6 0.6], 'LineWidth',0.8);

end

function s = g_sanitize(nameStr)

% Produce filesystem-safe filename tokens from titles.

s = regexprep(nameStr,'[^a-zA-Z0-9()]','_');

end

function y_out = f6_taper_right_to_zero(x, y, N_right)

% Gently fade the last N_right samples (toward x ≈ 0) to exactly zero.

% Assumes x is increasing and ends at/near 0 (as in Fig.5 grids).

y_out = y(:).';

x = x(:).';

n = numel(x);

N_right = max(0, min(n, round(N_right)));

if N_right < 1

y_out = y_out(:).'; return;

end

[~, iR] = min(abs(x - 0)); % index closest to x = 0

iR = max(1, min(n, iR));

iR1 = max(1, iR - max(1,N_right) + 1);

seg = iR1:iR;

t = linspace(0,1,numel(seg)); % fade factor from 1 -> 0

y_out(seg) = y_out(seg) .* (1 - t);

y_out(iR:end) = 0;

y_out = y_out(:).';

end

function y_out = f6_set_left_zero(y, N_left)

% Hard-set the first N_left samples exactly to zero (no taper).

y_out = y(:).';

n = numel(y_out);

N_left = max(0, min(n, round(N_left)));

if N_left > 0

y_out(1:N_left) = 0;

end

y_out = y_out(:).';

end

[Star Strider]

@Florian -- I'm not certain that I understand what you are doing. I do not have your data, and I cannot get it from the images.

If you wnat to calculate the derivatives, one option would be to use the gradient function. The problem is that the numerical derivative is by definition a highpass filter, and will amplify the noise. You could use a lowpass filter first, however that might not reduce the noise enough to provide a useful result.

Another option would be to use the polyfit function to get the closest fit possible to your data, then polyder to take the respective derivatives. You can use those coefficients with polyval to evaluate them.

The sgolayfilt function creates a multiband FIR 'notch' filter that filters out irregularitties in the signal. That works, however any filtering approach needs to be considered carefully. I would plot the Fourier transforms of the unfiltered and filtered signals to be certain I got the result I wanted.

----------

NOTE: The original question wanted the 'linear' portion of the curve, and by that I understood it to be the first part, as depicted in the original plot. Moein M wanted to define the breakpoint between that part and the rest of the curve. I understand that the curve has no truly llinear part. I used their definition to get the desired result.

[Florian]

thanks for the recommendations, @Star Strider!

yes, questions are different, but the identified point when it becomes non-linear (if we consider the start linear) would help to split my functions as well.

One example of my data (as mean of the multiple measurements) attached.

[Star Strider]

My pleasure!

I am not certain where you are going with the derivatives, however this would be my approach --

T1 = readtable('Fig5_Mean_centre_B.csv', VariableNamingRule='preserve');
x_mm = T1.(1);
F_kN = T1.(2);

B = polyfit(x_mm, F_kN, 14);
F_fit = polyval(B, x_mm);

figure
hp1 = plot(x_mm, F_kN, DisplayName='Original Data');
hold on
hp2 = plot(x_mm, F_fit, DisplayName=["Polynomial Fit (Order="+(numel(B)-1)+")"]);
hold off
grid
% legend([hp1 hp2], Location='best', Interpreter='LaTeX')
% legend("Original Data", ["Polynomial Fit (Order="+(numel(B)-1)+")"], "Roots of $\frac{d^2F}{dx^2}$", Location="best", Interpreter="LaTeX")
xlim([min(x_mm) max(x_mm)])

ax1 = gca;

dFdx = polyder(B);
d2Fdx2 = polyder(dFdx);

dF_eval = polyval(dFdx, x_mm);
d2F_eval = polyval(d2Fdx2, x_mm);

dF_rts = roots(dFdx)

dF_rts = 
  -6.9528 + 0.0000i
  -2.9805 + 0.1341i
  -2.9805 - 0.1341i
  -2.6188 + 0.3932i
  -2.6188 - 0.3932i
  -1.9937 + 0.4578i
  -1.9937 - 0.4578i
  -1.3149 + 0.5304i
  -1.3149 - 0.5304i
  -0.5090 + 0.2387i
  -0.5090 - 0.2387i
  -0.1527 + 0.0000i
  -0.0475 + 0.0000i

d2F_rts = roots(d2Fdx2)

d2F_rts = 
  -6.5550 + 0.0000i
  -2.8997 + 0.0000i
  -2.6647 + 0.0000i
  -2.5742 + 0.0000i
  -2.3402 + 0.0000i
  -2.0744 + 0.0000i
  -1.4781 + 0.2865i
  -1.4781 - 0.2865i
  -0.9803 + 0.0000i
  -0.5023 + 0.0000i
  -0.3507 + 0.0000i
  -0.0900 + 0.0000i

figure
tiledlayout(2,1)
nexttile
plot(x_mm, dF_eval, DisplayName="$\frac{dF}{dx}$")
grid
xline(real(d2F_rts))
xlim([min(x_mm) max(x_mm)])
title("$\frac{dF}{dx}$", Interpreter="LaTeX")
nexttile
plot(x_mm, d2F_eval, DisplayName="$\frac{d^2F}{dx^2}$")
grid
xline(real(d2F_rts))
xl1 = xline(ax1, real(d2F_rts), DisplayName="Roots of $\frac{d^2F}{dx^2}$");

xlim([min(x_mm) max(x_mm)])
ylim([-5 5])
title("$\frac{d^2F}{dx^2}$", Interpreter="LaTeX")

% legend(Location='best', Interpreter="LaTeX")

The xline lines are the roots of and (of course) correspond to the maxima and minima of .

Limit their number by thresholding them, for example --

use_roots = d2F_rts(d2F_rts >= -1.5)

use_roots = 
  -1.4781 + 0.2865i
  -1.4781 - 0.2865i
  -0.9803 + 0.0000i
  -0.5023 + 0.0000i
  -0.3507 + 0.0000i
  -0.0900 + 0.0000i

.

Answer 2

[John D'Errico]

You have already accepted Star's answer, so you may not even see this answer. But in fact, there is no linear part of such a curve. A perfectly elastic material should show an initially linear behavior, true. So is this force-displacement curve representative of a perfectly linearly elastic material?

plot(F,D,'-')

But now let us look more carefuly at that initial portion.

axis([0 .1 0 6e5])

So merely by viewing that portion where you think the curve is "linear", we see it does not look at all linear.

K = F<0.1;
mdl = fittype('a*F^2 + b*F','indep','F');
% Note that the model is linear in the parameters, but fit is too dumb to see that.
% so any starting values will eliminate the irritating warning message from fit
% about the lack of starting values.
mdlest = fit(F(K),D(K),mdl,'start',[1 1])
mdlest = 
     General model:
     mdlest(F) = a*F^2 + b*F
     Coefficients (with 95% confidence bounds):
       a =  -3.617e+07  (-3.623e+07, -3.611e+07)
       b =   8.393e+06  (8.389e+06, 8.398e+06)

I chose a model with no constant term, since we know that when the force is zero, so must be the displacement.

What is important to see in that fit is the fit is remarkably good over that region, and that the quadratic coefficient is actually more non-zero than the linear coefficient.

plot(mdlest)
hold on
plot(F(K),D(K),'b.')
xlabel 'Force'
ylabel 'Displacement'

And that leaves me with the claim this curve is not that of an perfectly elastic material (within bounds before it undergoes plastic deformation) but that of a nonlinearly elastic material. There is essentially NO linear section of the start of that curve.

At best, you may decide to approximate it with a linear section, but that would mistake the true behavior seen here. As well, if this really is the force-displacement curve (I wonder how much you smoothed the curve) I would look to see if the material undergoes plastic deformation, so after you stretch it by a small amount, and then release it, does that mateiral return to the original rest length?

Finally, there is a distinct possibility that whatever tool you used to "smooth" this curve was used inappropriately, causing the curve to APPEAR to be so strongly quadratic. I cannot know for certain if this is true, but it seems highly likely, as that section of the curve is clearly very strongly a quadratic polynomial. If that is the case, then your smoothing of the data made it impossible to truly understand the behavior of this material.

Answer 3

[Sam Chak]

Hi Moein M,

From the Excel file, you can check and verify that the first 90 points (up to 0.089) are covered by the circled region. Please note that the region up to 0.027 can be approximated by a linear function with a high (R-squared) value.

[Moein M]

Thank you for your feedback.

The ellipse was just to give a rough estimation of where the linear portion of the graph lies. It's not exact.

Anyway, thanks again.