Symbolic polynomial interpolation

Question

I'm trying to create polynomials from some symbolic points to discretize derivations. This means I'm having data like $(a, \phi(a)),\ (b, \phi(b) ) $and $(c, \phi(c))$ and want to fit a second order polynomial through these points.

I tried Matlabs polyfit but it can't handle syms. I also had a look at Maxima, but couldn't find any function there neither.

I'd be really happy if someone could point out how to do this in a CAS.

Answer 1

It annoys me when Matlab functions aren't overloaded for symbolic calculation. Looking at the code for polyfit (edit polyfit in the command window) it seems that doing what you need may be quite easy. I comes down to constructing a Vandermonde matrix and solving a least squares problem. The existing code can be adapted for symbolic math.

In fact, I couldn't resist and gave it try. You can find my code for polyfitsym. on my GitHub.

For an example, in your three point, second-order polynomial case, let $\phi(x) = \sin(x)$. First we'll solve for three arbitrary symbolic points $x_1 = a$, $x_2 = b$, and $x_3 = c$. Later we'll explicitly set $a = -3$, $b = -1$, and $c = 1$.

x = sym('x',[3 1]);   % Vector of three symbolic data points
y = sin(x);           % phi(x), apply function to vector, see help for other forms
p = polyfitsym(x,y,2) % Return coefficients for second order polynomial

which returns the symbolic result

p =

[ -(x1*sin(x2)-x2*sin(x1)-x1*sin(x3)+x3*sin(x1)+x2*sin(x3)-x3*sin(x2))/((x1-x2)*(x1*x2-x1*x3-x2*x3+x3^2)), ...
   (x1^2*sin(x2)-x2^2*sin(x1)-x1^2*sin(x3)+x3^2*sin(x1)+x2^2*sin(x3)-x3^2*sin(x2))/((x1-x2)*(x1*x2-x1*x3-x2*x3+x3^2)), ...
  -(x1*x2^2*sin(x3)-x1*x3^2*sin(x2)+x2*x3^2*sin(x1)-x1^2*x2*sin(x3)+x1^2*x3*sin(x2)-x2^2*x3*sin(x1))/((x1-x2)*(x1*x2-x1*x3-x2*x3+x3^2)) ]

This can be simplified a bit with simple or simplify or you can plug in numeric values with subs:

abc = [-3;-1;1];        % Three numeric data points
p2 = subs(p,x,sym(abc)) % Cast points to symbolic to keeps coefficients symbolic

which gives (still symbolic – use double to convert to floating point if you like)

p2 =

[ (3*sin(1))/8 - sin(3)/8, sin(1), sin(3)/8 - (3*sin(1))/8]

You might wish to plot this and see how good the fit is over a range:

x2 = -2*pi:0.2:2*pi;
f = polyval(double(p2),x2);  % polyval also only supports floating point
plot(x2,sin(x2),'o',x2,f,'-',abc,sin(abc),'*');
axis([x2(1) x2(end)  -2  2]);
xlabel('x'); ylabel('y(x)'); legend('sin(x)','fit','data'); legend boxoff;

which results in a figure like this:

There is another example in the help for the function. Since I don't know what exactly you're trying to do, the polyfitsym function is pretty general and works very much like polyfit. Feel free to modify, improve, and use it as you need. I can't guarantee that it won't have numeric issues (for some test cases all of the coefficients were Inf for some reason) or that it will be able to handle everything. It uses sym/linsolve rather that \ (mldivide) to solve the least squares linear system because it returns warnings if the the problem is ill-conditioned.