[C++] Implementacja interpolacji Newtona

[karpiuch]

No więc: dostałem za zadanie zaimplementować interpolację Newtona w języku C++, no i jak policzyć na kartce to potrafię, tak nie mam za wielkiego pomysłu na implementację. Dodam, że mam się obejść bez dwuwymiarowej tablicy.

Do wyliczania ilorazów różnicowych stworzyć tablicę o wielkości \(\displaystyle{ \frac{n \cdot \left( n-1\right) }{2}}\)? Czy zapisywać to jakoś umiejętnie do zmiennej? Nie mogę znaleźć jakieś zależności, która pomogłaby mi to jakoś wszystko zaimplementować.

Z góry dziękuję za pomoc.

[szw1710]

Z napisaniem programu sam sobie poradzisz - programistą nie jestem.

Wyliczenie ilorazów rzędu \(\displaystyle{ k}\) wymaga jedynie ilorazów rzędu \(\displaystyle{ k-1,}\) więc o reszcie można zapomnieć. Dlatego potrzebujesz tylko dwóch ciągów \(\displaystyle{ (x_1,\dots,x_n)}\) jako węzłów oraz \(\displaystyle{ (y_1,\dots,y_n)}\) jako wartości. W pierwszym kroku zapiszesz
\(\displaystyle{ \begin{aligned}
y_1&:=\frac{y_2-y_1}{x_2-x_1}\\
y_2&:=\frac{y_3-y_2}{x_3-x_2}\\
y_3&:=\frac{y_4-y_3}{x_4-x_3}
\end{aligned}}\)
itd. Zauważ, że nie będziesz tracił wartości, bo po podmianie \(\displaystyle{ y_1}\) starego \(\displaystyle{ y_1}\) już nie trzeba itp. Kolejne kroki podobnie. Działaj.

Drugi krok - początek

\(\displaystyle{ \begin{aligned}
y_1&:=\frac{y_2-y_1}{x_3-x_1}\\
y_2&:=\frac{y_3-y_2}{x_4-x_2}\\
y_3&:=\frac{y_4-y_3}{x_5-x_3}
\end{aligned}}\)

Na każdym kroku ubywa po jednym \(\displaystyle{ y}\)-ku. Pierwszy krok - zdefiniujesz do \(\displaystyle{ y_{n-1}}\), zaś \(\displaystyle{ y_n}\) (stare) zostaje jako śmieć. Na drugim kroku do \(\displaystyle{ y_{n-2}}\) itp.

Jeśli węzłami są kolejne liczby naturalne \(\displaystyle{ 1,2,\dots,n}\), to potrzeba nam jedynie znać ciąg \(\displaystyle{ y}\)-ków.

[karpiuch]

