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')

.

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.