How to find coefficients of polynomial equation?

Question

Given two points in the x, y plane:

x, f(x)
1, 3
2, 5

I can interpolate them using Lagrange and find f(1.5), which result in 4. Thinking a little I managed to find a way to discover the coefficients of the equation:

void l1Coefficients(const vector<double> &x, const vector<double> &y) {

    double a0 = y[0]/(x[0]-x[1]);
    double a1 = y[1]/(x[1]-x[0]);

    double b0 = (-x[1]*y[0])/(x[0]-x[1]);
    double b1 = (-x[0]*y[1])/(x[1]-x[0]);

    double a = a0 + a1;
    double b = b0 + b1;

    cout << "P1(x) = " << a << "x +" << b << endl;
}

That gives me P1(x) = 2x +1.

Thinking a little more I was able to extend that to 2nd order equations. So, given the points:

1, 1
2, 4
3, 9

I found the equation P2(x) = 1x^2 +0x +0 with the following:

void l2Coefficients(const vector<double> &x, const vector<double> &y) {

    double a0 =              y[0] / ((x[0]-x[1])*(x[0]-x[2]));
    double a1 =              y[1] / ((x[1]-x[0])*(x[1]-x[2]));
    double a2 =              y[2] / ((x[2]-x[0])*(x[2]-x[1]));

    double b0 = -(x[1]+x[2])*y[0] / ((x[0]-x[1])*(x[0]-x[2]));
    double b1 = -(x[0]+x[2])*y[1] / ((x[1]-x[0])*(x[1]-x[2]));
    double b2 = -(x[0]+x[1])*y[2] / ((x[2]-x[0])*(x[2]-x[1]));

    double c0 =  (x[1]*x[2])*y[0] / ((x[0]-x[1])*(x[0]-x[2]));
    double c1 =  (x[0]*x[2])*y[1] / ((x[1]-x[0])*(x[1]-x[2]));
    double c2 =  (x[0]*x[1])*y[2] / ((x[2]-x[0])*(x[2]-x[1]));

    double a = a0 + a1 + a2;
    double b = b0 + b1 + b2;
    double c = c0 + c1 + c2;

    cout << "P2(x) = " << a << "x^2 +" << b << "x +" << c << endl;
}

Working hard I actually was able to find the coefficients for equations of order up to 4th.

How to find the coefficients of order n equations?

Answer 1

Antonio already pointed you to the formula for Lagrange interpolation.

There is a slightly more general approach that you can take that might be useful: suppose your data consists of $(x_i, f(x_i))$ pairs, and you want to fit a degree $N$ polynomial $\sum_{i=0}^N \alpha_i x^i.$ You can set up a linear system to solve for the coefficients:

$$\left[\begin{array}{ccccc}1 & x_0 & x_0^2 & \cdots & x_0^N\\1 & x_1 & x_1^2 & \cdots & x_1^N\\\vdots & \vdots & \vdots & & \vdots\end{array}\right]\left[\begin{array}{c}\alpha_0\\\alpha_i\\ \vdots\\\alpha_N\end{array}\right] = \left[\begin{array}{c}f(x_0)\\f(x_1)\\\vdots \\f(x_N)\end{array}\right].$$

When you have $N+1$ pairs, you get a square system and solving for the $\alpha$s gives you exactly Lagrange interpolation. If you have more than $N+1$ pairs, the system is overdetermined, but you can solve it in the least squares sense, which gives you the best-fitting degree $N$ polynomial to your data.