Można interpolować na dwa sposoby, albo obliczać bezpośrednio z definicji wartość w każdym punkcie, albo obliczać od razu wielomian a później z niego korzystać. To drugie jest trudniejsze, więc zakładam że chodzi o pierwsze (i jest to to co Ty próbujesz najwyraźniej robić).
Czyli jednak o to drugie chodziło.
W takim razie z wikipedii:
Wielomian lagrange'a to suma (ściśle, liniowa kombinacja) wielomianów:
gdzie lj to:
(trochę wiecej na ten temat na wspomnianej wiki, nie ma co tego przepisywać.)
Teraz jeśli chcemy wyliczyć ten wielomian z definicji, możemy po prostu wyliczyć wszystkie lj a później je odpowiednio dodać.
Tym razem inne dane:
xs = [1, 2, 3, 4]
ys = [3, 1,-1, 2]
Potrzebna będzie pomocnicza funkcja do dodawania list:
def addv(vec_a, vec_b):
return [a + b for (a, b) in zip(vec_a, vec_b)]
Dodaje ona po prostu odpowiadające sobie elementy obu list:
>>> addv([1, 1, 1], [2, 3, 4])
[3, 4, 5]
Jeszcze rzecz o której nie napisałem - wielomian np. x^3 + 2x^2 + 5x - 2 jest reprezentowany najprościej jak się da, czyli listą [1, 2, 5, -2] (każdy element reprezentuje odpowiedni współczynnik).
I teraz faktyczny kod. Najpierw obliczanie lj:
def lagrange_basis(xs, j):
coeffs, denominator = [1.0], 1.0
for m in range(len(xs)):
if m != j:
coeffs = addv(coeffs + [0], [0] + [-xs[m] * v for v in coeffs])
denominator *= xs[j] - xs[m]
return [c / denominator for c in coeffs]
Patrz:
Obliczanymy tu niezależnie licznik (coeffs) i mianownik (denominator) (można by najpierw obliczyć mianownik i ustawić jako wartość początkową licznika, ale to IMO mniej czytelne), po czym je dzielimy żeby uzyskać wartość ułamka.
Mianownik tego ułamka to po prostu (xj - x0) * (xj - x1) * (xj - x2) ... (xj - xk). W obliczaniu go nie ma nic skomplikowanego chyba (po prostu dla każdego xm mnożymy go przez (xj - xm).
Licznik jest trochę bardziej skomplikowany, bo zawiera mnożenie wielomianów (mnożymy (x - x0) * (x - x1) ... (x - xk)). Na szczeście wszystkie wielomiany przez które mnożymy są postaci (x - n), więc korzystamy z:
(ax^4 + bx^3 + cx^2 + dx + e) * (x - n) = x*(ax^4 + bx^3 + cx^2 + dx + e) - n*(ax^4 + bx^3 + cx^2 + dx + e)
Niezbyt odkrywcze, wiem, ale jako że to ax^4 + bx^3 + cx^2 + dx + e to u nas tablica coeffs, to wynik takiego mnożenia wynosi właśnie addv(coeffs + [0], [0] + [-xs[m] * v for v in coeffs]).
([0] + w [0] + [-xs[m] * v ...] jest konieczne, inaczej dodawane tablice byłyby różnej długości).
Mając tą funkcję, obliczenie wielomianu sprowadza się do dodania pomnożonych przez odpowiedni współczynnik wielomianów:
def interpolate(xs, ys):
coeffs = [0 for x in xs]
for i in range(len(xs)):
basis = lagrange_basis(xs, i)
coeffs = addv(coeffs, [ys[i] * c for c in basis])
return coeffs
Wszystko chyba jest dość oczywiste, więc niezbyt jest co tłumaczyć - po prostu obliczamy dla każdego i lagrange_basis(xs, i), mnożymy przez ys[i] (jak we wzorze), i dodajemy.
Łącząc wszystko:
xs = [1, 2, 3, 4]
ys = [3, 1,-1, 2]
def addv(vec_a, vec_b):
return [a + b for (a, b) in zip(vec_a, vec_b)]
def lagrange_basis(xs, j):
coeffs, denominator = [1.0], 1.0
for m in range(len(xs)):
if m != j:
coeffs = addv(coeffs + [0], [0] + [-xs[m] * v for v in coeffs])
denominator *= xs[j] - xs[m]
return [c / denominator for c in coeffs]
def interpolate(xs, ys):
coeffs = [0 for x in xs]
for i in range(len(xs)):
basis = lagrange_basis(xs, i)
coeffs = addv(coeffs, [ys[i] * c for c in basis])
return coeffs
print interpolate(xs, ys)
Otrzymujemy wynik wyglądający na prawidłowy:
$ python test.py
[0.8333333333333333, -5.0, 7.166666666666666, 0.0]
Do porównania z:
Ofc można to wypisać czytelniej, czyli jako
0.8333*x^3 - 5*x^2 + 7.1666*x
Kwestia gustu/wymageń.
Sprawdź czy tym razem o to Ci chodziło.