Re: polynomial plot for spline interpolation

"Otavio " <otaviofm@xxxxxxxxx> wrote in message <ggerg2$n51$1@xxxxxxxxxxxxxxxxxx>...

Use mkpp and ppval was really a good hint, but I still can't find out how to convert all the segments in one. Instead, I tried to plot each one using 'hold', but it didn't work. I wrote the code as follows:

for i=1:20
cc = [a(i) b(i) c(i) d(i)];
pp = mkpp([x(i) x(i+1)],cc);
xx = x(i):0.1:x(i+1);
plot(xx,ppval(pp,xx))
hold on;
end

but the resulting plot was very far from what I expected.. The curve segments weren't even conected with each other. Any suggestions..?

Ohhh man, I found out what was going wrong in my plot! The plot wasn't really 'a(i)x^3 + b(i)x^2 + c(i)x + d1' ! Instead, I should do a(i)(x-x(i))^3 + b(i)(x-x(i))^2 + c(i)(x-x(i)) + d1. Too bad, since it means I can not use mkpp anymore... Or can I..? Now I'm really confused..

You are gradually getting there. Yes, the pp form
must be a set of polynomials in a specific form.
They will be evaluated, first subtracting the left
hand endpoint of the interval that the polynomial
lives in. This is done for a very good reason.

Suppose I have a spline that lived over the set of
knots 0:1:10. Now, suppose I had exactly the
same function with the same set of values it must
pass through, but assume that it lives on the knots
1e6+(0:1:10).

The polynomials in the second case will force us to
square and cube numbers on the order of 1e6, then
adding and subtracting the results. This will lead to
numerical garbage. So instead, we always create
the polynomials to live on the domains

[0, x(i+1) - x(i)]

In your case, you created a set of polynomials that
live on the wrong domains. This is still survivable.
Suppose we have the polynomial with coefficients

P = [a b c d]

representing the symbolic polynomial

a*x^3 + b*x^2 + c*x + d

Assume that P has been created to live on the
interval [x1,x2], but we need to convert P into
another polynomial that lives on the interval
[0,x2-x1], yet produces the same predictions
over these intervals.

Q = a*x^3 + (b + 3*a*x1)*x^2 + ...
(c+2*b*x1+3*a*x1^2)*x + ...
(d+c*x1+b*x1^2+a*x1^3)

In terms of the original polynomial P, this is
most easily written with just a matrix multiply
to create the new polynomial coefficients.

Q = P*[1 3*x1 3*x1^2 x1^3; ...
0 1 2*x1 x1^2; 0 0 1 x1;0 0 0 1];

Test it out on a random polynomial.

P = rand(1,4)
P =
0.94478 0.71447 0.67922 0.95938

x1 = 5;

Q = P*[1 3*x1 3*x1^2 x1^3; ...
0 1 2*x1 x1^2; 0 0 1 x1;0 0 0 1]

Q =
0.94478 14.886 78.683 140.32

evaluate P at a set of points starting from x = 5.

polyval(P,5:.1:5.5)
ans =
140.32 148.33 156.65 165.29 174.23 183.5

And evaluate Q at the corresponding points, from
zero up.

polyval(Q,0:.1:.5)
ans =
140.32 148.33 156.65 165.29 174.23 183.5

Yup, it worked. Must've got lucky.

HTH,
John
.

Relevant Pages

  • Polynomial Discriminant
    ... determinant of a matrix which consists of the polynomial coefficients ... the derivative of that polynoimial shifted according to a pattern. ... But I want to have a clue for polynomials of general degree n. ...
    (sci.math)
  • Re: Generate function
    ... Jens Gruschel wrote: ... >you to take a look at polynomials. ... Instead of the polynomial coefficients you get the ... >http://www.wikipedia.org (search for polynomial, spline, or interpolation). ...
    (borland.public.delphi.thirdpartytools.general)
  • Re: A motion profile problem or Mathcad sucks. Will some other CAS solve this?
    ... Bob Walton wrote: ... Correction to the Sylvester matrix description: ... matrix is an square matrix with the coefficients of the two polynomials arranged as follows (generally, the polynomial coefficients are zero-padded and shifted right on each successive row: ...
    (sci.math.symbolic)
  • Re: A motion profile problem or Mathcad sucks. Will some other CAS solve this?
    ... Bob Walton wrote: ... Correction to the Sylvester matrix description: ... matrix is an square matrix with the coefficients of the two polynomials arranged as follows (generally, the polynomial coefficients are zero-padded and shifted right on each successive row: ...
    (sci.math.symbolic)