function interpolation and root finding

Hello community,

i have 37 points in in the interval [-.5 .5]. and i'm looking for a minimum. I have written a script which uses polyfit for the polynomial and polyder with fzero to find the minium of the interpolated function. but I get some warnings..

points =[...
   0.500000000000000   0.500000000000000;...
   0.498097349045873   0.494452298903742;...
   0.492403876506104   0.477882019169956;...
   0.482962913144534   0.450598639000461;...
   0.469846310392954   0.413355783194506;...
   0.453153893518325   0.367259042487410;...
   0.433012701892219   0.313952225075596;...
   0.409576022144496   0.255588017440073;...
   0.383022221559489   0.194787975761779;...
   0.353553390593274   0.134485572953294;...
   0.321393804843270   0.077602649602957;...
   0.286788218175523   0.026729010250023;...
   0.250000000000000  -0.016230934343725;...
   0.211309130870350  -0.050225749043852;...
   0.171010071662834  -0.075118583665429;...
   0.129409522551260  -0.091527009978884;...
   0.086824088833465  -0.100617488371580;...
   0.043577871373829  -0.103815701636865;...
   0.000000000000000  -0.102599285091917;...
  -0.043577871373829  -0.098307643482893;...
  -0.086824088833465  -0.092067023864118;...
  -0.129409522551260  -0.084735622572768;...
  -0.171010071662834  -0.076929476902457;...
  -0.211309130870350  -0.069040079308481;...
  -0.250000000000000  -0.061295542051084;...
  -0.286788218175523  -0.053796351887970;...
  -0.321393804843270  -0.046573711545795;...
  -0.353553390593274  -0.039619979531691;...
  -0.383022221559489  -0.032931563548278;...
  -0.409576022144496  -0.026527753568514;...
  -0.433012701892219  -0.020473888300108;...
  -0.453153893518325  -0.014883984708655;...
  -0.469846310392954  -0.009918994629169;...
  -0.482962913144534  -0.005769046336078;...
  -0.492403876506104  -0.002628877217708;...
  -0.498097349045873  -0.000667477816139;...
  -0.500000000000000                   0];

I have written a script which give me the interpolated polynomial, but i get the Warning: Polynomial is badly conditioned. Add points with distinct X values or reduce the degree

I tested the function polyfit also with: [P,S,MU] = polyfit(X,Y,N) but i did not helped me.

Maybe the condition number of the Vandermod matrix is not low enough for the interpolation.

Does somebody know some alternatives? Maybe with an Aitken's Algorithm or Lagrange Interpolation ?

It would be nice to find a numerical accurate solution to interpolate to a polynomial degree of numel(x) -1.

Here is the code is used:

x=points(:,1);
y=points(:,2);
[P, S, mu]=polyfit(x,y,numel(x)-1);
[~,i_xmin]=min(x);
x=fzero(@(x) polyval(polyder(P),x,S,mu),x(i_xmin));
f=polyval(P,x,S,mu);

Best regards

emjay

Answer 1

[Star Strider]

Try this —

points =  [0.500000000000000   0.500000000000000
           0.498097349045873   0.494452298903742
           0.492403876506104   0.477882019169956
           0.482962913144534   0.450598639000461
           0.469846310392954   0.413355783194506
           0.453153893518325   0.367259042487410
           0.433012701892219   0.313952225075596
           0.409576022144496   0.255588017440073
           0.383022221559489   0.194787975761779
           0.353553390593274   0.134485572953294
           0.321393804843270   0.077602649602957
           0.286788218175523   0.026729010250023
           0.250000000000000  -0.016230934343725
           0.211309130870350  -0.050225749043852
           0.171010071662834  -0.075118583665429
           0.129409522551260  -0.091527009978884
           0.086824088833465  -0.100617488371580
           0.043577871373829  -0.103815701636865
           0.000000000000000  -0.102599285091917
          -0.043577871373829  -0.098307643482893
          -0.086824088833465  -0.092067023864118
          -0.129409522551260  -0.084735622572768
          -0.171010071662834  -0.076929476902457
          -0.211309130870350  -0.069040079308481
          -0.250000000000000  -0.061295542051084
          -0.286788218175523  -0.053796351887970
          -0.321393804843270  -0.046573711545795
          -0.353553390593274  -0.039619979531691
          -0.383022221559489  -0.032931563548278
          -0.409576022144496  -0.026527753568514
          -0.433012701892219  -0.020473888300108
          -0.453153893518325  -0.014883984708655
          -0.469846310392954  -0.009918994629169
          -0.482962913144534  -0.005769046336078
          -0.492403876506104  -0.002628877217708
          -0.498097349045873  -0.000667477816139
          -0.500000000000000                   0];


p = polyfit(points(:,1), points(:,2), 7)                                            % Fit 7th-Degree Polynomial

p = 1×8

-7.5119 -6.2764 0.8450 4.0360 2.5354 0.7934 -0.0694 -0.1026

dp = polyder(p)                                                                     % Calculate Derivative

dp = 1×7

-52.5832 -37.6587 4.2248 16.1441 7.6062 1.5868 -0.0694

dpr = roots(dp)                                                                     % Roots Of Derivative Polynomial

dpr = 
   0.6763 + 0.0000i
  -0.4904 + 0.2286i
  -0.4904 - 0.2286i
  -0.2242 + 0.3621i
  -0.2242 - 0.3621i
   0.0368 + 0.0000i

dpr_real = dpr(imag(dpr)==0)                                                        % Return Real Roots

dpr_real = 2×1

