Output doesn't display a value, just an empty space.

So, I've got 2 pulses (imp and S0) and their sum (Si). I've got to find minimum of sum (Si). And I approximate the sum with a polynome. I took derivatives (1-st and 2-nd order) of the sum and polynome.

I analyze derivatives, and decided to look for minimum of sum (Si), polynome and 2-nd derivatives.

But, in some cases the result is just empty, but size of variable is 1x999 double.

lb = -5;
ub = 15;                            
x_imp = linspace(lb, ub, 1001);    
x = linspace(lb, ub, 1001);

imp = sech(x_imp);                  
S0 = sech(x-10); 

imp1 = sech(x - 1);
S1 = 0.5 * (imp1 + S0);

xmin = -1;
xmax = 10;
% find inflection point Si
range_x = x(x>=xmin & x<=xmax);
range_S1 = S1(x>=xmin & x<=xmax);   
MinPtsS1 = islocalmin(range_S1);            
fprintf ('x_S1 = %f , y_S1 = %f .\n', range_x(MinPtsS1), range_S1(MinPtsS1));

% polynom approximation pv
[p_S1, ~, mu] = polyfit(x, S1, 26); 
pv_S1 = polyval (p_S1, x, [], mu); 

% find inflection point of polynom pv

range_x = x(x>=xmin & x<=xmax);
range_pvS1 = pv_S1(x>=xmin & x<=xmax);   
MinPts_pvS1 = islocalmin(range_pvS1); 
fprintf('x_pvS1 = %f , y_pvS1 = %f .\n', range_x(MinPts_pvS1), range_pvS1(MinPts_pvS1));

    % DIFF find inflection of diff2 Si, corresponding to inflection point of Si
        dS1 = 20*diff (S1);
        d2S1 = 10*diff(dS1);
        range_x = x(x>=xmin & x<=xmin);
        range_d2S1 = d2S1(x>=xmin & x<=xmin);   
        MinPts_d2S1 = islocalmin(range_d2S1); 
        fprintf('x_d2S1 = %f , y_d2S1 = %f .\n', range_x(MinPts_d2S1), range_d2S1(MinPts_d2S1));
    % DIFF find inflection of diff2 polynom Si, corresponding to inflection point of Si
        dpS1 = 20*diff(pv_S1);
        dp2S1 = 20*diff(dpS1);
        range_x = x(x>=xmin & x<=xmax);
        range_dp2S1 = dp2S1(x>=xmin & x<=xmax);   
        MinPts_dp2S1 = islocalmin(range_dp2S1); 
        fprintf('x_dp2S1 = %f , y_dp2S1 = %f .\n', range_x(MinPts_dp2S1), range_dp2S1(MinPts_dp2S1));

Answer 1

[Stephen23]

How many values fulfill your logical comparisons? Exactly one.

range_x = x(x>=xmin & x<=xmin);
range_d2S1 = d2S1(x>=xmin & x<=xmin);

should be

range_x = x(x>=xmin & x<=xmax);
range_d2S1 = d2S1(x>=xmin & x<=xmax);

Due to FPRINTF processing all input arguments in linear order you will also need to modify the FPRINTF input arguments (e.g. by defining a matrix), or use e.g. COMPOSE with column vectors.

There might be other bugs too...

[Stephen23]

I modified the FPRINTF calls:

lb = -5;
ub = 15;                            
x_imp = linspace(lb, ub, 1001);    
x = linspace(lb, ub, 1001);

imp = sech(x_imp);                  
S0 = sech(x-10); 

imp1 = sech(x - 1);
S1 = 0.5 * (imp1 + S0);

xmin = -1;
xmax = 10;
% find inflection point Si
range_x = x(x>=xmin & x<=xmax);
range_S1 = S1(x>=xmin & x<=xmax);   
MinPtsS1 = islocalmin(range_S1);
x_S1 = range_x(MinPtsS1);
y_S1 = range_S1(MinPtsS1);
fprintf ('x_S1 = %f , y_S1 = %+f .\n', [x_S1;y_S1]);

