Как найти коэффициенты полиномиального уравнения?
Учитывая два пункта в x, y самолет:
x, f(x)
1, 3
2, 5
Я могу интерполировать их с помощью Лагранжа и найти 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;
}
Работая усердно, я действительно смог найти коэффициенты для уравнений порядка до 4-го числа.
Как найти коэффициенты заказа n уравнения? куда
Pn(x) = c_2x^2 + c_1x^1 + c_0x^0 + ...
Решение
Это простая задача линейной алгебры.
У нас есть набор из N образцов вида xК -> f (xК) и мы знаем общий вид функции f (x), которая:
f (x) = c0Икс0 + с1Икс1 + … + сN-1ИксN-1
Мы хотим найти коэффициенты с0 … сN-1. Для этого мы строим систему из N уравнений вида:
с0ИксК0 + с1ИксК1 + … + сN-1ИксКN-1 = f (xК)
где К это номер образца. С хК и f (xК) являются константами, а не переменными, мы имеем линейную систему уравнений.
Выражаясь в терминах линейной алгебры, мы должны решить:
Ac = b
где А является Матрица Вандермонда полномочий Икс а также б является вектором f (xК) ценности.
Чтобы решить такую систему, вам нужна библиотека линейной алгебры, такая как Eigen, Увидеть Вот например код.
Единственное, что может пойти не так при таком подходе, — это то, что система линейных уравнений недоопределена, что произойдет, если ваши N выборок могут быть сопоставлены с многочленом степени меньше, чем N-1. В таком случае вы все еще можете решить эту систему с помощью псевдообращения Мура-Пенроуза следующим образом:
c = pinv (A) * b
К несчастью, Eigen не имеет pinv() реализация, хотя это довольно легко кодировать самостоятельно с точки зрения разложения по сингулярным значениям (SVD).
Другие решения
Я создал наивную реализацию матричного решения:
#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;
}