? polyfit.patch Index: polyfit.m =================================================================== RCS file: /cvs/octave/scripts/polynomial/polyfit.m,v retrieving revision 1.23 diff -u -r1.23 polyfit.m --- polyfit.m 23 Aug 2005 18:38:28 -0000 1.23 +++ polyfit.m 9 Aug 2006 19:24:50 -0000 @@ -78,10 +78,16 @@ x = reshape (x, l, 1); y = reshape (y, l, 1); + ## scale the x values + xmax = max(abs(x)); + x = x/xmax; + + X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0)); p = X \ y; + if (nargout > 1) yf = X*p; @@ -99,8 +105,29 @@ endif - ## Return value should be a row vector. - - p = p.'; + ## Scale the polynomial back and return a row vector + power = length(p)-1:-1:0; + power = (xmax .^ power)(:); + p = (p ./ power).'; endfunction + +## test simple polynomial; here we check the returned polynomial +%!test +%! x = [1,2,3]; +%! p = poly(x); +%! y = polyval(p,x); +%! pnew = polyfit([0,1,2,3]', [-6,0,0,0]', 3); +%! assert(p, pnew, 10^-6); + +## test difficult case where scaling is really needed +## here we check the values after the evaluation of the new polynomial +%!test +%! x = [ -1196.4, -1195.2, -1194, -1192.8, -1191.6, -1190.4, -1189.2, -1188, \ +%! -1186.8, -1185.6, -1184.4, -1183.2, -1182 ]; +%! y = [ 315571.7086, 315575.9618, 315579.4195, 315582.6206, 315585.4966, \ +%! 315588.3172, 315590.9326, 315593.5934, 315596.0455, 315598.4201, \ +%! 315600.7143, 315602.9508, 315605.1765 ]; +%! p = polyfit(x,y,10); +%! ynew = polyval(p,x); +%! assert(norm(y-ynew), 0.08, 0.01);