x_S1 = 5.500000 , y_S1 = +0.022215 .

% polynom approximation pv
[p_S1, ~, mu] = polyfit(x, S1, 26); 
pv_S1 = polyval (p_S1, x, [], mu); 

% find inflection point of polynom pv

range_x = x(x>=xmin & x<=xmax);
range_pvS1 = pv_S1(x>=xmin & x<=xmax);   
MinPts_pvS1 = islocalmin(range_pvS1);
x_pvS1 = range_x(MinPts_pvS1);
y_pvS1 = range_pvS1(MinPts_pvS1);
fprintf('x_pvS1 = %f , y_pvS1 = %+f .\n', [x_pvS1;y_pvS1]);

x_pvS1 = 5.500000 , y_pvS1 = +0.019493 .

% DIFF find inflection of diff2 Si, corresponding to inflection point of Si
dS1 = 20*diff (S1);
d2S1 = 10*diff(dS1);
range_x = x(x>=xmin & x<=xmax);
range_d2S1 = d2S1(x>=xmin & x<=xmax);   
MinPts_d2S1 = islocalmin(range_d2S1);
x_d2S1 = range_x(MinPts_d2S1);
y_d2S1 = range_d2S1(MinPts_d2S1);
fprintf('x_d2S1 = %f , y_d2S1 = %+f .\n', [x_d2S1;y_d2S1]);

x_d2S1 = 0.980000 , y_d2S1 = -0.039983 .
x_d2S1 = 5.480000 , y_d2S1 = +0.001776 .
x_d2S1 = 9.980000 , y_d2S1 = -0.039983 .

% DIFF find inflection of diff2 polynom Si, corresponding to inflection point of Si
dpS1 = 20*diff(pv_S1);
dp2S1 = 20*diff(dpS1);
range_x = x(x>=xmin & x<=xmax);
range_dp2S1 = dp2S1(x>=xmin & x<=xmax);   
MinPts_dp2S1 = islocalmin(range_dp2S1);
x_dp2S1 = range_x(MinPts_dp2S1);
y_dp2S1 = range_dp2S1(MinPts_dp2S1);
fprintf('x_dp2S1 = %f , y_dp2S1 = %+f .\n', [x_dp2S1;y_dp2S1]);

x_dp2S1 = 0.920000 , y_dp2S1 = -0.065128 .
x_dp2S1 = 4.380000 , y_dp2S1 = +0.002983 .
x_dp2S1 = 6.540000 , y_dp2S1 = +0.003272 .

Taking a look at that data:

plot(range_x,range_S1, range_x,range_pvS1, range_x,range_d2S1, range_x,range_dp2S1)
hold on
plot(x_S1,y_S1,'ok')
plot(x_pvS1,y_pvS1,'r+')
plot(x_d2S1,y_d2S1,'gd')
plot(x_dp2S1,y_dp2S1,'bx')

Are those the inflection points that you would expect to see?

[Anna_P]

@Stephen23 Thank you very much!!

This is so detailed and helpful!

[Paul]

The inflection points marked on the graph don't look like inflection points of the function that's plotted. Is the logic for determining inflection points correct?

[Anna_P]

@Paul Yes, thank you for noticing this! Actually, the polynomial is not plotted here, and those points mark its inflections, because the order is not enough. I increased the order of polynomial, and the minimum appeared in corresponding to min(S1) point

[Anna_P]

@Stephen23 and @Paul, maybe you could show an "easy" way to find zeros of 1-st derivative? Or there is a "100%" working function with gradient?

I've tried several ways, with fzero, roots and solve, but got confused and cannot grasp a concept so far..

[Paul]

The red "inflection S1" points are based on the original data, not the polynomial, and appear to be incorrect. So the fundamental issue is not with polynomial fit, but with the criteria for finding an inflection point.

