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

As part of my algorithm, I need the arc length of a cubic function. The cubic function is generated by polyfit based off of my input data (in this case, the centroids of words in a line of skewed text).

For some (but not all) of my data, polyfit returns a warning stating: "Polynomial is badly conditioned." Based on the polyfit documentation I want to use the centering and scaling factors returned in MU to deal with this source of error. However, based on some of the methods used here, the best way to get arc length would be to differentiate my cubic function and then use the integral function.

Is there a way that I can easily rework the centering and scaling factors back into the polynomial coefficients so that I can pass them to polyder? Or will it be necessary to manually solve for the new coefficient values?

Answer 1

[John D'Errico]

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.

[Alex Burden]

Thanks very much for the quick response, I haven't really used Matlab's symbolic tools much but this makes sense. I assume that since my polynomial will always be cubic I can just used the final result you reached, and will only need to carry out the expansion should my algorithm necessitate a different order function.

[John D'Errico]

If it is always a cubic, then yes, simplest is just to use those coefficients as the derivative.

Then integrate the square of that function over the domain of your curve. So it might look like...

integral(@(x) sqrt(1 + polyval(dpdx,x).^2),x1,x2)

where dpdx is the vector of derivative coefficients I derived using the symbolic toolbox.

Were you to change the order pf the polynomial in the future, I'd use the alternative chain rule derivation. Either will work as long as you are careful.

Finally, test it out!!!!!!! Verify that I did not get a typo in there, and that you got it all right. ALWAYS test your results to make sure they make sense.