Oblicz

13

Funkcja ma osobliwość zbliżoną do x = 0 . Osobliwość tę można jednak znieść: dla x = 1 należy mieć f ( x ) = 1 , ponieważ e x = k = 0 x kf:x(ex1)/xx=0x=1f(x)=1 a zatem (ex-1)/x=k=1x k - 1

ex=k=0xkk!
Jednak forma(ex-1)/xnie jest nie tylko zdefiniowana przyx=0, ale jest również niestabilna numerycznie w pobliżu tego punktu; w celu ocenyf(x)dla bardzo małychxliczbowo, można zastosować rozszerzenie Taylora, tj. obcięcie wyżej wymienionych szeregów mocy.
(ex1)/x=k=1xk1k!
(ex1)/xx=0f(x)x

P : Czy funkcja ma nazwę? Innymi słowy, czy jest to powszechny problem?f

P : Czy ktoś wie o bibliotece C / C ++, która ładnie radzi sobie z tą sytuacją, tj. Używa rozszerzenia Taylora o odpowiednim stopniu w pobliżu 0 i innej reprezentacji od zera?

anonimowy
źródło

Odpowiedzi:

19

expm1ex1x=0

n00b
źródło
17

Jest to przypadek błędu anulowania. Biblioteka standardowa C (od C99) zawiera funkcję o nazwie, expm1która pozwala uniknąć tego problemu. Jeśli użyjesz expm1(x) / xzamiast (exp(x) - 1.0) / x, nie wystąpi ten problem (patrz wykres poniżej). <code> fabs (expm1 (x) / x - (exp (x) - 1.0) / x) </code>

Szczegóły i rozwiązanie tego konkretnego problemu zostały szczegółowo omówione w rozdziale 1.14.1 Dokładności i stabilności algorytmów numerycznych . To samo rozwiązanie wyjaśniono również na stronie 19 artykułu W. Kahana pt. „ Jak daremne są bezmyślne oceny zaokrągleń w obliczeniach zmiennoprzecinkowych? . Rzeczywista implementacja expm1w bibliotece GNU C różni się od podejścia opisanego w odnośnikach powyżej i jest dokładnie udokumentowana w kodzie źródłowym .

Juan M. Bello-Rivas
źródło
1
Dziękuję, właśnie tego potrzebowałem! Niestety mogę zaakceptować tylko jedną odpowiedź ...
anonimowy
Oczywiście! Nie ma problemu :-)
Juan M. Bello-Rivas
3

Aby odpowiedzieć na twoje pierwsze pytanie, nie, funkcja nie ma nazwy (a przynajmniej takiej, która jest powszechnie znana).

Jak wspomnieli inni, najlepszym sposobem na obliczenie tej funkcji jest leczenie kilku szczególnych przypadków. W ten sposób każda biblioteka obliczy tę funkcję.

  1. Przypadek 0: x = 0, zwraca 1.
  2. |x|<δ1+x/2δdouble2e-85e-4
  3. W innym przypadku: powrót expm1(x)/x.

Możesz być bardziej wyrafinowany i specjalny przypadek więcej rzeczy ze ściętą serią Taylor, ale prawdopodobnie nie jest tego wart. W rzeczywistości nie jest do końca jasne, że przypadek 1 należy rozpatrywać osobno, ponieważ, jak wskazano w k20, anulowanie jest bezpieczne. Jednak osobne obchodzenie się z nim pozwoliłoby mi być bardziej pewnym siebie.

Victor Liu
źródło
2

Pamiętam to pytanie, które zostało zadane wcześniej na tej stronie, i zaskakująco odpowiedź brzmi, że potrzebujesz tylko równości równej zeru w specjalnym przypadku. Błędy anulują się w pobliżu zera. Nie mam linku

Tak, ta odpowiedź była całkowicie błędna. Nie jestem pewien, dlaczego został tak wysoko oceniony, prawdopodobnie dlatego, że został napisany tak autorytatywnie. Znalazłem łącze, które miałem na myśli. Było to tutaj wymiana stosu matematyki , a nie na wymianie stosu Scicomp. expm1-Darmowy formuła anulowanie błędu podano w odpowiedzi JM i wykorzystuje u = exp(x)transformację.

k20
źródło
xdx(edx1)/dx(1+dx1)/dx1
1
dx1+dx=1
0

Aby odpowiedzieć na pierwsze pytanie i podać (prawdopodobnie liczbowo nieefektywną) metodę drugiego, zauważ, że jest to odwrotność funkcji generującej liczby Bernoulliego .

Nikolaj-K
źródło
To interesujące połączenie, dzięki za zwrócenie na to uwagi. Niestety, sądzę, że potrójna suma sprawi, że będzie to zbyt drogie. Co więcej, nie jest od razu jasne, gdzie należy obciąć każdą sumę, aby uzyskać pożądaną dokładność.
anonimowy
@anonymous: Jaką sumę potrójną masz na myśli? Nie potrzebujesz wielomianów Bernoulliego, tylko liczby Bernoulliego i możesz je wymienić z góry. Ale tak, wciąż nie może być lepiej niż seria Taylor.
Nikolaj-K,
Możesz je obliczyć z wyprzedzeniem, jeśli jest oczywiste, że potrzebujesz tylko skończonej liczby dla dowolnego wejścia.
anonimowy
@anonymous: Cóż, tak, jakbyś z góry podał współczynniki Taylora.
Nikolaj-K