1

x, y给定平面上的两点：

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

我可以使用 Lagrange 和 find 对它们进行插值f(1.5)，从而得到4. 想了想，我设法找到了一种方法来发现方程的系数：

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;
}

这给了我 P1(x) = 2x +1。

再想一想，我能够将其扩展到2nd排序方程。因此，鉴于以下几点：

1, 1
2, 4
3, 9

我找到了以下等式P2(x) = 1x^2 +0x +0：

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;
}

努力工作，我实际上能够找到四阶方程的系数。

如何求阶n方程的系数？在哪里

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

2 回答 2

1

这是一个简单的线性代数问题。

我们有一组形式为 x k -> f(x k ) 的 N 个样本，我们知道函数 f(x) 的一般形式，即：

f(x) = c 0 x 0 + c 1 x 1 + ... + c N-1 x N-1

我们想要找到系数 c 0 ... c N-1。为了实现这一点，我们构建了一个由 N 个方程组成的系统，其形式为：

c 0 x k 0 + c 1 x k 1 + ... + c N-1 x k N-1 = f(x k )

其中k是样本数。因为 x k和 f(x k ) 是常数而不是变量，所以我们有一个线性方程组。

用线性代数表示，我们必须解决：

Ac = b

其中 A 是x的幂的Vandermonde 矩阵，b是 f(x k ) 值的向量。

要解决这样的系统，您需要一个线性代数库，例如Eigen. 有关示例代码，请参见此处。

这种方法唯一可能出错的是线性方程组未确定，如果您的 N 个样本可以与次数小于 N-1 的多项式拟合，就会发生这种情况。在这种情况下，您仍然可以使用 Moore-Penrose 伪逆求解这个系统，如下所示：

c = pinv(A)*b

不幸的是，Eigen没有pinv()实现，尽管根据奇异值分解 (SVD) 自己编写代码非常容易。

于 2017-04-23T11:02:45.650 回答
0

我创建了一个简单的矩阵解决方案实现：

#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;
}
于 2017-04-23T19:54:39.873 回答