Bezpośrednia konwersja wielomianu interpolacyjnego w postaci Newtona na postać Beziera

Jednym z zadań jakie miałem do zrobienia wraz z przyjacielem  na pracownię (a potem była na ćwiczeniach) z analizy numerycznej dotyczyła konwersji wielomianu w postaci Newtona na postać Beziera. Ciężko coś znaleźć na ten temat w internecie i trochę ciężko było, dlatego postanowiłem napisać tego posta ;)
Napisaliśmy konwersję bezpośrednią jak i pośrednią (z wykorzystaniem bazy potęgowej) - opiszę w tym poście tylko tą pierwszą. Ładnie sformatowany pdf (który zawiera również kod) można pobrać tutaj.

Opis problematyki

Niech dany będzie wielomian interpolacyjny w postaci Newtona $L_n$, przy czym $b_k$ oznacza stałe współczynniki a $p_k$ iloraz różnicowy. Poszukujemy jego zapisu w postaci Béziera $B_n$. Wprowadźmy oznaczenia:
$$\begin{align*} L_n(x) = \sum_{k=0}^{n}b_k p_k(x) \\ B_n(x) = \sum^n_{k=0}c_kB^n_k(x) \end{align*}$$
Chcemy by:
$$\begin{equation} \label{eq:rownosc} L_n(x) = B_n(x) \end{equation}$$
Zatem problem sprowadza się to do problemu dobrania takich stałych $c_k$, by ta równość zachodziła.

Szkic rozwiązania

Bazą rozwiązania będzie uogólniony schemat Hornera, który definiujemy w ten sposób:
$$\begin{cases} \label{eq:horner} w_n(x) = b_n \\ w_k(x) = (x - x_k)w_{k+1}(x) + b_{k} \\ L_n(x) = w_0(x) \end{cases}$$
Korzystając z równości wielomianów oraz powołując się na uogólniony schemat Hornera otrzymujemy, że:
$$\begin{align} B^i_k &= (x - x_k)B^{i+1}_k(x) + b_{k} \label{eq:bezhorner} \end{align}$$
Przy czym tutaj indeks górny $i \in [n, n - 1,\dots,0]$ w $B^i_k$ oznacza, że $B^i_k \in \Pi_{n-i}$ (tj. dany wielomian Bernsteina jest $(n - i)$-tego stopnia).

Idea algorytmu polega na tym, by stopniowo budować wielomian $k$-tego stopnia w postaci Béziera na podstawie obliczonych wartościach współczynników wielomianu $(k-1)$-stopnia znalezionego w poprzednim kroku aż do otrzymania ostatecznie wielomianu $n$-tego stopnia.


Wyznaczanie współczynników

