Kiedy chce się obliczyć pochodne numeryczne, metoda przedstawiona przez Bengt Fornberg tutaj (i zgłaszane tutaj ) jest bardzo wygodne (zarówno precyzyjne i proste do wykonania). Jako, że oryginalny artykuł pochodzi z 1988 roku, chciałbym wiedzieć, czy istnieje dziś lepsza alternatywa (jako (lub prawie tak) prosta i bardziej precyzyjna)?

Przegląd

Dobre pytanie. Istnieje artykuł zatytułowany „Poprawa dokładności metody różnicowania macierzy dla dowolnych punktów kolokacji” autorstwa R. Baltenspergera. Moim zdaniem nie jest to nic wielkiego, ale ma sens (który był znany już przed pojawieniem się w 2000 r.): Podkreśla znaczenie dokładnego przedstawienia faktu, że pochodna funkcji stałej $f(x)=1$ powinna być równe zero (dotyczy to dokładnie matematyki, ale niekoniecznie w postaci liczbowej).

Łatwo zauważyć, że wymaga to, aby sumy wierszy n-tej macierzy pochodnych $D^{(n)}$ były równe zero. Często wymusza się to ograniczenie, dostosowując przekątną, tj. Ustawiając

$$\tag 1 D^{(n)}_{jj} := -\sum_{\substack{i=1\\i \neq j }}^N D_{ij}\,.$$ Oczywiste jest, że ta funkcja nie działa dokładnie podczas pracy na komputerze z powodu błędów zaokrągleń w obliczeniach zmiennoprzecinkowych. Bardziej zaskakujące jest to, że błędy te są jeszcze poważniejsze, gdy stosuje się wzory analityczne dla macierzy pochodnych (które są dostępne dla wielu klasycznych punktów kolokacji, np. Gaussa-Lobatto).

Teraz artykuł (i zawarte w nim odnośniki) stwierdza, że ​​błąd pochodnej jest rzędu odchylenia sum wierszy od zera. Dlatego celem jest, aby były one jak najmniejsze numerycznie.

Testy numeryczne

Chodzi o to, że procedura Fornberga wydaje się być pod tym względem całkiem dobra. Na poniższym obrazku porównałem zachowanie dokładnej, tj. Analitycznej, pierwszej macierzy pochodnej i tej uzyskanej przez algorytm Fornberga, dla różnej liczby punktów kolokacji Czebyszewa-Lobatta.

Ponownie, wierząc w stwierdzenie w cytowanym artykule, oznacza to, że algorytm Fornberga da bardziej dokładne wyniki dla pochodnej.

Aby to udowodnić, użyję tej samej funkcji, co w pracy,

$$\tag 2 f(x) = \frac{1}{1+x^2}\,.$$ $$\tag 3 E_n = \max_{i\in \{0,\ldots,n\}} \bigg| f^\prime(x_i) - \sum_{j=1}^n D_{ij} f(x_j) \bigg|\,.$$ $$\tag 4 \tilde D_{jj} = D_{jj} - \left(\sum_{i=1}^n D_{ji} \right), \qquad \text{for all}\ j\,.$$

Wniosek

Podsumowując, metoda Fornberga wydaje się dość dokładna, w przypadku nawet o około 3 rzędy wielkości dokładniejsza niż wyniki z wzorów analitycznych. To powinno być wystarczająco dokładne dla większości aplikacji. Co więcej, jest to niezwykłe, ponieważ Fornberg nie wydaje się wyraźnie uwzględniać tego faktu w swojej metodzie (przynajmniej nie ma wzmianki w dwóch artykułach Fornberga).$N=512$

Kolejny rząd wielkości można uzyskać dla tego przykładu poprzez proste włączenie równania (4). Ponieważ jest to dość proste podejście i zastosowane tylko raz dla każdej pochodnej, nie widzę powodu, aby go nie używać.

Metoda z pracy Baltenspergera - która wykorzystuje bardziej wyrafinowane podejście do oceny sumy w równaniu (1) w celu zmniejszenia błędów zaokrągleń - daje w przybliżeniu ten sam rząd wielkości błędu. Przynajmniej w tym przykładzie jest to mniej więcej odpowiednik powyższej metody „Dostosowany Fornberg”.

