Interpolation

#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>


// Нужно для использования метода прогонки при интерполяции сплайнами
#include "D:\#Cpp_Projects\Interpolation\SlaeSolving\SlaeSolver.h"
#include "D:\#Cpp_Projects\Interpolation\SlaeSolving\SquareMatrix.cpp"
#include "D:\#Cpp_Projects\Interpolation\SlaeSolving\TridiagonalMatrix.cpp"
#include "D:\#Cpp_Projects\Interpolation\SlaeSolving\ExtendedMatrix.cpp"
#include "D:\#Cpp_Projects\Interpolation\SlaeSolving\SlaeSolver.cpp"
#include "D:\#Cpp_Projects\Interpolation\SlaeSolving\Tester.cpp"


using namespace std;


double sqr(double x) { return x * x; };
double cb(double x) { return x * x * x; };


// Список доступных в программе методов интерполяции
enum class Method { Lagranz, Splines };



// Класс, описывающий функцию y = f(x), заданную в виде таблицы
class Table
{
private:
	vector<double> x;
	vector<double> y;
	size_t n = 0;

public:
	Table(vector<double> x, vector<double> y) : x(x), y(y)
	{
		if (x.size() != y.size())
			throw invalid_argument("Incorrect function table initialization!");
		this->n = x.size();
	}

	vector<double> getX() { return x; }
	vector<double> getY() { return y; }
	size_t getN() { return n; }

	bool argsAreOrederedAscending() // Упорядочен ли столбец иксов? (по возрастанию)
	{
		for (vector<double>::iterator i = x.begin(); i != x.end() - 1; ++i)
			if (!(*(i + 1) >= *i))
				return false;
		return true;
	}
};



// Класс, содержащий методы для интерполяции таблично заданных функций различными способами
class Interpolator
{
public:
	static double lagranzMethod(Table f, double x) // Методом Лагранжа
	{
		if (f.getN() == 0)
			throw invalid_argument("Table must not be empty!");

		const int n = f.getN();
		vector<double> X = f.getX();
		vector<double> Y = f.getY();

		if(x > *std::max_element(X.begin(), X.end()) || x < *std::min_element(X.begin(), X.end()))
			throw invalid_argument("X must be inside interpolation section!");

		double p = 1.0; // Произведение элементов главной диагонали таблицы разностей
		for (int i = 0; i < n; ++i)
			p *= (x - X[i]);

		double sum = 0.0; // Сумма частных (y / k), где y - i-ое значение функции, а k - произведение элементов i-ой строки таблицы разностей
		for (int i = 0; i < n; ++i)
		{
			double k = 1.0; // Произведение элементов i-ой строки таблицы разностей
			for (int j = 0; j < n; ++j)
				if (i != j)
					k *= (X[i] - X[j]);
			k *= x - X[i]; // Домножаем k на элемент главной диагонали таблицы разностей

			sum += Y[i] / k; // Находим сумму частных
		}

		return p * sum; // Ответ - произведение суммы частных и произведения эл-тов главной диагонали таблицы разностей
	}

