function interpolation and root finding
조회 수: 11 (최근 30일)
이전 댓글 표시
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
채택된 답변
Star Strider
2021년 7월 25일
편집: Star Strider
2021년 7월 25일
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
Star Strider
2021년 7월 25일
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.
.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Polynomials에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!