176 views (last 30 days)

Hi everyone,

I'm currently working on a fonction that need to find the inflexion point of a curve that I drew with only discrete datas (I have a abs vector with the abscissa and a vector y with ordinate). The typical shape of the curve look like this :

I don't use matlab a lot so maybe there's a very easy way to get it.

Here is an example for this plot :

x = [ 0 0.05 0.1 0.17 0.25 0.37 0.5 0.72 1 1.42 2 2.83 4 5.67 8 11.32 15 21.22 30 42.44 60 84.87 120 169.72 339.42 480 720 960 1200 1440];

y = [ 17.836 17.742 17.716 17.693 17.669 17.642 17.615 17.578 17.536 17.484 17.422 17.347 17.26 17.159 17.048 16.934 16.849 16.765 16.706 16.669 16.644 16.626 16.612 16.6 16.58 16.572 16.563 16.557 16.552 16.547];

The abcissa are first computed into their log10 form before being plotted?

I hope someone can help me but anyway, have a good day !

PS : I don't know if it's useful but I use matlab R2012a.

Azzi Abdelmalek
on 11 Aug 2016

x = [ 0 0.05 0.1 0.17 0.25 0.37 0.5 0.72 1 1.42 2 2.83 4 5.67 8 11.32 15 21.22 30 42.44 60 84.87 120 169.72 339.42 480 720 960 1200 1440];

y = [ 17.836 17.742 17.716 17.693 17.669 17.642 17.615 17.578 17.536 17.484 17.422 17.347 17.26 17.159 17.048 16.934 16.849 16.765 16.706 16.669 16.644 16.626 16.612 16.6 16.58 16.572 16.563 16.557 16.552 16.547];

xx=log10(x)

df1=gradient(y,xx)

df2=gradient(df1,xx)

id=sign(df2)

idx=strfind(id,[-1 1])

inflexionP=x(idx+1)

Azzi Abdelmalek
on 11 Aug 2016

Like I said, you can interpolate to find your 10000 points, you can use interp1 function

x = [ 0 0.05 0.1 0.17 0.25 0.37 0.5 0.72 1 1.42 2 2.83 4 5.67 8 11.32 15 21.22 30 42.44 60 84.87 120 169.72 339.42 480 720 960 1200 1440];

y = [ 17.836 17.742 17.716 17.693 17.669 17.642 17.615 17.578 17.536 17.484 17.422 17.347 17.26 17.159 17.048 16.934 16.849 16.765 16.706 16.669 16.644 16.626 16.612 16.6 16.58 16.572 16.563 16.557 16.552 16.547];

xi=linspace(min(x),max(x),10000);

yi=interp1(x,y,xi)

John D'Errico
on 11 Aug 2016

Edited: John D'Errico
on 11 Aug 2016

You have a nice smooth curve. So you can use an interpolating spline.

xlog = log10(x);

pp = spline(xlog(2:end),y(2:end));

fnplt(pp)

hold on

plot(xlog,y,'o')

Note that I fit the curve from the second point on, because the log function is singular at 0, approaching -inf as a limit.

The problem is, finding an inflection point can be difficult to pin down exactly, because it is equivalent to estimating the second derivative of a curve, only given data. This is one variety of something called an ill-posed problem. Differentiation of data tends to amplify any noise in the data, and do so greatly. Thus, a good reason for why is called ill-posed.

So, we can now differentiate that cubic spline twice, then look for where the second derivative crosses zero.

ppd2 = fnder(pp,2);

fnplt(ppd2)

Ok, as I said, it gets a bit nasty looking at the ends, because as I said, differentiation amplifies noise. Two derivatives make it worse. We can find all locations where the curve crosses zero easily enough though.

Here, I'll just employ a tool from my SLM toolbox.

d2zerolocs = slmsolve(ppd2,0)

d2zerolocs =

-1.0711 -0.63129 -0.58148 0.93178 2.6653 2.7129 2.929

slmsolve uses interpolation on the actual function to find all locations on the curve that satisfy the given y value. In this case, I'm using it to find all locations that satisfy the given second derivative value.

We know from the plot that the location of interest is near 1, on the log scale, so log(x) is approximately 0.93. (In fact, I'd not even trust that location to more than about 1 significant digit.) If you want that in the original x domain, just exponentiate using 10^().

10.^d2zerolocs(4)

ans =

8.5462

I suppose, if you are going to use my slmtoolbox anyway, to do the solve, then you could afford to use it to do the fit too.

slm = slmengine(xlog(2:end),y(2:end),'plot','on','decreasing','on','knots',8,'result','pp');

So maybe a little smoother.

ppd2 = fnder(pp,2);

fnplt(ppd2)

grid on

ppd2 = fnder(slm,2);

fnplt(ppd2)

d2zerolocs = slmsolve(ppd2,0)

d2zerolocs =

-1.0979 0.87456 2.4792 2.6467

This time, it came out a wee bit below 0.9. As I said before, I'd not trust ANY estimator of the location of an inflection point to more than about ONE significant digit. So be careful. Computers return 16 digits. But here, almost all of those digits are meaningless, virtually random garbage.

10.^d2zerolocs(2)

ans =

7.4914

I just noticed that the posted version did not have slmsolve in it yet. So I've posted a new release, that does include slmsolve. You may need to wait a few minutes before it becomes visible to the world.

RAYDEL RAMOS
on 16 Jun 2019

The "slmsolve" function is available in matlab 2019a?

I would like to get all the inflexion points as explained in several comments but "slmsolve" is not defined in matlab and then I have not been able to get the desired results. Whenever I run the script I get the following alert message

"Indefinite or variable function 'slmsolve'"

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.