Interpolacja

Interpolacja Newtona jest kolejną z metod interpolacji wielomianowej. Polega ona na tworzeniu wielomianu, który przechodzi przez wszystkie węzły interpolacji. Przy budowie tego wielomianu interpolacja Newtona wykorzystuje rekurencję i tzw. różnice dzielone, które tworzą proces rekursywnego dzielenia. Różnice dzielone wykorzystywano dawniej do wyliczania tablic logarytmów i funkcji trygonometrycznych. Algorytm różnic dzielonych zastosował Charles Babage (zwany Ojcem Komputerów) w swojej maszynie różnicowej w roku 1822:

W algorytmie interpolacji Newtona wykorzystywane są tzw. ilorazy różnicowe tworzone rekurencyjnie. Na prostym przypadku zobaczmy, jaka jest idea tworzenia różnic. Utwórzmy tabelę, w której pierwszej kolumnie znajdują się kwadraty kolejnych liczb naturalnych. Kolumnę nazwijmy różnicami stopnia 0, r₀.

Następna kolumna powstaje z różnic wyrazów kolumny poprzedniej: od wyrazu leżącego w tym samym wierszu w sąsiedniej kolumnie odejmujemy wyraz leżący nad nim.

(?) oznacza wartość nieokreśloną (pomyśl dlaczego?). Zwróć uwagę, iż z kwadratów zrobił nam się ciąg kolejnych liczb nieparzystych. Liczymy dalej:

Różnice stopnia drugiego to wartość stała, równa 2 (bo różnica dwóch kolejnych liczb nieparzystych zawsze jest równa 2). Dwie pierwsze wartości są nieokreślone, bo dla nich nie da się policzyć różnic.

Różnice stopnia trzeciego są zerowe. Dalej nie ma już potrzeby liczyć, bo otrzymamy tylko same zera.

Gdy znamy już ideę tablicowania różnic, możemy wrócić do interpolacji Newtona. Wykorzystuje ona tzw. ilorazy różnicowe. Prześledźmy sposób postępowania. W przedziale interpolacyjnym mamy cztery węzły interpolacyjne (muszą być różne, inaczej otrzymamy dzielenie przez 0) o współrzędnych:

Węzły powinny być równomiernie rozłożone.

Tworzymy ilorazy różnicowe rzędu 0:

Wartość ilorazu różnicowego rzędu 0 f[x_i], to po prostu wartość współrzędnej y_i węzła v_i:

Iloraz różnicowy rzędu 1 powstaje jako różnica ilorazów różnicowych stopnia 0 podzielonych przez różnicę współrzędnych x:

(?) - tego nie liczymy, bo nie da się utworzyć ilorazu różnicowego.

Liczymy ilorazy różnicowe stopnia 2. Powstają one identycznie, jak ilorazy różnicowe stopnia 1:

Została nam ostatnia kolumna:

Elementy przekątnej są współczynnikami wielomianu interpolacyjnego Newtona.

Dla czterech węzłów wielomian interpolacyjny Newtona ma postać:

Program interpoluje funkcję:

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

Przedział <-^π/₂;^π/₂>. W przedziale program generuje NN równoodległych współrzędnych x, wylicza dla nich współrzędną y wg przepisu funkcji i dostaje współrzędne NN węzłów interpolacyjnych. Na ich podstawie wyznacza współczynniki wielomianu interpolacyjnego Newtona. Następnie generuje N pseudolosowych punktów w przedziale <XP;XE> i   przy pomocy wielomianu interpolacyjnego Newtona wylicza wartości współrzędnych y tych punktów, po czym wyświetla wyniki.

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

#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 interpolacji
const int NN = 6;
// Liczba interpolowanych
// punktów
const int N = 10;