	static double splinesMethod(Table f, double x) // Через сплайны
	{
		if (f.getN() == 0)
			throw invalid_argument("Table must not be empty!");

		if(!f.argsAreOrederedAscending())
			throw invalid_argument("Arguments row in table must be ordered ascending!");


		// Получаем данные о табличной функции
		const int n = f.getN();
		const vector<double> X = f.getX();
		const vector<double> Y = f.getY();


		// Вычисляем h i-ые шаги между соседними точками в таблице
		vector<double> h(n - 1);
		for (int i = 0; i != n - 1; ++i)
			h[i] = X[i+1] - X[i];

		vector<double> s(n); // Наклоны сплайна

		// Матрица и столбец свободных членов для трёхдиагональной слау
		TridiagonalMatrix matrix(n);
		vector<double> column(n);

		// Первая строка для слау
		matrix.at(0, 0) = -4.0 / h[0];
		matrix.at(0, 1) = -2.0 / h[0];
		column.at(0) = -6.0 * (Y[1] - Y[0]) / (h[0] * h[0]);

		// Последняя строка для слау
		matrix.at(n - 1, n - 2) = 2.0 / h[n - 2];
		matrix.at(n - 1, n - 1) = 4.0 / h[n - 2];
		column.at(n - 1) = 6.0 * (Y[n - 1] - Y[n - 2]) / (h[n - 2] * h[n - 2]);

		// i-ая строка слау, 1 <= i <= n-2
		for (int i = 1; i != n - 1; ++i)
		{
			matrix.at(i, i - 1) = 1.0 / h[i - 1];
			matrix.at(i, i) = 2 * (1.0 / h[i - 1] + 1.0 / h[i]);
			matrix.at(i, i + 1) = 1.0 / h[i];
			column.at(i) = 3 * ((Y[i] - Y[i - 1]) / (h[i - 1] * h[i - 1]) + (Y[i + 1] - Y[i]) / (h[i] * h[i]));
		}

		// Передаём нашу слау в класс, который её решит
		SlaeSolver solver(ExtendedMatrix(matrix.getMatrix(), column));

		// Вычисляем наклоны сплайна методом прогонки
		s = solver.sweepMethod();

		// Вычисляем многочлен Эрмита для заданного промежутка
		for (int i = 1; i < n; ++i)
			if (x >= X[i - 1] && x <= X[i])
				return Y[i - 1] * sqr(x - X[i]) * (2 * (x - X[i - 1]) + h[i - 1]) / cb(h[i - 1]) +
				s[i - 1] * sqr(x - X[i]) * (x - X[i - 1]) / sqr(h[i - 1]) +
				Y[i] * sqr(x - X[i - 1]) * (2 * (X[i] - x) + h[i - 1]) / cb(h[i - 1]) +
				s[i] * sqr(x - X[i - 1]) * (x - X[i]) / sqr(h[i - 1]);

		throw invalid_argument("X must be inside interpolation section!");
	}
};



// Класс, описывающий интерполяционный многочлен, полученный при интерполяции функции тем или иным способом
class Polynomial
{
private:
	Table f;
	Method method;

public:
	Polynomial(Table f, Method method) : f(f)
	{
		if (f.getN() == 0)
			throw invalid_argument("Table must not be empty!");

		// Для метода сплайнов нужно проверить, что иксы упорядочены по возрастанию
		if(method == Method::Splines && !f.argsAreOrederedAscending())
			throw invalid_argument("Arguments row in table must be ordered ascending!");

		this->method = method;
	}

	// Считаем значения P(x) в заданной точке
	double at(double x)
	{
		try
		{
			switch (method)
			{
			case Method::Lagranz: // P(x) сформирован методом Лагранжа
				return Interpolator::lagranzMethod(this->f, x);
			case Method::Splines: // P(x) сформирован через сплайны
				return Interpolator::splinesMethod(this->f, x);
			default:
				throw invalid_argument("You must specify method of interpolation correctly!");
			}
		}
		catch (invalid_argument & e)
		{
			cerr << e.what() << endl;
		}
	}
};



int main()
{
	vector<double> x = { 0.5, 1.0, 7.0, 11.0, 12.0, 15.0, 20.0}; // Столбец иксов для таблицы
	vector<double> y = { 23.4, 2.6, -16.0, -186.5, 58.9, 132.0, 300.0}; // Столбец игреков для таблицы

	Table f(x, y); // Таблично заданная функция

	// Многочлен, полученный путём интерполяции функции методом Лагранжа
	Polynomial p(f, Method::Lagranz);

	// Вывод значений многочлена
	cout << " ================== Lagranz method ================== " << endl;
	const int width = 15;
	cout << setw(width) << "P (0) = " << setw(width) << p.at(0.0) << endl;
	cout << setw(width) << "P (2) = " << setw(width) << p.at(2.0) << endl;
	cout << setw(width) << "P (11.5) = " << setw(width) << p.at(11.5) << endl;
	cout << setw(width) << "P (16) = " << setw(width) << p.at(16.0) << endl;
	cout << setw(width) << "P (35) = " << setw(width) << p.at(35.0) << endl;
	cout << endl << endl;



	// Многочлен, полученный путём интерполяции функции через сплайны
	Polynomial p1(f, Method::Splines);

	// Вывод значений многочлена
	cout << " ================== Splines method ================== " << endl;
	cout << setw(width) << "P (0) = " << setw(width) << p1.at(0.0) << endl;
	cout << setw(width) << "P (2) = " << setw(width) << p1.at(2.0) << endl;
	cout << setw(width) << "P (11.5) = " << setw(width) << p1.at(11.5) << endl;
	cout << setw(width) << "P (16) = " << setw(width) << p1.at(16.0) << endl;
	cout << setw(width) << "P (35) = " << setw(width) << p1.at(35.0) << endl;



	getchar();
	return 0;
}