0.6763 0.0368

dpr_in_interval = dpr_real(dpr_real>=min(points(:,1)) & dpr_real<=max(points(:,1))) % Return Real Roots Within Interval

dpr_in_interval = 0.0368

dpryval = polyval(p, dpr_in_interval)                                               % Evaluate Original Polynomial At That Value

dpryval = -0.1039

dpv = polyval(dp,points(:,1));
v = polyval(p, points(:,1));

figure
plot(points(:,1), points(:,2),'.')
hold on
plot(points(:,1), v)
plot(dpr_in_interval, dpryval, 'sg')
hold off
grid
legend('Original Data',' Fitted Polynomial', 'Calculated Minimum', 'Location','best')

EDIT — (25 Jun 2021 at 17:55)

It is not necessary to use fzero with polyder. Just use the roots function to find the roots of the derivative of the polynomial. It is likely more accurate than fzero, and likely much more efficient.

If you absolutely must use fzero:

pointsmin = fzero(@(x)polyval(dp,x), -0.5)

pointsmin = 0.0368

With the identical result.

.

[Marko]

Hello Star Strider, thank you for the information to use roots, instead of fzero.

Do you know if there is possible increase the polynomial order to numel(x)-1 without this warning:

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

> In polyfit (line 83)

I belive that polyfit is not capable of this high order, because of the high condition number of the Vandermonde matrix. So an another algorithm is needed.

[Star Strider]

My pleasure!

The problem with the polynomial order is that any order more than 7 will not produce a more accurate fit. That is simply the nature of polynomilal approximations. (Even increasing the order from 5 originally to 7 in the edited version did not change the result beyond a few non-signifcant decimal places.) I also douobt that another algorithm is necessary, simply because a higer order will not increase the accuracy of the fit.

It is relatively easy to do that experiment.

Example —

points =  [0.500000000000000   0.500000000000000
           0.498097349045873   0.494452298903742
           0.492403876506104   0.477882019169956
           0.482962913144534   0.450598639000461
           0.469846310392954   0.413355783194506
           0.453153893518325   0.367259042487410
           0.433012701892219   0.313952225075596
           0.409576022144496   0.255588017440073
           0.383022221559489   0.194787975761779
           0.353553390593274   0.134485572953294
           0.321393804843270   0.077602649602957
           0.286788218175523   0.026729010250023
           0.250000000000000  -0.016230934343725
           0.211309130870350  -0.050225749043852
           0.171010071662834  -0.075118583665429
           0.129409522551260  -0.091527009978884
           0.086824088833465  -0.100617488371580
           0.043577871373829  -0.103815701636865
           0.000000000000000  -0.102599285091917
          -0.043577871373829  -0.098307643482893
          -0.086824088833465  -0.092067023864118
          -0.129409522551260  -0.084735622572768
          -0.171010071662834  -0.076929476902457
          -0.211309130870350  -0.069040079308481
          -0.250000000000000  -0.061295542051084
          -0.286788218175523  -0.053796351887970
          -0.321393804843270  -0.046573711545795
          -0.353553390593274  -0.039619979531691
          -0.383022221559489  -0.032931563548278
          -0.409576022144496  -0.026527753568514
          -0.433012701892219  -0.020473888300108
          -0.453153893518325  -0.014883984708655
          -0.469846310392954  -0.009918994629169
          -0.482962913144534  -0.005769046336078
          -0.492403876506104  -0.002628877217708
          -0.498097349045873  -0.000667477816139
          -0.500000000000000                   0];

size(points,1)

ans = 37

dpr_in_interval = num2cell(NaN(size(points,1)-1,1));
for k = 1:size(points,1)-1
    order(k,:) = k;
    p = polyfit(points(:,1), points(:,2), k);                                                   % Fit 7th-Degree Polynomial
    dp = polyder(p);                                                                            % Calculate Derivative
    dpr = roots(dp);                                                                            % Roots Of Derivative Polynomial
    dpr_real = dpr(imag(dpr)==0);                                                               % Return Real Roots
    select = dpr_real(dpr_real>=min(points(:,1)) & dpr_real<=max(points(:,1)))                 % Return Real Roots Within Interval
    if ~isempty(select)
        dpr_in_interval{k,:} = select;
    end
end

select =

  0×1 empty double column vector

select = -0.1325

select = 2×1

-0.4525 0.0113

select = 0.0156

select = 0.0351

select = 0.0416

select = 0.0368

select = 0.0368

select = 0.0362

select = 0.0359

select = 0.0359

select = 0.0359

select = 0.0359

select = 0.0360

select = 0.0359

select = 0.0359

select = 0.0359

select = 0.0360

select = 0.0359

select = 0.0359

select = 0.0359

select = 0.0359

select = 0.0360

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0360

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

select = 0.0359

Results = table(order,dpr_in_interval)

Results = 36×2 table
    order    dpr_in_interval
    _____    _______________

      1       {[     NaN]}  
      2       {[ -0.1325]}  
      3       {2×1 double}  
      4       {[  0.0156]}  
      5       {[  0.0351]}  
      6       {[  0.0416]}  
      7       {[  0.0368]}  
      8       {[  0.0368]}  
      9       {[  0.0362]}  
     10       {[  0.0359]}  
     11       {[  0.0359]}  
     12       {[  0.0359]}  
     13       {[  0.0359]}  
     14       {[  0.0360]}  
     15       {[  0.0359]}  
     16       {[  0.0359]}  

So the result does not change significantly (and then only in the third decimal place) beyond the 7th-order polynomial.

.