Getting the linear portion of a non-linear curve

Hello,

I am trying to extract the linear section of a non-linear curve. (See picture attached - I would like to get only the section between the red lines). I also attach the data for that curve - it is data obtained from mechanical testing.

I have tried taking the 2nd derivative but so far I have not managed to do this.

Would you have any suggestions?

Thank you.

Salomé

[Moein M]

Hello Salomé

Could you please post the full code so that others can use it too?

That would be awesome.

Answer 1

[Steven Lord]

The ischange function may be of use to you. Depending on what you're planning to do with this data the detrend function may also be of interest.

[Salomé Sanchez]

Hello,

Thanks for your suggestion.

I used ischange to find the bounds where the curve is linear and then I only plotted the values between the bounds.

TF=ischange(y,'linear'); % Using ischange to find the "abrupt changes" in the graph
brkpt=x(TF==1); % Getting the breakpoint values of x where the changes happen
 
linear_x= x(x>=brkpt(1) & x<=brkpt(2)); % In my case, the linear section seemed to be between these 2 breaks
linear_y= fct(linear_x); % Putting the x values of the linear section into my function

I attach a picture of the results.

Thank you!

Salomé

Answer 2

[John D'Errico]

As is often the case on a question like this, you have been too vague for us to comprehensively understand what you do "have". (That is why there are so many conjectures on how to solve your problem. But none of us really know what you have or need to learn.)

That is, is the picture you show a picture of some points taken from some known function? Probably not, since then you would know what the "linear" portion would be. So then you must have some relationhip, sampled at a list of points. You probably have (x,y) pairs only. Are they sampled at some uniform spacing in x? Or is this a question where you have some general nonlinear function, where you do know the function (i.e., it can be written down on paper) and you want to find some interval in x where the function deviates from linearity by the least amount? Of course, then you should recognize that no nonlinear function is truly linear, except over some short interval, where essentially any continuously differentiable function is locally linear everywhere.

Piecewise is irrelevant in any case, since piecewise is not a tool that can find anything.

But if you want to find some interval where a known differentiable nonlinear function is most linear, then you might just compute the second derivative of said function. Then search for the longest interval where that second derivative is less than some tolerance over the entire interval.

If what you have is a set of sampled points only from some relationship, do they have noise in them? So is the picture you show just some generic example you have cooked up? If you do have sampled data, then you can use the gradient functino to compute a second deritive estimate along the curve. Now again, just search for the longest interval where the absolute value of the second derivative is less than some tolerance.

Another possibility is to interpolate your data with a spline interpolant. Now take the second derivative of the cubic spline, which will now be a piecewise linear function. Again, find the longest interval where said second derivative function is less than some tolerance in absolute value.

If there is noise in your sampled data, then you need to smooth it first. In that case, a smoothing spline probably makes sense. Once you have the smoothing spline, again, differentiate twice.

In the end, I can see a zillion variations on what you might really have, and what you really need in the end.

[Salomé Sanchez]

Hello,

Thanks for your answer.

I have now attached my data in my original post. There is no noise in my data - I already smoothed it.

Answer 3

[Matt J]

I'm not sure why a 2nd derivative test wouldn't have worked, as long as your points are noiseless:

i=find( abs(diff(x,2))>something_small   ,1);

xlinear=x(1:i);

[Image Analyst]

To build on Matt's idea with additional code:

% Create sample data - an exponential.
x = linspace(1, 10, 1000);
y = 0.20 * exp(x);
% Make the part from 1 to 7 linear
y(1:700) = linspace(y(1), y(700), 700);
% Now we have sample data.
subplot(3, 1, 1);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
xlabel('x', 'FontSize', 20);
ylabel('y', 'FontSize', 20);
title('y vs. x', 'FontSize', 20);
xl = xlim(); % Remember range of x axis for later plotting of linear portion.
% Matt's code:
deriv2 = abs(diff(x,2));
subplot(3, 1, 2);
plot(x(1:end-2), deriv2, 'b-', 'LineWidth', 2);
xlabel('x', 'FontSize', 20);
ylabel('Second Derivative', 'FontSize', 20);
grid on;
title('Second Derivative', 'FontSize', 20);
something_small = 1.0e-15;
yline(something_small);
index = find(deriv2 > something_small , 1);
xlinear = x(1:index);
% Plot the linear section.
subplot(3, 1, 3);
plot(xlinear, y(1:index), 'r-', 'LineWidth', 2);
xlabel('x', 'FontSize', 20);
ylabel('y', 'FontSize', 20);
grid on;
title('Linear Portion', 'FontSize', 20);
xlim(xl); % Make sure x limits are the same as the original data.

Attach your actual data if you want more help.

[Salomé Sanchez]

Hello!

Thank you for your answer!

I think the problem with the 2nd derivative is that there are so many datapoints that from one to the other, the 2nd derivative doesn't change much - maybe my bounds were not adequate though.

(I have now attached my data in my original post.)

I will give a go to what you suggested, thanks!

Answer 4

[Adam Danz]

If you have the x values of the red lines, it's as easy as

bounds = [x1,x2];   % [lower, upper] bounds (red lines)
keepIdx = x >= bounds(1) && x <= bounds(2); 
x(~keepIdx) = []; 
y(~keepIdx) = []; 

If you're trying to find the upper bound of the linear portion, where the curve starts to become nonlinear, you could fit a line to the first part of the curve, measure the error between the fit and the curve, and detect when the error increases beyond some threshold close to 0. Note that Steven Lord's suggestion to use ischange or detrend is a better first-attempt than this.

To fit a line to the first part of the curve, isolate the a section of the linear part of the line and use p = polyfit(x,y,n) to get the slope and intercept.

bounds = [x1,x2];   % [lower, upper] bounds of a linear section of the line
keepIdx = x >= bounds(1) && x <= bounds(2); 
p = polyfit(x(keepIdx), y(keepIdx),1)

Then compute the error between the fit and your curve

yhat = p(1)*x + p(2); 
err = abs(y-yhat); 
threshold = 0.05; % ???
nonlinear_starting_index = find(err,1,'first');

[Salomé Sanchez]

Hello,

Thank you for your answer!

I do not know the position of the red lines - it's what I need to find.

I will try the ischange first then try your suggestion.

Thank you, I will let you know if it works!

[Salomé Sanchez]

Thanks for your help!

I used ischange to find the bounds and then used your method to get the values between the bounds.

Thanks a lot!

[Adam Danz]

Glad I could help!

Answer 5

[Salomé Sanchez]

Hello,

Using Steven Lord and Adam Danz's suggestions, I used ischange to finds the bounds where the curve is linear and then I only plotted the values between the bounds.

TF=ischange(y,'linear'); % Using ischange to find the "abrupt changes" in the graph
brkpt=x(TF==1); % Getting the breakpoint values of x where the changes happen
 
linear_x= x(x>=brkpt(1) & x<=brkpt(2)); % In my case, the linear section seemed to be between these 2 breaks
linear_y= fct(linear_x); % Putting the x values of the linear section into my function

I attach a picture of the results.

Thank you!

Salomé

[Md Mezbah Uddin]

Could you please tell me what is 'fct' in the last line?