Jak uniknąć katastrofalnego anulowania funkcji Pythona?

13

Mam problem z implementacją funkcji numerycznie. Cierpi na tym fakt, że przy dużych wartościach wejściowych wynik jest bardzo dużą liczbą razy bardzo małą liczbą. Nie jestem pewien, czy katastrofalne anulowanie jest właściwym terminem, więc proszę mnie poprawić, jeśli tak jest. Dowody na coś nie tak:

wprowadź opis zdjęcia tutaj

Jak mogę uniknąć oscylacji i przypisania wartości 0,0 dla większych danych wejściowych wynoszących 6?

Oto moja funkcja: 

import numpy as np

def func(x):
    t = np.exp(-np.pi*x)
    return 1/t*(1-np.sqrt(1-t**2))
python numerical-analysis floating-point Dipol
źródło

Odpowiedzi:

31

Jest to rzeczywiście nazywane katastrofalnym anulowaniem. W rzeczywistości ten szczególny przypadek jest bardzo łatwy: przepisać funkcję, używając równoważnego, stabilnego liczbowo wyrażenia Ponieważ prawdopodobnie potrzebujesz odniesienia, jest to omówione w większości podręczników metod numerycznych w odniesieniu do wzoru na rozwiązywanie równań kwadratowych (ten wzór w postaci standardowej jest niestabilny numerycznie dla niektórych wartości parametrów). Aby pozbyć się należy pomnożyć i podzielić przez . Oto porównanie:1-√

t1+1−t2.
1+ √1−1−t21+1−t2

wprowadź opis zdjęcia tutaj

Cyryl
źródło
Fantastyczny! Czy możesz polecić jedną z takich książek, czy te techniki zostały przedstawione?
Dipole
2
@Jack „Dokładność i stabilność algorytmów numerycznych” to dobra książka wyższego poziomu. Każdy podręcznik wprowadzający również to omówi.
Kirill,
Chciałbym wiedzieć, czy użyłeś Wolfram Mathematica do narysowania tego wykresu. THX :)
xyz
Czy znasz jakieś odniesienia gromadzące i / lub omawiające podobne sztuczki w celu przepisania wyrażeń matematycznych w matematycznie równoważny sposób, które zmniejszają utratę znaczenia? Czytam książkę Highama, ale dyskusja jest ogólna, a wszystkie dalsze rozdziały dotyczą algebry liniowej (która obecnie nie jest moim tematem).
becko
@becko Z mojego doświadczenia wynika, że ​​jest to dość ad hoc. Jest to o wiele łatwiejsze, jeśli masz sposób na przetestowanie formuły z poprawnymi odpowiedziami (nawet jeśli po prostu generujesz je z bardzo precyzyjną arytmetyką), abyś nie szukał niestabilności numerycznej bez wcześniejszych nieudanych przypadków testowych. A jeśli działa dla wszystkich znanych danych wejściowych, nie ma prawdziwego problemu, czy niestabilność numeryczna występuje gdziekolwiek, czy nie.
Kirill,