Interpolacja

Podstawową wadą interpolacji wielomianowej jest zjawisko Rungego, które powoduje duże oscylacje wielomianu interpolacyjnego o wysokim stopniu na krańcach przedziału interpolacji. Metoda interpolacji Chebyszewa eliminuje tę wadę. Osiąga to poprzez nierównomierne rozmieszczenie węzłów w przedziale interpolacji - węzły są rozmieszczane gęściej na krańcach przedziału, co stabilizuje wielomian interpolacyjny i zmniejsza błędy interpolacji.

Na rozwiązanie problemu interpolacji wpadł genialny matematyk rosyjski Pafnutij Chebyszew (Пафнутий Львович Чебышёв – ur. 4 maja?/16 maja 1821 w Okatowie w guberni  kałuskiej, zm. 26 listopada?/8 grudnia 1894 w Sankt Petersburgu), który wymyślił prosty sposób: zamiast węzłów równoodległych, należy zagęścić je na krańcach przedziału interpolacyjnego.

Algorytm interpolacji Chebyszewa składa się z kilku kroków:

Chcemy interpolować wartości pewnej funkcji \({y = f(x)}\) w przedziale \([XS;XE]\).

1. Określamy liczbę węzłów NN.

2. Wyznaczamy węzły Chebyszewa w przedziale bazowym \([-1;1]\).

Nad przedziałem \([-1;1]\) tworzymy półokrąg jednostkowy (o promieniu równym 1):

Na obwodzie półokręgu rozmieszczamy równomiernie \(NN\) węzłów \({h_0, h_1,...,h_{NN-1}}\) w kolejności od prawej do lewej:

Węzły rzutujemy prostopadle na oś OX. W tym celu wystarczy wyliczyć ich współrzędne \({t_k, \; k = 0,1, \dots NN - 1}\):

Pierwszy węzeł wypada w pobliżu 1:

Ostatni węzeł wypada w pobliżu -1:

3. Wyznaczone węzły \(t_k, \; k = 0,1, \dots ,NN-1\) przeliczamy na docelowy przedział interpolacji \([XS;XE]\)

4. Obliczamy:

W ten sposób otrzymujemy współrzędne NN węzłów Chebyszewa \({(x_k;y_k)}\) w przedziale docelowym \({[XS;XE]}\)

5. Teraz konstruujemy wielomian interpolacyjny Chebyszewa stopnia NN-1. Wielomian ten powstaje z wielomianów bazowych, które definiujemy rekurencyjnie:

Oto przykładowe pięć kolejnych wielomianów bazowych Chebyszewa:

Wielomiany bazowe Chebyszewa odnoszą się do przedziału \({[-1;1]}\). Parametr \(t\) jest wartością w przedziale \({[-1;1]}\), która odpowiada wartości \(x\) w przedziale aproksymacji \({[XS;XE]}\). Mając wartość \(x\) musimy ją przeliczyć na \(t\):

Wielomian interpolacyjny Chebyszewa jest sumą tych wielomianów bazowych \({T_j (t), \; j = 0,1, \dots ,NN-1}\) pomnożonych przez odpowiednie współczynniki \(c_j\). Musimy zacząć od wyznaczenia tych współczynników. Dla danego zestawu węzłów wyznaczamy je tylko raz:

Mamy wielomiany bazowe \(T_j(t)\) oraz współczynniki \(c_j\), możemy napisać wzór na wielomian interpolacyjny Chebyszewa:

6. Mając wielomian interpolacyjny Czebyszewa \(P_{NN-1}(t)\) używamy go następująco. Gdy potrzebujemy znaleźć przybliżenie \(f(x)\) w przedziale interpolacyjnym \({[XS;XE]}\), współrzędną \(x\) przeliczamy na \(t\) i wyliczamy wartość wielomianu Czebyszewa \({P_{NN-1}(t)}\). Wynik jest przybliżeniem wartości \(f(x)\).

Schemat Clenshawa

Okazuje się, że wielomianów bazowych Chebyszewa \(T_j (t)\) wcale nie musimy wyliczać. Potrzebne nam są jedynie współczynniki \(c_j\). Dla wielomianów Chebyszewa Schemat Clenshawa jest odpowiednikiem schematu Hornera. Gdy mamy wyliczone współczynniki \({c_j ; j = 0,1, \dots ,NN-1}\) i chcemy wyznaczyć wartość wielomianu Chebyszewa \({P_{NN - 1} (x)}\) dla \({x \in [XS;XE]}\) postępujemy następująco:

1. Wyliczamy położenie \(x\) w przedziale bazowym \({[-1;1]}\):

2. Tworzymy zmienne pomocnicze:

3. Idąc wstecz od \(k = NN - 1, NN-2, \dots ,1\), obliczamy:

4. Na koniec obliczamy wartość wielomianu interpolacyjnego Czebyszewa jako:

i otrzymujemy przybliżenie funkcji \(f(x)\).

Program interpoluje funkcję:

  \[y = \sin(x); x \in \left[-\dfrac{\pi}{2}; \dfrac{\pi}{2} \right]\]

Przedział interpolacji: \(\left[-\frac{π}{2};\frac{π}{2} \right]\). Program wylicza \(NN\) węzłów Czebyszewa w przedziale bazowym \([-1;1]\). Następnie przelicza współrzędne \(t\) węzłów na przedział interpolacji i wyznacza dla nich wartości funkcji \(f(x)\). W kolejnym kroku wylicza współczynniki wielomianu Czebyszewa. Na koniec generuje \(N\) współrzędnych \(x\) w przedziale interpolacji i wylicza dla nich wartość wielomianu interpolacyjnego Chebyszewa przy pomocy schematu Clenshawa. Wynik jest przybliżoną wartością interpolowanej funkcji w punkcie \(f(x)\). Wyświetla wynik i porównuje go do wartości funkcji.

