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? Where

Pn(x) = c_2x^2 + c_1x^1 + c_0x^0 + ...

Answer 1

It's a simple linear algebra problem.

We have a set of N samples of the form x_k -> f(x_k) and we know the general form of function f(x), which is:

f(x) = c₀x⁰ + c₁x¹ + ... + c_N-1x^N-1

We want to find the coefficients c₀ ... c_N-1. To achieve that, we build a system of N equations of the form:

c₀x_k⁰ + c₁x_k¹ + ... + c_N-1x_k^N-1 = f(x_k)

where k is the sample number. Since x_k and f(x_k) are constants rather than variables, we have a linear system of equations.

Expressed in terms of linear algebra, we have to solve:

Ac = b

where A is a Vandermonde matrix of powers of x and b is a vector of f(x_k) values.

To solve such a system, you need a linear algebra library, such as Eigen. See here for example code.

The only thing that can go wrong with such an approach is the system of linear equations being under-determined, which will happen if your N samples can be fit with with a polynomial of degree less than N-1. In such a case you can still solve this system with Moore-Penrose pseudo inverse like this:

c = pinv(A)*b

Unfortunately, Eigen doesn't have a pinv() implementation, though it's pretty easy to code it by yourself in terms of Singular Value Decomposition (SVD).

Answer 2

I created a naive implementation of the matrix solution:

#include <iostream>
#include <vector>
#include <stdexcept>

class Matrix
{

private:

    class RowIterator
    {
    public:
        RowIterator(Matrix* mat, int rowNum) :_mat(mat), _rowNum(rowNum) {}
        double& operator[] (int colNum) { return _mat->_data[_rowNum*_mat->_sizeX + colNum]; }
    private:
        Matrix* _mat;
        int _rowNum;
    };

    int _sizeY, _sizeX;
    std::vector<double> _data;

public:

    Matrix(int sizeY, int sizeX) : _sizeY(sizeY), _sizeX(sizeX), _data(_sizeY*_sizeX){}
    Matrix(std::vector<std::vector<double> > initList) : _sizeY(initList.size()), _sizeX(_sizeY>0 ? initList.begin()->size() : 0), _data()
    { 
        _data.reserve(_sizeY*_sizeX);
        for (const std::vector<double>& list : initList)
        {
            _data.insert(_data.end(), list.begin(), list.end());
        }
    }

    RowIterator operator[] (int rowNum) { return RowIterator(this, rowNum); }

    int getSize() { return _sizeX*_sizeY; }
    int getSizeX() { return _sizeX; }
    int getSizeY() { return _sizeY; }

    Matrix reduce(int rowNum, int colNum)
    {
        Matrix mat(_sizeY-1, _sizeX-1);
        int rowRem = 0;
        for (int y = 0; y < _sizeY; y++)
        {
            if (rowNum == y)
            {
                rowRem = 1;
                continue;
            }
            int colRem = 0;
            for (int x = 0; x < _sizeX; x++)
            {
                if (colNum == x)
                {
                    colRem = 1;
                    continue;
                }
                mat[y - rowRem][x - colRem] = (*this)[y][x];
            }
        }
        return mat;
    }

    Matrix replaceCol(int colNum, std::vector<double> newCol)
    {
        Matrix mat = *this;
        for (int y = 0; y < _sizeY; y++)
        {
            mat[y][colNum] = newCol[y];
        }
        return mat;
    }

};

double solveMatrix(Matrix mat)
{
    if (mat.getSizeX() != mat.getSizeY()) throw std::invalid_argument("Not square matrix");
    if (mat.getSize() > 1)
    {
        double sum = 0.0;
        int sign = 1;
        for (int x = 0; x < mat.getSizeX(); x++)
        {
            sum += sign * mat[0][x] * solveMatrix(mat.reduce(0, x));
            sign = -sign;
        }
        return sum;
    }

    return mat[0][0];
}

std::vector<double> solveEq(std::vector< std::pair<double, double> > points)
{
    std::vector<std::vector<double> > xes(points.size());
    for (int i = 0; i<points.size(); i++)
    {
        xes[i].push_back(1);
        for (int j = 1; j<points.size(); j++)
        {
            xes[i].push_back(xes[i].back() * points[i].first);
        }
    }

    Matrix mat(xes);

    std::vector<double> ys(points.size());

    for (int i = 0; i < points.size(); i++)
    {
        ys[i] = points[i].second;
    }

    double w = solveMatrix(mat);

    std::vector<double> result(points.size(), 0.0);

    if(w!=0)
        for (int i = 0; i < ys.size(); i++)
        {
            result[i] = solveMatrix(mat.replaceCol(i, ys));
            result[i] /= w;
        }

    return result;
}

void printCoe(std::vector<double> coe)
{
    std::cout << "f(x)=";
    bool notFirstSign = false;
    for (int i = coe.size() - 1; i >= 0; i--)
    {
        if (coe[i] != 0.0)
        {
            if (coe[i] >= 0.0 && notFirstSign)
                std::cout << "+";
            notFirstSign = true;
            if (coe[i] != 1.0)
                if (coe[i] == -1.0)
                    std::cout << "-";
                else
                    std::cout << coe[i];
            if (i == 1)
                std::cout << "x";
            if (i>1)
                std::cout << "x^" << i;
        }
    }
    std::cout << std::endl;
}

int main()
{
    std::vector< std::pair<double, double> > points1 = { {3,31}, {6,94}, {4,48}, {0,4} };
    std::vector<double> coe = solveEq(points1);
    printCoe(coe);

    std::vector< std::pair<double, double> > points2 = { { 0,0 },{ 1,-1 },{ 2,-16 },{ 3,-81 },{ 4,-256 } };
    printCoe(solveEq(points2));

    printCoe(solveEq({ { 0,0 },{ 1,1 },{ 2,8 },{ 3,27 } }));

    std::cin.ignore();
    return 0;
}