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

archlengthintegralMATLABpolyderpolyfitpolynomial

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?

Best Answer

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.
Related Question