// Interpolacja Chebyszewa
// (C)2026 mgr Jerzy Wałaszek
// Metody numeryczne 0078
//---------------------------

#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;
// Przyrównanie do zera
const double EPS = 1e-12;

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

// Współrzędne węzłów
double nx[NN];
double ny[NN];
// Współczczynniki
double c[NN];
double b[NN+2];

// 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 k;
  double t,x;
  for(k = 0; k < NN; k++)
  {
    t = cos(double(k+k+1) * M_PI / 2 / NN);
    x = (XE - XS) * t / 2 + (XS + XE) / 2;
    nx[k] = x;
    ny[k] = f(x);
  }
}

// Funkcja wylicza współczynniki
// wielomianu Chebyszewa
//------------------------------
void set_c(void)
{
  int j,k;
  double sm;

  for(j = 0; j < NN; j++)
  {
    sm = 0;
    for(k = 0; k < NN; k++)
      sm += ny[k] *
            cos((j*(k+k+1)*M_PI)/2/NN);
    c[j] = (sm + sm) / NN;
  }
}

// Wylicza wartość wielomianu
// interpolacyjnego P(x)
// schematem Clenshawa
//---------------------------
double f_i(double x)
{
  int k;
  double t = (x + x - XS -XE) / (XE - XS);
  b[NN + 1] = b[NN] = 0;
  for(k = NN - 1; k > 0; k--)
    b[k] = c[k] +
           2 * t * b[k + 1] - b[k + 2];
  return c[0] / 2 + t * b[1] - b[2];
}

// 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
  set_c();

  cout << "Interpolacja Chebyszewa\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;
  }

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

Interpolacja Chebyszewa
-----------------------

x = -1.3727, f(x) = -0.9804, interpolacja f(x) = -0.9804
x =  0.3957, f(x) =  0.3854, interpolacja f(x) =  0.3854
x = -1.5372, f(x) = -0.9994, interpolacja f(x) = -0.9994
x =  0.3152, f(x) =  0.3100, interpolacja f(x) =  0.3100
x = -0.9216, f(x) = -0.7966, interpolacja f(x) = -0.7966
x =  1.3860, f(x) =  0.9830, interpolacja f(x) =  0.9830
x =  0.1252, f(x) =  0.1249, interpolacja f(x) =  0.1249
x = -0.1579, f(x) = -0.1572, interpolacja f(x) = -0.1572
x =  0.8288, f(x) =  0.7371, interpolacja f(x) =  0.7371
x =  0.6460, f(x) =  0.6020, interpolacja f(x) =  0.6020

# Interpolacja Chebyszewa
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0078
#---------------------------

from math import sin, cos, pi
from random import uniform

# Definicje stałych
#------------------

# Przedział interpolacji od XS do XE
XS = -pi / 2
XE =  pi / 2
# Ilość węzłów interpolacji
NN = 10
# Ilość punktów interpolowanych
N = 10
# Przyrównanie do zera
EPS = 1e-12

# Tablice
#--------

# Współrzędne węzłów
nx = [0.0] * NN
ny = [0.0] * NN
# Współczczynniki
c = [0.0] * NN
b = [0.0] * (NN + 2)

# Funkcje
#--------

# Funkcja interpolowana
#----------------------
def f(x):
    return sin(x)

# Wyznacza węzły interpolacji
#----------------------------
def set_nodes():
    for k in range(NN):
        t = cos((k + k + 1)*pi/2/NN)
        x = (XE - XS)*t/2 + (XS + XE)/2
        nx[k] = x
        ny[k] = f(x)

# Funkcja wylicza współczynniki
# wielomianu Chebyszewa
#------------------------------
def set_c():
    for j in range(NN):
        sm = 0
        for k in range(NN):
            sm += (ny[k] *
              cos((j*(k+k+1)*pi)/2/NN))
        c[j] = (sm + sm) / NN

# Wylicza wartość P(x)
#---------------------
def f_i(x):
    t = (x + x - XS -XE)/(XE - XS)
    b[NN + 1] = 0
    b[NN] = 0
    for k in reversed(range(1,NN)):
        b[k] = (c[k] +
                2*t*b[k + 1] - b[k + 2])
    return c[0]/2 + t*b[1] - b[2]

# Program główny

# Generujemy węzły
set_nodes()

# Obliczamy współczynniki
set_c()

print("Interpolacja Chebyszewa\n"
      "-----------------------\n")

# Generujemy N punktów x,
# liczymy dla nich y
# i wyświetlamy wynik
for i in range(N):
    x = uniform(XS,XE)
    y = f_i(x);
    print(f"x = {x:7.4f}, "
          f"f(x) = {f(x):7.4f}, "
          f"interpolacja f(x) = {y:7.4f}")

print()
input("Naciśnij Enter...")

Sprawdź działanie algorytmu Chebyszewa dla bardziej skomplikowanych funkcji i dla wielomianu o wyższym stopniu. Zwróć uwagę na dokładność wyników i brak zjawiska Rungego.

Zastanów się nad usprawnieniem programu poprzez eliminację tablicy \({b[\; ]}\) w algorytmie Clenshawa. Co tak naprawdę jest ci potrzebne do wyliczenia wartości wielomianu Czebyszewa?