Symbolic polynomial interpolation

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.

• For three points, why not do it by hand symbolically? The Lagrange interpolation isn't hard for that. – Daniel Fischer Jul 28 '13 at 11:06
• Yeah, the thing is that I need to to it for higher orders as well, so by hand gets really frustrating. – Paul Jul 28 '13 at 11:08
• Ah, yes, for more points, it quickly becomes cumbersome by hand. You might have some success checking the docs for "Lagrange" or "interpolation". – Daniel Fischer Jul 28 '13 at 11:11
• What are your points? – NoChance Oct 16 '15 at 22:31

Mathematica can find an interpolating polynomial symbolically. For example,

InterpolatingPolynomial[{{a, w}, {b, x}, {c, y}, {d, z}}, t]

$$(t-a) \left((t-b) \left(\frac{(t-c) \left(\frac{\frac{z-y}{d-c}-\frac{y-x}{c-b}}{d-b}-\frac{\frac{y-x}{c-b}-\frac{x-w}{b-a}}{c-a}\right)}{d-a}+\frac{\frac{y-x}{c-b}-\frac{x-w}{b-a}}{c-a}\right)+\frac{x-w}{b-a}\right)+w.$$

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.

• Thanks a lot for the great example! – Paul Oct 10 '13 at 19:44
• @Paul: Did this or another answer help you with your problem? Consider accepting an answer if so. – horchler Apr 7 '16 at 21:01