function interpolation and root finding
12 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Marko
el 25 de Jul. de 2021
Comentada: Star Strider
el 25 de Jul. de 2021
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
0 comentarios
Respuesta aceptada
Star Strider
el 25 de Jul. de 2021
Editada: Star Strider
el 25 de Jul. de 2021
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
dp = polyder(p) % Calculate Derivative
dpr = roots(dp) % Roots Of Derivative Polynomial
dpr_real = dpr(imag(dpr)==0) % Return Real Roots
dpr_in_interval = dpr_real(dpr_real>=min(points(:,1)) & dpr_real<=max(points(:,1))) % Return Real Roots Within Interval
dpryval = polyval(p, dpr_in_interval) % Evaluate Original Polynomial At That Value
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)
With the identical result.
.
2 comentarios
Star Strider
el 25 de Jul. de 2021
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)
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
Results = table(order,dpr_in_interval)
So the result does not change significantly (and then only in the third decimal place) beyond the 7th-order polynomial.
.
Más respuestas (0)
Ver también
Categorías
Más información sobre Polynomials en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!