// 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()
{
  // Wyznaczamy NN węzłów
  // w równych odległościach
  int i;
  double dx = (XE - XS) / (NN - 1);
  for(i = 0; i < NN; i++)
     nx[i] = XS + dx * i;

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

// Oblicza wartość
// interpolowaną za pomocą
// schematu Hornera
//------------------------
double f_i(double x)
{
  int i;
  double r = a[NN-1];
  for(i = NN - 2; i >= 0; i--)
    r = r * (x - nx[i]) + a[i];
  return r;
}

// Wylicza współczynniki wielomianu
// interpolacyjnego Newtona
//---------------------------------
void set_a()
{
  int i, j;

  // Kopiujemy wartości początkowe
  // do tablicy współczynników
  for(i = 0; i < NN; i++)
    a[i] = ny[i];

  // Liczymy w miejscu od końca
  for(j = 1; j < NN; j++)
    for(i = NN - 1; i >= j; i--)
      a[i] = (a[i] - a[i-1]) /
             (nx[i] - nx[i-j]);
}

// 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
  // wielomianu Newtona
  set_a();

  cout << "Interpolacja wielomianowa "
          "Newtona\n"
          "--------------------------"
          "-------\n\n";

  // Generujemy N punktów x
  // i liczymy dla nich y
  // przy pomocy wielomianu
  // Newtona
  // 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 wielomianowa Newtona
---------------------------------

x = -1.5117, f(x) = -0.9983, interpolacja f(x) = -0.9981
x =  1.2545, f(x) =  0.9504, interpolacja f(x) =  0.9502
x = -0.6943, f(x) = -0.6399, interpolacja f(x) = -0.6399
x = -0.0012, f(x) = -0.0012, interpolacja f(x) = -0.0012
x =  1.5358, f(x) =  0.9994, interpolacja f(x) =  0.9993
x =  1.1873, f(x) =  0.9274, interpolacja f(x) =  0.9272
x = -0.9111, f(x) = -0.7902, interpolacja f(x) = -0.7902
x =  0.0206, f(x) =  0.0206, interpolacja f(x) =  0.0206
x = -0.4180, f(x) = -0.4059, interpolacja f(x) = -0.4060
x = -1.0989, f(x) = -0.8907, interpolacja f(x) = -0.8906

# Interpolacja wielomianowa
# Newtona
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0076
#---------------------------

from math import sin, pi
from random import uniform

# STAŁE
#------

# Przedział interpolacji
# od XS do XE
XS = -pi / 2
XE =  pi / 2
# liczba węzłów interpolacji
NN = 6
# Liczba interpolowanych
# punktów
N = 10

# Tablice
#--------

# Współrzędne x węzłów
nx = []
# Współrzędne y węzłów
ny = []
# Współczynniki wielomianu
a = []

# Funkcje
#--------

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

# Funkcja wyznacza węzły
# dla wielomianu
#-----------------------
def set_nodes():
    # Wyznaczamy NN węzłów
    # w równych odległościach
    dx = (XE - XS) / (NN - 1)
    for i in range(NN):
         nx.append(XS + dx * i)

    # Dla współrzędnych x
    # liczymy współrzędne y
    for i in nx:
        ny.append(f(i))

# Oblicza wartość
# interpolowaną za pomocą
# schematu Hornera
#------------------------
def f_i(x):
    r = a[NN-1]
    for i in reversed(range(NN - 1)):
        r = r * (x - nx[i]) + a[i]
    return r

# Wylicza współczynniki wielomianu
# interpolacyjnego Newtona
#---------------------------------
def set_a():
    # Kopiujemy wartości początkowe
    # do tablicy współczynników
    for i in ny:
        a.append(i)

    # Liczymy w miejscu od końca
    for j in range(1,NN):
        for i in reversed(range(j,NN)):
            a[i] = ((a[i] - a[i-1]) /
                    (nx[i] - nx[i-j]))

# Program główny
#---------------

# Generujemy węzły
set_nodes()

# Obliczamy współczynniki
# wielomianu Newtona
set_a()

print("Interpolacja wielomianowa "
      "Newtona\n"
      "--------------------------"
      "-------\n")

# Generujemy N punktów x
# i liczymy dla nich y
# przy pomocy wielomianu
# Newtona
# 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) = "
          f"{y:7.4f}")

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