Zakładając, że próbujesz rozróżnić numeryczną implementację funkcji ciągłej, istnieje wiele metod:

1) Automatyczne różnicowanie. Najdokładniejsza i ogólna metoda. Bolesne dla kodu, wymagające przeciążenia operatora i wyszukiwania zależnego od argumentów. Obciąża użytkownika zrozumieniem tych pojęć. Zmaga się również z wymiennymi osobliwościami, takimi jak różnicowanie cynku przy .$x= 0$

2) Przekształcenie Czebyszewa. Rzuć swoją funkcję na rozpiętość wielomianów Czebyszewa i rozróżnij trzyokresowy nawrót. Super szybki, bardzo dokładny. Wymaga jednak posiadania krótko obsługiwanej domeny zainteresowań; poza wybraną domeną trzyterminowa rekurencja jest niestabilna.$[a,b]$

3) Różnicowanie skończone. Niedoceniany w 1D; zobacz Wskazówki i triki Nicka Highama w dziedzinie obliczeń numerycznych . Chodzi o to, że jeśli zrównoważysz błąd obcięcia i błąd zaokrąglenia, nie musisz wybierać wielkości kroku; można go wybrać automatycznie. W trybie Boost ten pomysł służy do odzyskania (domyślnie) 6/7 poprawnych cyfr dla typu. (Higham pokazuje tylko pomysł na prostszy przypadek 1/2 poprawnych cyfr, ale pomysł można łatwo rozszerzyć.) Współczynniki pochodzą z wyrównanej tabeli Fornberga, ale wielkość kroku jest wybierana przy założeniu, że funkcję można ocenić na 1ULP precyzja. Wadą jest to, że wymaga 2 ocen funkcji, aby odzyskać połowę cyfr tego typu, 4, aby odzyskać 3/4 cyfr itd. W 1D nie jest to zły interes. W wyższych wymiarach jest katastrofalny.

4) Złożona pochodna krokowa. Użyj . Poświęć aby być jednostką zaokrągleń, a to odzyska prawie każdą poprawność. Jest to jednak trochę oszustwo, ponieważ generalnie trudniej jest zaimplementować funkcję w złożonej płaszczyźnie niż ręcznie kodować jej prawdziwą pochodną. Nadal fajny pomysł i przydatny w pewnych okolicznościach.$f'(x) \approx \Im(f(x+ih))$$h$

Nie znam nikogo, kto poprawiłby algorytm Fornberga (patrz także jego nieco nowszy artykuł ). Nawiasem mówiąc, wydaje mi się, że postrzeganie jego algorytmu jako sposobu obliczania pochodnych numerycznych jest niewłaściwe. Wszystko, co zrobił, to opracowanie wydajnego algorytmu do obliczania wag dla metod różnic skończonych. Zaletą tej metody jest to, że daje Ci wagi wszystkich pochodnych do żądanej pochodnej za jednym razem.

Prostszy schemat

Oprócz mojej innej odpowiedzi, która dotyczy bardziej rozszerzenia metody Fornberga, zajmę się tutaj pytaniem o prostsze alternatywy.

W tym celu szkicuję alternatywny schemat, który bardziej bezpośrednio wytwarza współczynniki pochodne interpolacji Lagrangiana. Jego implementacja wymaga tylko kilku wierszy kodu, działa dla dowolnych siatek i według moich pierwszych eksperymentów jest równie dokładna jak Fornberga.

Podstawą implementacji jest wyimaginowana pochodna gdzie jest zmienną rzędu precyzji maszyny. Pochodna wyimaginowanego kroku jest znana z tego, że stabilnie wytwarza wartości pochodnych i nie cierpi na niestabilność liczbową implementacji różnicy skończonej z .

$$f'(x) \ = \ \frac{1}{h}\,\text{Im}\big(f(x+ih)\big)\,,$$ $h$$h\rightarrow 0$

