Interpolacja

// Aproksymacja sklejana
// (C)2026 mgr Jerzy Wałaszek
// Metody numeryczne 0077
//---------------------------

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

using namespace std;

// Definicje stałych
//------------------

// Przedział interpolacji
// od XS do XE
const double XS = -M_PI / 2;
const double XE =  M_PI / 2;
// Ilość węzłów interpolacji
const int NN = 10;
// Ilość punktów interpolowanych
const int  N = 10;
// Szerokość segmentu
const double DX = (XE-XS)/(NN-1);
// Przyrównanie do zera
const double EPS = 1e-12;
// Tablice
//--------

// Współrzędne węzłów
double nx[NN];
double ny[NN];
// Współczynniki wielomianów
// w kolejnych segmentach
double A[NN-1][4];

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

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

// Zwraca wartość pseudolosową
// w przedziale [0;1)
//----------------------------
double random()
{
  return rand() /
         ((double)RAND_MAX + 1);
}

// Wyznacza węzły interpolacji
//----------------------------
void set_nodes(void)
{
  int i;

  // wyliczamy NN punktów x,y
  for(i = 0; i < NN; i++)
  {
    nx[i] = XS + i * DX;
    ny[i] = f(nx[i]);
  }
}

// Funkcja wylicza współczynniki
// dla wielomianów kubicznych.
// Współczynniki są w tablicy
// A[NN-1][4]. Współczynniki
// są czwórkami: a b c d dla
// kolejnych wielomianów.
// Zwraca true, gdy operacja
// się powiedzie.
bool set_a(void)
{
  // Najpierw ustawiamy macierz
  // rozszerzoną układu równań
  int n = 4 * (NN - 1); // wiersze
  double AB[n][n+1];
  int i,j;
  // Zerowanie macierzy AB
  for(i = 0; i < n; i++)
    for(j = 0; j <= n; j++)
      AB[i][j] = 0;

  // Wiersze z punktu 1 (wartości)
  int row = 0;
  for(i = 0; i < NN - 1; i++)
  {
    // Początek segmentu i
    AB[row][i*4] = 1; // a_i
    AB[row][n] = ny[i];
    row++;

    // Koniec segmentu i
    AB[row][i*4] = 1;          // a_i
    AB[row][i*4+1] = DX;       // b_i
    AB[row][i*4+2] = DX*DX;    // c_i
    AB[row][i*4+3] = DX*DX*DX; // d_i
    AB[row][n] = ny[i+1];
    row++;
  }

  // Wiersze z punktu 2 (pochodne S')
  for(i = 0; i < NN - 2; i++)
  {
    AB[row][i*4+1] = 1;      //  b_i
    AB[row][i*4+2] = 2*DX;   //  c_i
    AB[row][i*4+3] = 3*DX*DX;//  d_i
    AB[row][(i+1)*4+1] = -1; // -b_{i+1}
    row++;
  }

  // Wiersze z punktu 3 (pochodne S'')
  for(i = 0; i < NN - 2; i++)
  {
    AB[row][i*4+2] = 2;      // 2c_i
    AB[row][i*4+3] = 6*DX;   // 6*DX*d_i
    AB[row][(i+1)*4+2] = -2; // -2c_{i+1}
    row++;
  }

  // Wiersze z punktu 4 (brzegi)
  AB[row][2] = 2; // 2c_0 = 0
  row++;

  int last = NN - 2;
  AB[row][last*4+2] = 2;    // 2c_{last}
  AB[row][last*4+3] = 6*DX; // 6*DX*d_{last}
  row++;

  // Eliminacja Gaussa
  int k;
  for(i = 0; i < n; i++)
  {
    // Szukanie elementu osiowego
    for(j = i + 1; j < n; j++)
    {
      if(fabs(AB[j][i]) >
         fabs(AB[i][i]))
      {
        for(k = i; k <= n; k++)
        {
          double tmp = AB[i][k];
          AB[i][k] = AB[j][k];
          AB[j][k] = tmp;
        }
      }
    }

    // Sprawdzenie osobliwosci
    if(fabs(AB[i][i]) < EPS)
      return false;

    // Zerowanie pod przekątną
    for(j = i + 1; j < n; j++)
    {
      double factor =
        AB[j][i] / AB[i][i];
      for(k = i; k <= n; k++)
        AB[j][k] -=
          factor * AB[i][k];
    }
  }

  // Postępowanie wsteczne
  double X[n];
  for(i = n - 1; i >= 0; i--)
  {
    double sum = 0;
    for(j = i + 1; j < n; j++)
      sum += AB[i][j] * X[j];
    X[i] = (AB[i][n] - sum) /
           AB[i][i];
  }

  // Przepisanie do tablicy A
  for(i = 0; i < NN - 1; i++)
  {
    A[i][0] = X[i*4];   // a
    A[i][1] = X[i*4+1]; // b
    A[i][2] = X[i*4+2]; // c
    A[i][3] = X[i*4+3]; // d
  }

  return true;
}

// Funkcja zwraca numer
// segmentu dla x
int find_s(double x)
{
  return floor((x - XS) / DX);
}

// Wylicza wartość S_i(x)
//-----------------------
double f_i(double x)
{
  int i = find_s(x);
  double dx = x - nx[i];

  // Schemat Hornera
  return A[i][0] + dx * (
         A[i][1] + dx * (
         A[i][2] + dx *
         A[i][3]));
}

// Program główny

int main()
{
  int i;

  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(4)
       << fixed;

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

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

  // Obliczamy współczynniki
  if(set_a())
  {
    cout << "Interpolacja sklejana\n"
            "---------------------\n\n";

    // Generujemy N punktów x,
    // liczymy dla nich y
    // 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 obliczeniach\n";

  cout << endl;
  system("pause");
  return 0;
}