MATLAB: Poly and polyval for large vectors…

polypolyval

Hello,

I need to get the equation of the polynomial through 50 points… a 49th order should make a perfect fit. This doesn't do it:

(load x and y vectors with unique monotonically increasing values…then:)

            [p,S,mu] = polyfit(x,y,49);
            yp = polyval(p,x);

There is some advice to use "centering and scaling" from the warning message. Using HELP I see:

    [P,S,MU] = polyfit(X,Y,N) finds the coefficients of a polynomial in
    XHAT = (X-MU(1))/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). This
    centering and scaling transformation improves the numerical properties
    of both the polynomial and the fitting algorithm.

I'm not sure how to follow that advice. Is the an example where a large number of points is used with a poly that should capture them perfectly?

Thanks for helping.

Best Answer

The documentation of polyfit has an example of how to use the centering and scaling option. You need to use that returned mu in polyval for it to make sense:

[p, ~, mu] = polyfit(x, y, 49);
yp = polyval(p, x, [], mu);

This is way outside my area of expertise but I believe that fitting a 49 degree polynomial to 50 points is not an appropriate way to fit data. You probably would be better off fitting splines.

Related Solutions

MATLAB: How does the automatic scaling factor from the polyfit function interact with polyder

Sigh. First of all, what Star Strider showed you is NOT necessarily the best way to solve the problem, since trapz is simply a trapezoidal rule. Don't use trapezoidal rule when there is no need to do so, since it introduces errors into the result that can be significant.

Yes, you can use integral to solve the arclength problem, arguably a better choice than trapz here. But differentiation of a polynomial is pretty trivial. NO. You do NOT need to redo the estimation.

From the help for polyfit...

[P,S,MU] = polyfit(X,Y,N) finds the coefficients of a polynomial in
    XHAT = (X-MU(1))/MU(2)

So in the vector P, we have the coefficients of a polynomial in the form

P(1)*X^3 + P(2)*X^2 + P(3)*X + P(4)

where X = (x-mu(1))/mu(2)

Simplest is often just to apply brute force. The sledge hammer will suffice here.

syms P1 P2 P3 P4 x mu1 mu2 X
P = P1*X^3 + P2*X^2 + P3*X + P4;
dpdx = expand(diff(subs(P,X,(x-mu1)/mu2)))
dpdx =
P3/mu2 + (3*P1*mu1^2)/mu2^3 + (3*P1*x^2)/mu2^3 - (2*P2*mu1)/mu2^2 + (2*P2*x)/mu2^2 - (6*P1*mu1*x)/mu2^3

Now, extract the coefficient of the derivative as:

flip(coeffs(dpdx,x))
ans =
[ (3*P1)/mu2^3, (2*P2)/mu2^2 - (6*P1*mu1)/mu2^3, P3/mu2 + (3*P1*mu1^2)/mu2^3 - (2*P2*mu1)/mu2^2]

The result is the list of coefficients of the quadratic (derivative) polynomial, with highest order term first.

But you can use that in a call to integral.

Could you have done it without brute force? Of course. Basic calculus. Given a function P(X), where we have X=(x-mu(1))/mu(2), what is the derivative of P here?

P((x-mu(1))/mu(2))

So chain rule tells us that

dP/dx = dP/dX * DX/dx

And you do have X as a function of x, mu(1), and mu(2), so we know that

dX/dx = 1/mu(2)

So just apply polyder to P directly. Evaluate it at X=(x-mu(1))/mu(2), then divide that derivative by mu(2).

In either case, just use integral to integrate the kernel

sqrt(1 + dPdX(x).^2)

over the support of your function.

Best Answer

Related Solutions

MATLAB: How does the automatic scaling factor from the polyfit function interact with polyder

Related Question