212 lines
6.9 KiB
C++
212 lines
6.9 KiB
C++
|
#include <algorithm>
|
||
|
|
#include <cmath>
|
||
|
|
#include <cstdlib>
|
||
|
|
#include <cstring>
|
||
|
|
#include <fstream>
|
||
|
|
#include <iomanip>
|
||
|
|
#include <ios>
|
||
|
|
#include <iostream>
|
||
|
|
#include <numeric>
|
||
|
|
#include <string>
|
||
|
|
#include <utility>
|
||
|
|
#include <vector>
|
||
|
|
|
||
|
|
using namespace std::literals;
|
||
|
|
|
||
|
|
using points_t = std::pair<std::vector<double>, std::vector<double>>;
|
||
|
|
using matrix_t = std::vector<std::vector<double>>;
|
||
|
|
|
||
|
|
std::vector<double> newton_differences(points_t points) {
|
||
|
|
std::vector<double> factors{};
|
||
|
|
auto &[x, fx] = points;
|
||
|
|
int n = points.first.size();
|
||
|
|
for (int i = 1; i <= n; ++i) {
|
||
|
|
factors.push_back(fx[0]);
|
||
|
|
for (int j = 0; j < (n - i); ++j) {
|
||
|
|
fx[j] = (fx[j + 1] - fx[j]) / (x[j + i] - x[j]);
|
||
|
|
}
|
||
|
|
|
||
|
|
}
|
||
|
|
return factors;
|
||
|
|
}
|
||
|
|
|
||
|
|
double newton_polynomial(points_t points, std::vector<double> factors, double x) {
|
||
|
|
auto xs = points.first;
|
||
|
|
int n = points.second.size() - 1;
|
||
|
|
double val = 0;
|
||
|
|
for (int k = 0; k <= n; ++k) {
|
||
|
|
double acc = 1;
|
||
|
|
for(int i = 0; i < k; ++i) {
|
||
|
|
acc *= (x - xs[i]);
|
||
|
|
}
|
||
|
|
val += acc * factors[k];
|
||
|
|
}
|
||
|
|
return val;
|
||
|
|
}
|
||
|
|
|
||
|
|
void calculate_natural_cubic_spline_matrix(points_t points, matrix_t &mat) {
|
||
|
|
auto &[xs, fx] = points;
|
||
|
|
int n = xs.size();
|
||
|
|
|
||
|
|
// Construção da matriz
|
||
|
|
for (int i = 1; i < n - 1; ++i) {
|
||
|
|
mat[i][i - 1] = (xs[i] - xs[i - 1])/6;
|
||
|
|
mat[i][i] = (xs[i + 1] - xs[i - 1])/3;
|
||
|
|
mat[i][i + 1] = (xs[i + 1] - xs[i])/6;
|
||
|
|
mat[i][n] = (fx[i + 1] - fx[i])/(xs[i + 1] - xs[i]) - (fx[i] - fx[i - 1])/(xs[i] - xs[i - 1]);
|
||
|
|
}
|
||
|
|
mat[0][0] = 1;
|
||
|
|
mat[0][n - 1] = 0;
|
||
|
|
mat[n - 1][n - 1] = 1;
|
||
|
|
mat[n - 1][n] = 0;
|
||
|
|
|
||
|
|
std::ofstream matriz("matriz.txt", std::ios_base::app);
|
||
|
|
|
||
|
|
matriz << "mat:\n";
|
||
|
|
for (auto i : mat){
|
||
|
|
for (auto j = i.begin(); j != i.end() - 1; ++ j)
|
||
|
|
matriz << *j << ",";
|
||
|
|
matriz << "\n";
|
||
|
|
}
|
||
|
|
|
||
|
|
matriz << "b:\n";
|
||
|
|
for (auto i : mat){
|
||
|
|
matriz << i[n] << '\n';
|
||
|
|
}
|
||
|
|
|
||
|
|
// Passar para a forma triangular
|
||
|
|
for (int k = 0; k < n; ++k) {
|
||
|
|
for (int i = k + 1; i < n; ++i) {
|
||
|
|
if (mat[k][k] != 0) {
|
||
|
|
double mul = mat[i][k]/mat[k][k];
|
||
|
|
for (int j = k; j < n; ++j) {
|
||
|
|
mat[i][j] -= mul * mat[k][j];
|
||
|
|
}
|
||
|
|
mat[i][n] -= mul * mat[k][n];
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// Resolução da matriz
|
||
|
|
for (int i = n - 1; i > 0; --i) {
|
||
|
|
if (mat[i][i] != 0) {
|
||
|
|
double mul = mat[i - 1][i]/mat[i][i];
|
||
|
|
for (int j = 1; j < n + 1; ++j) {
|
||
|
|
mat[i - 1][j] -= mul * mat[i][j];
|
||
|
|
}
|
||
|
|
mat[i - 1][i] = 0;
|
||
|
|
mat[i][n] /= mat[i][i];
|
||
|
|
mat[i][i] = 1;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
matriz << "x:\n";
|
||
|
|
for (auto i : mat){
|
||
|
|
matriz << i[n] << '\n';
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
double natural_cubic_spline(points_t points, matrix_t &mat, double x) {
|
||
|
|
|
||
|
|
auto &[xs, fx] = points;
|
||
|
|
int n = xs.size();
|
||
|
|
|
||
|
|
int i = 0;
|
||
|
|
for (int i_ = 0; i_ < n; ++i_) {
|
||
|
|
if (xs[i_] > x){
|
||
|
|
i = i_;
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
double hi = xs[i] - xs[i - 1];
|
||
|
|
return mat[i - 1][n] * std::pow((xs[i] - x), 3)/(6 * hi) +
|
||
|
|
mat[i][n] * std::pow((x - xs[i - 1]), 3)/(6 * hi) +
|
||
|
|
(fx[i - 1] - mat[i - 1][n] * (hi * hi)/6)*(xs[i] - x)/hi +
|
||
|
|
(fx[i] - mat[i][n] * (hi * hi)/6)*(x - xs[i - 1])/hi;
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
void exercise_a(points_t points) {
|
||
|
|
auto &[xs, fx] = points;
|
||
|
|
auto factors = newton_differences(points);
|
||
|
|
std::ofstream poly{"a_polynomial.txt"};
|
||
|
|
for(double x = 0; x < 4; x += 0.001)
|
||
|
|
poly << x << " " << newton_polynomial(points, factors, x) << '\n';
|
||
|
|
|
||
|
|
std::ofstream spline{"a_spline.txt"};
|
||
|
|
unsigned long n = xs.size();
|
||
|
|
matrix_t mat(n, std::vector<double>(n + 1));
|
||
|
|
calculate_natural_cubic_spline_matrix(points, mat);
|
||
|
|
for(double x = 0; x < 4; x += 0.001)
|
||
|
|
spline << x << " " << natural_cubic_spline(points, mat, x) << '\n';
|
||
|
|
}
|
||
|
|
|
||
|
|
void exercise_b() {
|
||
|
|
points_t points;
|
||
|
|
auto &[xs, fx] = points;
|
||
|
|
auto f = [](double x) { return 4 * std::pow(x, 2) + std::sin(9 * x); };
|
||
|
|
|
||
|
|
for (double x = -1; x <= 1; x += (1 - -1)/8.0) {
|
||
|
|
xs.push_back(x);
|
||
|
|
fx.push_back(f(x));
|
||
|
|
}
|
||
|
|
|
||
|
|
std::ofstream poly{"b_polynomial.txt"};
|
||
|
|
auto factors = newton_differences(points);
|
||
|
|
for(double x = -1; x < 1; x += 0.001)
|
||
|
|
poly << x << " " << newton_polynomial(points, factors, x) << '\n';
|
||
|
|
|
||
|
|
std::ofstream spline{"b_spline.txt"};
|
||
|
|
unsigned long n = xs.size();
|
||
|
|
matrix_t mat(n, std::vector<double>(n + 1));
|
||
|
|
calculate_natural_cubic_spline_matrix(points, mat);
|
||
|
|
for(double x = -1; x < 1; x += 0.001)
|
||
|
|
spline << x << " " << natural_cubic_spline(points, mat, x) << '\n';
|
||
|
|
|
||
|
|
std::cout << "Erros poly: "
|
||
|
|
<< std::abs(std::round((f(0.3) - newton_polynomial(points, factors, 0.3)) * 100.0) / 100.0) << " "
|
||
|
|
<< std::abs(std::round((f(0.83) - newton_polynomial(points, factors, 0.83)) * 100.0) / 100.0) << '\n';
|
||
|
|
|
||
|
|
std::cout << "Erros spline: "
|
||
|
|
<< std::abs(std::round((f(0.3) - natural_cubic_spline(points, mat, 0.3))*100.0)/100.0) << " "
|
||
|
|
<< std::abs(std::round((f(0.83) - natural_cubic_spline(points, mat,0.83))*100.0)/100.0) << '\n';
|
||
|
|
|
||
|
|
std::cout << "p(0.3) = " << newton_polynomial(points, factors, 0.3) << ' '
|
||
|
|
<< "p(0.83) = " << newton_polynomial(points, factors, 0.83) << '\n';
|
||
|
|
|
||
|
|
std::cout << "p(0.3) = " << natural_cubic_spline(points, mat, 0.3) << ' '
|
||
|
|
<< "p(0.83) = " << natural_cubic_spline(points, mat, 0.83) << '\n';
|
||
|
|
}
|
||
|
|
|
||
|
|
int main(int argc, char **argv) {
|
||
|
|
std::istream *s;
|
||
|
|
std::ifstream file;
|
||
|
|
if (argc == 2) {
|
||
|
|
if (argv[1] == "-"s) {
|
||
|
|
s = &std::cin;
|
||
|
|
} else {
|
||
|
|
file.open(argv[1]);
|
||
|
|
s = &file;
|
||
|
|
}
|
||
|
|
} else{
|
||
|
|
std::cerr << "An argument required\n";
|
||
|
|
return 1;
|
||
|
|
}
|
||
|
|
|
||
|
|
points_t points;
|
||
|
|
std::pair<double, double> p;
|
||
|
|
while(true) {
|
||
|
|
*s >> p.first >> p.second;
|
||
|
|
if(s->eof())
|
||
|
|
break;
|
||
|
|
|
||
|
|
points.first.push_back(p.first);
|
||
|
|
points.second.push_back(p.second);
|
||
|
|
}
|
||
|
|
|
||
|
|
exercise_a(points);
|
||
|
|
exercise_b();
|
||
|
|
|
||
|
|
}
|