Drugim składnikiem jest wielomian interpolacyjny Lagrange'a na siatce oceniany za pomocą jednej z form barycentrycznych, np. gdzie Aby użyć pochodnej kroku złożonego, należy upewnić się, że te formuły działają również dla złożonych argumenty . Ponadto, dla danej funkcji f (x) i wektora współczynników , oznaczamy wielomian interpolacji poprzez przez $L_i(x)$$\{x_1,\ldots, x_N\}$

$$L_i(z) \ = \ \begin{cases}1 & \text{if }z = x_i\\\ \frac{\frac {\mu_i}{z-x_k}}{ \sum_k\frac {\mu_k}{z-x_k} } & \text{otherwise}\,. \end{cases}$$ $\mu_i = \frac{1}{\prod_{k\neq i} (x_i - x_k)}$$z$$f_i = f(x_i)$$(x_1,f_i)$ $$P(x;f) = \sum_{i=1}^N f_i L_i(x)\,.$$

---

Algorytm

Algorytm został naszkicowany poniżej. Ma te same parametry wejściowe i wyjściowe co w Fornbergu, ale jest znacznie bardziej oczywiste.

Wejście:

• $x$ : siatka z N odrębnymi punktami siatki
• $ord$ : rząd pochodnych
• z: punkt, w którym pochodna ma być oceniana
• Ewentualnie: funkcja lub jego wartości funciton w gridpoints (tylko wymagane dla wariantu wyjściowego 2)$f(x)$$f_i$

Inicjalizacja

• Zainicjuj wielomiany Lagrange'a za pomocą interpolacji barycentrycznej.$L$
• Zainicjuj tablicę -matrices , dla$N\times N$$D^{(o)}$$o=0, \ldots, ord$
• Ustaw , tj. Macierz jednostek wymiaru$D^{(0)} := I_N$$N\times N$
• Ustaw$o := 0$

Algorytm

Podczas gdy :$o < ord$

• Oblicz przez pochodną dla dla wszystko i . Tutaj oznacza ty rząd .$D^{(o+1)}_{ik} = CSD \bigg( P\big(x_k; D^{(o)}_{i}\big) \bigg)$$CSD$$i$$k$
  $D^{(o)}_{i}$$i$$D^{(o)}$

• Ustaw o = o + 1;

Zdecyduj, co wydrukować :

1. Wektor współczynników różnic skończonych w punkcie , gdzie . To właśnie robi Fornberg.$d^{(ord)}$$z$$d_i = P(z;D^{(ord)}_i)$

2. Funkcja interpolacji do pochodnej zamówienia . W tym celu musisz wprowadzić funkcję lub. wartość funkcji w algorytmu.$p^{(ord)}(x) = \sum_{i=1}^N f(x_i) P(x;D^{(ord)}_i)$$f^{(ord)}(x)$$ord$$f_i$$x_i$

3. Meta-funkcja która zwraca funkcję interpolacji wariantu 2., ale dla dowolnej funkcji która ma być interpolowana w punktach siatki.$\text{diff}(\text{int } ord, \text{function } g)$$g$

Osobiście najbardziej podoba mi się wariant 3..

---

Analiza algorytmu

Algorytmem Fornberga jest . Jeśli znajdę czas, opublikuję więcej wyników empirycznych dotyczących dokładności, stabilności itp.$\mathcal O(ord \cdot N^2)$

Aby zwiększyć precyzję różnicowania numerycznego, wykonaj następujące czynności:

1) Wybierz swoją ulubioną wysoce precyzyjną „standardową” metodę opartą na rozmiarze kroku h .

2) Oblicz wartość pochodnej metodą wybraną w 1) wiele razy przy różnych, ale rozsądnych rozmiarach kroku h . Za każdym razem możesz wybrać h jako liczbę losową z przedziału (0,5 * H / 10, 1,5 * H / 10), gdzie H jest odpowiednim rozmiarem kroku dla używanej metody.

3) Średnia wyników.

Twój wynik może zyskać 2-3 rzędy wielkości w bezwzględnym błędzie wrt. wynik nie uśredniony.

https://arxiv.org/abs/1706.10219