Going to a polynomial order higher than 26 will likely cause even more problems in finding inflection points as compared to the original curve (at least in this example), and there likely is already problems with such a high order polynomial (can't say for sure because I did not plot pv_S1) in finding its inflection points as compared to the underlying data (at least for this example).

[Stephen23]

"maybe you could show an "easy" way to find zeros of 1-st derivative? Or there is a "100%" working function with gradient?"

Here are some approaches to find the inflection points of a set of data, which are zeros of the 2nd derivative. We assume that the data represents some reasonably well-behaved curve, e.g.:

X = linspace(-5, 15, 1001);
Y = 0.5 * (sech(X-1) + sech(X-10));

Method one: DIFF

The most straightforward approach since inflection points occur where the second derivative changes sign.

d2y = diff(Y,2);
idx = 1+find(diff(sign(d2y)) ~= 0);
X(idx)

ans = 1×4

0.1000 1.8800 9.1000 10.8800

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Y(idx)

ans = 1×4

0.3489 0.3542 0.3492 0.3539

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Method two: GRADIENT

d2y = gradient(gradient(Y, X), X);
idx = find(diff(sign(d2y)) ~= 0);
X(idx)

ans = 1×4

0.1000 1.8800 9.1000 10.8800

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Y(idx)

ans = 1×4

0.3489 0.3542 0.3492 0.3539

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Method 3: FZERO and INTERP1

Returns inflection points interpolated from the original X and Y data:

d2y = gradient(gradient(Y, X), X);
f_d2y = @(x) interp1(X, d2y, x, 'spline');

sc = find(diff(sign(d2y)) ~= 0);
Xi = nan(size(sc));
for k = 1:numel(sc)
    idx = sc(k);
    x1 = X(idx);
    x2 = X(idx + 1);
    Xi(k) = fzero(f_d2y, [x1, x2]);
end
Xi

Xi = 1×4

0.1184 1.8811 9.1189 10.8816

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Yi = interp1(X,Y,Xi)

Yi = 1×4

0.3536 0.3539 0.3539 0.3536

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Method 4: FZERO & SPLINE

Returns inflection points interpolated from the original X and Y data:

pp = spline(X, Y);
pp2 = fnder(pp,2); % Requires Curve Fitting Toolbox
breaks = pp.breaks;
Xj = [];
for k = 1:numel(breaks)-1
    segment = @(x) ppval(pp2, x);
    try
        zpt = fzero(segment, breaks([k,k+1]));
    catch
        continue
    end
    Xj(end+1) = zpt;
end
Xj

Xj = 1×4

0.1188 1.8807 9.1193 10.8812

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Yj = interp1(X,Y,Xj)

Yj = 1×4

0.3536 0.3540 0.3540 0.3536

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Lets take a look at the inflection points from method 4:

plot(X,Y)
hold on
plot(Xj,Yj,'rx')

●

[Anna_P]

@Stephen23, thank you very much, this is really a lot!!

And yes, I was looking for those new inflection points)

But, the problem is that I need the min point of polynomial, around x = 5, which is not presented in all 4 methods..

[Stephen23]

"... I need the min point of polynomial, around x = 5, which is not presented in all 4 methods.. "

X = linspace(-5, 15, 1001);
Y = 0.5 * (sech(X-1) + sech(X-10));
[p_S1, ~, mu] = polyfit(X, Y, 26); 
pv_S1 = polyval(p_S1, X, [], mu);
dp_S1 = polyder(p_S1);
dpoly = @(xq) polyval(dp_S1, (xq-mu(1))./mu(2));
x_min = fzero(dpoly, 5)

x_min = 5.5017

y_min = interp1(X,Y,x_min)

y_min = 0.0222

[Anna_P]

@Stephen23 THANK YOU!)

I'm not desperate anymore! Going to explore and study HELP section with more clarity now!