Numerycznie stabilna ocena sqrt (x + a) - sqrt (x)

Question

Czy istnieje elegancki sposób numerycznego stabilnego oceniania następującego wyrażenia dla pełnego zakresu parametrów x, a> = 0?

f(x,a) = sqrt(x+a) - sqrt(x) 

jest również tam dowolny język programowania lub biblioteki, które nie zapewniają tego rodzaju funkcji? Jeśli tak, pod jaką nazwą? Nie mam konkretnego problemu z użyciem powyższego wyrażenia już teraz, ale spotkałem się z nim wiele razy w przeszłości i zawsze uważałem, że ten problem musiał zostać rozwiązany wcześniej!

Answer 1

Tak, jest! Pod warunkiem, że co najmniej jedna z x i a jest pozytywna, można użyć:

f(x, a) = a/(sqrt(x + a) + sqrt(x)) 

który jest doskonale numerycznie stabilny, ale nie warto funkcją biblioteki w sobie. Oczywiście, gdy x = a = 0, wynik powinien być 0.

Objaśnienie: sqrt(x + a) - sqrt(x) jest równa (sqrt(x + a) - sqrt(x)) * (sqrt(x + a) + sqrt(x))/(sqrt(x + a) + sqrt(x)). Teraz pomnóż pierwsze dwa terminy, aby uzyskać sqrt(x+a)^2 - sqrt(x)^2, co upraszcza do a.

Oto przykład wykazania trwałości: kłopotliwe przypadku pierwotnego wyrażenia gdzie x + a i x są bardzo blisko wartości (lub równoważnie, gdy a jest znacznie mniejszy pod względem wielkości niż x). Na przykład, jeśli x = 1 i a jest mały, wiemy, że z rozszerzenia Taylor około 1, że sqrt(1 + a) powinien być 1 + a/2 - a^2/8 + O(a^3), więc sqrt(1 + a) - sqrt(1) powinien być bliski a/2 - a^2/8. Spróbujmy tego dla konkretnego wyboru małego a. Oto oryginalna funkcja (napisany w Pythonie, w tym przypadku, ale można traktować go jako Pseudokod):

def f(x, a): 
    return sqrt(x + a) - sqrt(x) 

i tu jest stabilna wersja:

def g(x, a): 
    if a == 0: 
     return 0.0 
    else: 
     return a/((sqrt(x + a) + sqrt(x)) 

Teraz zobaczmy, co mamy z x = 1 i a = 2e-10:

>>> a = 2e-10 
>>> f(1, a) 
1.000000082740371e-10 
>>> g(1, a) 
9.999999999500001e-11 

wartość powinna mamy jest (do dokładności maszyny): a/2 - a^2/8 - na ten szczególny a , warunki sześcienne i wyższe są nieistotne w kontekście podwójnej precyzji IEEE 754, które zapewniają tylko około 16 cyfr dziesiętnych. Załóżmy obliczyć tę wartość dla porównania:

>>> a/2 - a**2/8 
9.999999999500001e-11