Czyli lepszym wyjściem jest nadpisywanie tablicy \(\displaystyle{ y}\), tak? Po każdym kroku zmniejszam "wielkość" tablicy o 1. Tylko jak to uwzględnić sensownie w obliczeniach? Właśnie tą utratę jednego \(\displaystyle{ y'ka}\) przy każdym kroku?

[szw1710]

Nie tyle utratę, ile śmieć. Na drugim kroku nie korzystasz już z \(\displaystyle{ y_n}\). Na ostatnim kroku szukanym ilorazem jest \(\displaystyle{ y_1}\), a cała reszta to ponadpisywane śmiecie. Do policzenia samych ilorazów wystarczą więc dwa ciągi. Aby jednak zastosować wzór interpolacyjny, musimy pamiętać \(\displaystyle{ y_1}\) z każdego kroku, bo te będą się nam gubić w krokach kolejnych. Dlatego potrzebny będzie jeszcze jeden ciąg. Albo jeszcze lepiej. Bierzesz zmienną pomocniczą \(\displaystyle{ z:=y_1}\) i po pierwszym kroku przypisujesz \(\displaystyle{ y_n:=z}\), bo stare \(\displaystyle{ y_n}\) to śmieć i nie jest potrzebne, można je wykorzystać do szczytnych celów. Drugi krok tak samo: \(\displaystyle{ z:=y_1}\) i po zakończeniu \(\displaystyle{ y_{n-1}:=z}\) Tak więc wystarczy w sumie jeden ciąg i nie będzie żadnych śmieci. Wzór interpolacyjny Newtona zastosujesz tak:

\(\displaystyle{ \begin{aligned}w(x)&=y_n+y_{n-1}(x-x_1)+y_{n-2}(x-x_1)(x-x_2)+\dots\\&+y_1(x-x_1)(x-x_2)\dots(x-x_{n-1}).\end{aligned}}\)

Przećwicz to sobie na kartce interpolując funkcję (dla przykładu) \(\displaystyle{ f(x)=2^x}\) w węzłach \(\displaystyle{ -1,0,1,2}\). Każdorazowo zapisuj ciąg \(\displaystyle{ y}\)-ków.

[karpiuch]

Zabrałem się obecnie za funkcję liczącą iloraz różnicowy,

Kod: Zaznacz cały

for (i = 0; i < n; i++)

{

for (j = n; j > i; j--)

{		

y[j] = (y[j] - y[j - 1]) / (x[j] - x[j - i - 1]);

n--;

}
}

Czy to w ogóle ma prawo dobrze liczyć te ilorazy? Wydaje mi się, że gdzieś popełniłem tutaj błąd, ale nawet nie mam pojęcia gdzie.

e; o matko, brak umiejętności obsługi icode, przepraszam.

[szw1710]

To jest dopiero pierwszy krok. Ale zapomniałeś o zmiennej pomocniczej \(\displaystyle{ z}\). Definiujesz ją przed tą pętlą, a przypisanie \(\displaystyle{ y_n:=z}\) robisz po jej zakończeniu.

Za programowanie się nie biorę. Język C++ to dla mnie egipskie hieroglify. Poza Pascala nigdy nie wyszedłem (może to nie do końca prawda, bo ostatnio piszę proste skrypty w \(\displaystyle{ R}\)). Ale ponieważ dobrze znam się na teorii interpolacji Newtona, wymyśliłem algorytm, który - moim zdaniem - powinien zadziałać. Nie bierz się od razu do programowania, nie siadaj do komputera. To błąd programistów stary jak świat światem. Rozpisz to sobie na kartce jak prosiłem. Zobaczysz czy mój opis działa. A jestem niemal pewien, że jest poprawny. Ale zawsze jest jakiś margines błędu.

[karpiuch]

Dobrze, spróbuję to sobie teraz rozpisać, ale nie rozumiem tego etapu z użyciem zmiennej pomocniczej: w pierwszym kroku \(\displaystyle{ y _{n} := z}\), a w ostatnim \(\displaystyle{ y_{n-1} := z}\), czemu?
+ rozumiem, że przy każdym kroku zmienną pomocniczą ustawiam na pierwszy wyraz ciągu, tak? Następnie tą zmienną pomocniczą ustawiam wartość \(\displaystyle{ y}\), który w następnym kroku już potrzebny nie będzie, a w poprzednim go używałem, tak?

[szw1710]

Tak. Inaczej dane poznikają. Kolejność jest bardzo ważna.

Zmienna pomocnicza jest potrzebna, jeśli chcesz interpolować. Musimy znać po pierwszym ilorazie z każdego kroku. Tak będzie najlepiej. Nie ma śmieci, wszystko jest wykorzystane i masz współczynniki wielomianu raz na zawsze. Umiałbym jednak pominąć i zmienną pomocniczą. Np. po pierwszym kroku definiujesz \(\displaystyle{ w(x):=y_1}\). Po drugim kroku \(\displaystyle{ w(x):=w(x)+y_1(x-x_1)}\). Po trzecim \(\displaystyle{ w(x):=w(x)+y_1(x-x_1)(x-x_2)}\) itp. Ale chyba zdecydowałbym się na tę zmienną pomocniczą, a nie na definiowanie wielomianu interpolacyjnego na raty.

[karpiuch]

Na kartce to wygląda super, a siadając ponownie do kodu nadal bez pomysłu jak zapętlić te ilorazy różnicowe.

[szw1710]

Pętla będzie podwójna. Jeden indeks na krok, drugi na ilorazy w kroku.

Kod: Zaznacz cały

for i:=1 to n-1
 z:=y_1
 for j:=1 to n-i+1
  {Tu liczenie ilorazów}
 y_{n-i+1}:=z

Nie dbam o poprawne indeksy na końcu ale mniej więcej o to chodzi.

A teraz programuj. Idę na wieczorny wykład.

[karpiuch]

Próbowałem robić to na przykładzie, który Pan podał (ten dotyczący \(\displaystyle{ 2^{x}}\).. No i dla sprawdzenia postanowiłem sobie wypisać tablicę końcową z ilorazami różnicowymi, niestety wypisuje mi całą tablice wartościami \(\displaystyle{ 0.5}\), a tak to nie powinno wyglądać.
Trochę edytowałem kod, no i wydaje mi się, że wszystko jest zapisane dobrze, ale jednak coś nie tak...

Kod: Zaznacz cały

for i := 0 to i := n
  z = y[0]
  for j = 1 to j := n - i + 1
    y[j] = (y[j] - y[j - 1]) / (x[i + j] - x[j - 1])
  y[n - i - 1] = z

Nie mam już pojęcia gdzie popełniam błąd, sprawdzałem to kilka dobrych razy z pomocą kartki...

e; działają ilorazy, błędem było:

Kod: Zaznacz cały

 y[j] = (y[j] - y[j - 1]) / (x[i + j] - x[j - 1])

Przypisywałem to do wyrazu \(\displaystyle{ j}\) zamiast \(\displaystyle{ j-1}\).

[szw1710]

I co: działa?

[karpiuch]

Tak jak ilorazy działają, tak teraz problemem z lekka jest sama interpolacja.

[szw1710]

Podałem Ci jak to trzeba zrobić. Nawet dwukrotnie. Samego wzoru na symbolach nie ma sensu generować. Wybieramy jakieś \(\displaystyle{ x}\) i piszemy program na obliczenie

\(\displaystyle{ \begin{aligned}w(x)&=y_n+y_{n-1}(x-x_1)+y_{n-2}(x-x_1)(x-x_2)+\dots\\&+y_1(x-x_1)(x-x_2)\dots(x-x_{n-1}).\end{aligned}}\)

U ciebie indeksy są chyba o jeden mniejsze - od \(\displaystyle{ 0}\) do \(\displaystyle{ n-1}\) albo do \(\displaystyle{ n}\). Ale zasada jest jasna - idziemy od końcowego \(\displaystyle{ y}\) i to jest stała. Bierzemy poprzednie i mnożymy przez \(\displaystyle{ (x-x_1)}\) (węzły numerujesz od \(\displaystyle{ x_1}\) czy od \(\displaystyle{ x_0}\), jeśli od \(\displaystyle{ x_0}\), to \(\displaystyle{ (x-x_0)}\)?). Itd. itp. jak mówią początkowe wyrazy.

[karpiuch]

Wydaje mi się, że zrobiłem:

Kod: Zaznacz cały

for i:=n - 1 to i:=0
  w:=1
  for j:=0 to j:=k
    w*= (arg - x[j])
  w*= p[i];
  sum += w;
  k++;

Gdzie \(\displaystyle{ k}\) iteruję od zera, a \(\displaystyle{ arg}\) to argument dla którego mam znaleźć wartość.
Czy to prawidłowy algorytm? Na razie dla wyników wychodzi dobrze.

Bardzo dziękuję za pomoc, przy okazji zrozumiałem dość dobrze tą Interpolację.