Interpolacja

// Interpolacja wielomianowa
// (C)2026 mgr Jerzy Wałaszek
// Metody numeryczne 0074
//---------------------------

#include <iostream>
#include <windows.h>
#include <iomanip>
#define _USE_MATH_DEFINES
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;

// STAŁE
//------

// Przedział interpolacji
// od XS do XE
const double XS = -M_PI / 2;
const double XE =  M_PI / 2;
// liczba węzłów
const int NN = 6;
// Liczba interpolowanych
// punktów
const int N = 10;
// Dokładność przyrównania do 0
const double EPS = 0.00001;

// Tablice
//--------

// Współrzędne x węzłów
double nx[NN];
// Współrzędne y węzłów
double ny[NN];
// Współczynniki wielomianu
double a[NN];

// Funkcje
//--------

// Funkcja interpolowana
//----------------------
double f(double x)
{
  return sin(x);
}

// Funkcja zwraca wartość
// pseudolosową w przedziale
// <0,1)
//--------------------------
double random()
{
  return rand() /
         ((double)RAND_MAX + 1);
}

// Funkcja wyznacza węzły
// dla wielomianu
//-----------------------
void set_nodes()
{
  // Dwa pierwsze węzły
  // definiują przedział
  // interpolacji
  nx[0] = XS;
  nx[1] = XE;

  // Pozostałe punkty są
  // rozsiane po przedziale
  int i = 2,j;
  double xx,dx;

  // minimalna odległość w osi OX
  // pomiędzy węzłami
  dx = (XE - XS) / 2 / NN;
  while(i < NN)
  {
    xx = XS + random() * (XE - XS);
    bool t = false;
    for(j = 0; j < i; j++)
      // Sprawdzamy odległość
      // węzłów w osi OX
      if(fabs(xx - nx[j]) < dx)
      {
        t = true;
        break;
      }
    // Węzeł za blisko innego,
    // losujemy go ponownie
    if(t) continue;
    nx[i++] = xx;
  }

  // Dla współrzędnych x
  // liczymy współrzędne y
  for(i = 0; i < NN; i++)
    ny[i] = f(nx[i]);
}

// Oblicza wartość wielomianu
//---------------------------
double f_i(double x)
{
  // Współczynnik przy x^(NN-1)
  double w = a[NN - 1];
  int i;
  for(i = NN-1; i > 0; i--)
  {
    w *= x;
    w += a[i-1];
  }
  return w;
}

// Funkcja rozwiązuje układ równań
// metodą eliminacji Gaussa
// NN - stopień macierzy
// AB - macierz rozszerzona NN x (NN + 1)
// X  - macierz rozwiązań NN x 1
// Jeśli układ został rozwiązany,
// to wynikiem jest true.
// Inaczej wynikiem jest false.
//--------------------------------
bool GaussCrout(double AB[][NN+1],
                double X[NN])
{
  int i,j,k;
  double m,s;
  int WK[NN + 1];

  // Przygotowujemy wektor
  // kolumnowy
  for(i = 0; i <= NN; i++)
    WK[i] = i;

  // Eliminacja
  for(i = 0; i < NN - 1; i++)
  {
    for(k = i, j = i + 1; j < NN; j++)
      if(fabs(AB[i][WK[k]]) <
         fabs(AB[i][WK[j]]))
        k = j;

    if(k != i)
      swap(WK[k],WK[i]);

    if(fabs(AB[i][WK[i]]) < EPS)
    {
      return false;
    }

    for(j = i + 1; j < NN; j++)
    {
      m = -AB[j][WK[i]] / AB[i][WK[i]];
      for(k = i + 1; k <= NN; k++)
        AB[j][WK[k]] += m * AB[i][WK[k]];
    }
  }

  // Obliczanie wartości niewiadomych
  for(i = NN - 1; i >= 0; i--)
  {
    if(fabs(AB[i][WK[i]]) < EPS)
    {
      return false;
    }
    s = AB[i][NN];
    for(j = NN - 1; j > i; j--)
      s -= AB[i][WK[j]] * X[WK[j]];
    X[WK[i]] = s / AB[i][WK[i]];
  }

  return true;
}

// Program główny
//---------------
int main()
{
  int i,j;

  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(4)
       << fixed;

  // Inicjujemy generator
  // pseudolosowy
  srand(time(nullptr));

  // Generujemy węzły
  set_nodes();

  // Tworzymy macierz
  // układu równań
  double ab[NN][NN+1];

  for(i = 0; i < NN; i++)
  {
    // W ostatniej kolumnie
    // są współrzędne y
    ab[i][NN] = ny[i];
    // W pierwszej są 1
    ab[i][0] = 1;
    // Pozostałe są potęgami
    for(j = 1; j < NN; j++)
      ab[i][j] = nx[i] * ab[i][j-1];
  }

  // Wyznaczamy współczynniki
  // wielomianu
  if(GaussCrout(ab,a))
  {
    cout <<
    "Interpolacja wielomianowa\n"
    "-------------------------\n\n";

    // Generujemy N punktów x
    // i liczymy dla nich y
    // przy pomocy wielomianu
    // i wyświetlamy wynik
    double x,y;
    for(i = 0; i < N; i++)
    {
      x = XS + random() * (XE - XS);
      y = f_i(x);
      cout << "x = "
           << setw(7) << x
           << ", f(x) = "
           << setw(7) << f(x)
           << ", interpolacja f(x) = "
           << setw(7) << y
           << endl;
    }
  }
  else
    cout <<
    "Błąd w trakcie obliczeń!!!\n";
  cout << endl;
  system("pause");
  return 0;
}