98 lines
3.4 KiB
C++
98 lines
3.4 KiB
C++
//
|
|
// Created by quentin on 9/13/22.
|
|
//
|
|
|
|
#include "math.h"
|
|
|
|
//Polynomial Fit
|
|
#include <iostream>
|
|
#include <cmath>
|
|
|
|
using namespace Budget::Utilities;
|
|
|
|
PolynomialFunction Math::polynomialFit(int degree, int numberOfPoints, const double x[], const double y[]) {
|
|
//Array that will store the values of sigma(xi),sigma(xi^2),sigma(xi^3)....sigma(xi^2numberOfPoints)
|
|
double X[2 * degree + 1];
|
|
for (int i = 0; i < 2 * degree + 1; i++) {
|
|
X[i] = 0;
|
|
for (int j = 0; j < numberOfPoints; j++) {
|
|
//consecutive positions of the array will store numberOfPoints,sigma(xi),sigma(xi^2),sigma(xi^3)....sigma(xi^2n)
|
|
X[i] = X[i] + pow(x[j], i);
|
|
}
|
|
}
|
|
//B is the Normal matrix(augmented) that will store the equations, 'a' is for value of the final coefficients
|
|
double B[degree + 1][degree + 2], a[degree + 1];
|
|
for (int i = 0; i <= degree; i++) {
|
|
for (int j = 0; j <= degree; j++) {
|
|
//Build the Normal matrix by storing the corresponding coefficients at the right positions except the last column of the matrix
|
|
B[i][j] = X[i + j];
|
|
}
|
|
}
|
|
|
|
//Array to store the values of sigma(yi),sigma(xi*yi),sigma(xi^2*yi)...sigma(xi^degree*yi)
|
|
double Y[degree + 1];
|
|
for (int i = 0; i < degree + 1; i++) {
|
|
Y[i] = 0;
|
|
for (int j = 0; j < numberOfPoints; j++) {
|
|
//consecutive positions will store sigma(yi),sigma(xi*yi),sigma(xi^2*yi)...sigma(xi^degree*yi)
|
|
Y[i] = Y[i] + pow(x[j], i) * y[j];
|
|
}
|
|
}
|
|
for (int i = 0; i <= degree; i++) {
|
|
//load the values of Y as the last column of B(Normal Matrix but augmented)
|
|
B[i][degree + 1] = Y[i];
|
|
}
|
|
|
|
//degree is made degree+1 because the Gaussian Elimination part below was for degree equations, but here degree is the degree of polynomial and for degree degree we get degree+1 equations
|
|
degree = degree + 1;
|
|
//From now Gaussian Elimination starts(can be ignored) to solve the set of linear equations (Pivotisation)
|
|
for (int i = 0; i < degree; i++) {
|
|
for (int k = i + 1; k < degree; k++) {
|
|
if (B[i][i] < B[k][i]) {
|
|
for (int j = 0; j <= degree; j++) {
|
|
double temp = B[i][j];
|
|
B[i][j] = B[k][j];
|
|
B[k][j] = temp;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//loop to perform the gauss elimination
|
|
for (int i = 0; i < degree - 1; i++) {
|
|
for (int k = i + 1; k < degree; k++) {
|
|
double t = B[k][i] / B[i][i];
|
|
for (int j = 0; j <= degree; j++) {
|
|
//make the elements below the pivot elements equal to zero or elimnate the variables
|
|
B[k][j] = B[k][j] - t * B[i][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
//back-substitution
|
|
for (int i = degree - 1; i >= 0; i--) {
|
|
//x is an array whose values correspond to the values of x,y,z..
|
|
//make the variable to be calculated equal to the rhs of the last equation
|
|
a[i] = B[i][degree];
|
|
for (int j = 0; j < degree; j++) {
|
|
//then subtract all the lhs values except the coefficient of the variable whose value is being calculated
|
|
if (j != i) {
|
|
a[i] = a[i] - B[i][j] * a[j];
|
|
}
|
|
}
|
|
//now finally divide the rhs by the coefficient of the variable to be calculated
|
|
a[i] = a[i] / B[i][i];
|
|
}
|
|
std::cout << "\nThe values of the coefficients are as follows:\n";
|
|
for (int i = 0; i < degree; i++)
|
|
std::cout << "x^" << i << "=" << a[i] << std::endl; // Print the values of x^0,x^1,x^2,x^3,....
|
|
std::cout << "\nHence the fitted Polynomial is given by:\ny=";
|
|
for (int i = 0; i < degree; i++)
|
|
std::cout << " + (" << a[i] << ")" << "x^" << i;
|
|
std::cout << "\n";
|
|
|
|
std::vector<double> va(a, a + degree);
|
|
return PolynomialFunction(va);
|
|
}
|
|
|