Znając już zarys pomysłu, znajdźmy sposób wyrażania współczynników w danym kroku. Wprowadźmy oznaczenie $c_i^j$ oznaczające $i$-ty współczynnik wielomianu Béziera przy $(n-j)$-tej iteracji algorytmu.
Ze schematu hornera otrzymujemy, że $c_n^n = b_n = c_n$, oraz korzystając z nowego oznaczenia zapiszemy teraz  poprzednią równość jako:
$$\begin{equation} \label{eq:sumy} \sum_{i=k-1}^{n}c_i^{n-k-1}B_{i-k+1}^{n-k+1} = (x-x_{k-1})(\sum_{i=k}^{n} c_i^{n-k} B_{i-k}^{n-k}) + b_{k-1} \end{equation}$$
Musimy rozpisać prawą stronę tak, by dopasować zakres sumowania oraz stopnie wielomianów względem lewej strony:
$$\begin{equation*} (x-x_{k-1})(\sum_{i=k}^{n} c_i^{n-k} B_{i-k}^{n-k}) + b_{k-1} = x\sum_{i=k}^{n} c_i^{n-k} B_{i-k}^{n-k} - x_{k-1}\sum_{i=k}^{n} c_i^{n-k} B_{i-k}^{n-k} + b_{k-1} \end{equation*}$$
Wielomiany Bernsteina stanowią podział jedności, tj $\sum_{k=0}^{n} B^n_k(x) = 1$ co umożliwia zapisanie $b_{k-1}$ jako sumę:
$$\begin{equation*} = \sum_{i=k}^{n} c_i^{n-k} (x B_{i-k}^{n-k}) - x_{k-1} \sum_{i=k}^{n} c_i^{n-k} B_{i-k}^{n-k} + \sum_{i=k-1}^{n} b_{k-1} B_{i-k+1}^{n-k+1} \end{equation*}$$
Korzystając z tego, że $xB^n_k = \frac{k+1}{n+1}B^{n+1}_{k+1}$ oraz zauważając, że dla $i = k - 1$ ułamek się wyzeruje, możemy również bezkarnie zmienić zakres sumowania:
$$\begin{equation*} = \sum_{i=k-1}^{n} c_i^{n-k} \frac{i-k+1}{n-k+1}B_{i-k+1}^{n-k+1} - x_{k-1} \sum_{i=k}^{n} c_i^{n-k} B_{i-k}^{n-k} + \sum_{i=k-1}^{n} b_{k-1} B_{i-k+1}^{n-k+1} \end{equation*}$$
Pozostaje jedynie zająć się środkową sumą - podniesiemy stopień wielomianu korzystając z tego, że $B(x) = (1-x)B(x) + xB(x)$ i stosując wcześniej używane metody:
$$\begin{align*} -x_{k-1}\sum_{i=k}^{n} c_i^{n-k} B^{n-k}_{i-k} &= -x_{k-1} \left[ x \sum_{i=k}^n c_i^{n-k} B^{n-k}_{i-k} + (1 - x) \sum_{i=k}^n c_i^{n-k} B^{n-k}_{i-k} \right] \\ &= -x_{k-1} \sum_{i=k-1}^{n} c_i^{n-k} \frac{i-k+1}{n-k+1}B_{i-k+1}^{n-k+1} - x_{k-1}\sum_{i=k-1}^n \frac{n-i}{n-k+1}c_{i+1}^{n-k} B_{i-k+1}^{n-k+1} \end{align*}$$
czyli mamy:
$$\begin{align*} &= (1 - x_{k-1})\sum_{i=k-1}^{n} c_i^{n-k} \frac{i-k+1}{n-k+1}B_{i-k+1}^{n-k+1} - x_{k-1} \sum_{i=k-1}^{n} \frac{n-i}{n-k+1}c_{i+1}^{n-k} B_{i-k+1}^{n-k+1} + \sum_{i=k-1}^{n} b_{k-1} B_{i-k+1}^{n-k+1} \\ &= \sum_{i=k-1}^{n} (1 - x_{k-1}) c_i^{n-k} \frac{i-k+1}{n-k+1}B_{i-k+1}^{n-k+1} - x_{k-1} \frac{n-i}{n-k+1}c_{i+1}^{n-k} B_{i-k+1}^{n-k+1} + b_{k-1} B_{i-k+1}^{n-k+1} \end{align*}$$
podstawiając otrzymaną równość i wyciągając wspólny czynnik:
$$\begin{equation*} \sum_{i=k-1}^{n}c_i^{n-k-1}B_{i-k+1}^{n-k+1} = \sum_{i=k-1}^{n}\left[(1-x_{k-1})\frac{i-k+1}{n-k+1}c_i^{n-k} - x_{k-1}\frac{n-i}{n-k+1}c_{i+1}^{n-k+1} + b_{k-1}\right]B_{i-k+1}^{n-k+1} \end{equation*}$$
z tego wynika, że:
$$\begin{equation*} c_i^{n-k-1} = (1-x_{k-1})\frac{i-k+1}{n-k+1}c_i^{n-k} - x_{k-1}\frac{n-i}{n-k+1}c_{i+1}^{n-k+1} + b_{k-1} \end{equation*}$$
Można to ostatecznie dużo ładniej zapisać w taki sposób:
$$\begin{equation*} c_i^{n-k} = \frac{(1-x_k)(i-k)c_i^{n-k-1} - x_{k}(n-i)c_{i+1}^{n-k-1}}{n-k} + b_{k} \end{equation*}$$
Teraz bardzo łatwo jest napisać program, który obliczy poszukiwane współczynniki $c_i = c_i^0$, za pomocą uzyskanego wzoru ;)