Układy równań nieliniowych - prosta metoda

Jakie są metody rozwiązywania układów równań nieliniowych?

Przykładowy układ 4 równań:
(x-xi)^2 + (y-yi)^2 - (z-zi) - a*(ti-q)^2 = 0; i = 1..4

szukane: x,y,z,q;
te: xi,yi,zi,ti oraz 'a' są oczywiście dane - znane.

Jak to rozwiązać numerycznie?

Mam na myśli prostą metodę, np. coś w stylu metody bisekcji albo siecznych, tyle że dla układu równań zamiast jednego.

Metodę Newtona można użyć do takiego rozwiązania, ma ona swoją wersję dla układu równań.

To nie jest zbyt prosta metoda, bo tam trzeba obliczać pochodne, co byłoby kłopotliwe algorytmicznie.

x_next = x - F(x)*A^-1

gdzie: A to macierz pochodnych; jakobian chyba:
http://pl.wikipedia.org/wiki/Macierz_Jacobiego

potem trzeba to jeszcze odwracać, ewentualnie rozwiązać układ równań..
dx = F(x)*A^-1
czyli można rozwiązać układ równań (liniowych): Adx = F(x);

Pochodne należy obliczać symbolicznie, tj. przekształcić wzór symbolicznie, uprościć i podstawiać. Jeśli postać równania jest stała, to można to zrobić na etapie programowania tj. ręcznie i wklepać wzory na pochodne wprost w program. Tak to zrobiono np. w symulatorze SPICE. Jeśli chcemy mieć solver uniwersalny, to niestety trochę więcej roboty i trzeba pisać coś co umie parsować równania i je przekształcać symbolicznie. Obliczanie numeryczne wprost z ilorazu różnicowego jest powolniejsze i mniej dokładne i ma tę brzydką cechę, że jeśli weźmiesz za mały krok, to rosną błędy numeryczne, a z kolei przy dużym kroku, zwiększają się błędy związane z nieliniowością funkcji.

Odwracania macierzy nie zaleca się, bo jest mało dokładne. Znacznie lepiej układ liniowy rozwiązywać przez rozkład LU.

A co do tego, czy jest to prosta metoda? W kilkudziesięciu linijkach C / Javy da radę, zakładjąc użycie tylko wbudowanych bibliotek. Zależy jak bardzo chcemy to zoptymalizować. Przy macierzach rzadkich może być więcej kodu, ale przy gęstych i naiwnej implementacji, kodu dużo raczej nie ma.

Dlaczego muszę obliczać pochodne?

Jeśli obliczę ilorazy różnicowe typu centralnego: [f(x+h)-f(x-h)] / 2h,
wynik będzie chyba taki sam dla dostatecznie małych h.

W najgorszym przypadku wyjdzie pewnie coś w stylu metody siecznych dla w wersji n-wymiarowe,
czyli niedużo gorsza od Newtona, ale pochodnych nie muszę już wprost kodować.

Wielki Ogórek napisał(a):

Dlaczego muszę obliczać pochodne?

Nie musisz, ale jest to szybsza i dokładniejsza metoda.

Jeśli obliczę ilorazy różnicowe typu centralnego: [f(x+h)-f(x-h)] / 2h,
wynik będzie chyba taki sam dla dostatecznie małych h.

Właśnie nie będzie taki sam, a jak dobierzesz zbyt małe h, będzie tak daleki od wartości prawidłowej, że całkowicie bezużyteczny.
Ogólnie dobór h jest tutaj najbardziej problematyczny, bo optymalne h będzie zależeć od przebiegu funkcji, i dla jednej funkcji będzie trzeba dobrać h o rzędy wielkości inne niż dla drugiej funkcji.

Przy h dążącym do zera, wynik odejmowania w liczniku będzie też dążył do zera, natomiast błąd bezwzględny odejmowania w liczniku będzie zależny od wielkości argumentów odejmowania. Czyli przy małym h, będziesz mieć w liczniku odejmowane dwie duże liczby, dające w wyniku bardzo małą. Trudno sobie wyobrazić gorszą sytuację jeśli chodzi o błąd numeryczny. Przy dostatecznie małym h, błąd numeryczny przekroczy nawet moduł wartości wyniku, co spowoduje, że Twój wynik stanie się bezużyteczny, tj. może nawet mieć inny znak, co zniszczy zbieżność metody Newtona. Dodam, że odpowiedź na pytanie "jakie h jest już za małe a jakie nie" nie jest trywialna i praktycznie nie da się na nie odpowiedzieć traktując funkcję jako "black box", tj. bez wnikania w strukturę wzoru i analizę błędów numerycznych w obliczaniu funkcji.

Można chyba zastosować do odejmowania schemat typu:

a-b = (sqrt(a) - sqrt(b))*((sqrt(a) + sqrt(b));
co chyba eliminuje problem redukcji cyfr?

ewentualnie obliczamy te różnice w większej precyzji, np. w poczwórnej,
co łatwo zrobić za pomocą par double.

A co ten schemat niby daje? Dalej masz odejmowanie dwóch liczb o zbliżonych modułach.

To chyba eliminuje problem z utratą cyfr.

Weźmy przykład - trzy cyfry precyzji:
a = 1.000, b = 1.001
a-b = -0.001, i tu jest OK.

problem pojawia się dopiero gdy są różne cechy:
a = 1.000, b = 0.999
w tym przypadku cechy są różne zatem trzeba wyrównać, a wtedy: 0.999 = 1.000,
co daje: a-b = 0, zamiast 0.001

Tak mi się przynajmniej wydaje.
Należałoby to sprawdzić empirycznie, czyli wyliczyć te f(x+h) = a, oraz f(x-h) = b, itd.

Poza tym zapominasz o tym, że masz również błąd obliczenia samego a oraz b. Poczwórna precyzja - oczywiście, ale to tylko zmniejsza problem, a nie go rozwiązuje, kosztem znacznego wydłużenia obliczeń.

Odejmowanie na podwójnej precyzji (w stosunku do double) raczej niewiele kosztuje -
z 4 razy tyle co na doble chyba, co jest praktycznie nic, biorąc po uwagę że obliczenia samych równań kosztują tysiące operacji.

problem pojawia się dopiero gdy są różne cechy:
a = 1.000, b = 0.999
w tym przypadku cechy są różne zatem trzeba wyrównać, a wtedy: 0.999 = 1.000,
co daje: a-b = 0, zamiast 0.001

Tak mi się przynajmniej wydaje.
Należałoby to sprawdzić empirycznie, czyli wyliczyć te f(x+h) = a, oraz f(x-h) = b, itd.

Chyba się zgadza z moimi przypuszczeniami:
http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
*
The answer is exactly the same as if the difference had been computed exactly and then rounded. Take another example: 10.1 - 9.93. This becomes

x = 1.01 × 101
y = 0.99 × 101
x - y = .02 × 101

The correct answer is .17, so the computed difference is off by 30 ulps and is wrong in every digit! How bad can the error be? *

Wielki Ogórek napisał(a):

A co ten schemat niby daje? Dalej masz odejmowanie dwóch liczb o zbliżonych modułach.

To chyba eliminuje problem z utratą cyfr.

Cały czas koncentrujesz się na samej operacji odejmowania, a nie bierzesz pod uwagę, że składniki wchodzące do odejmowania są obarczone błędem. Te cyfry utracone zostały już wcześniej, na etapie liczenia a i b. Więc jeśli chcesz to liczyć w poczwórnej precyzji, to samo policzenie odejmowania na poczwórnej precyzji nic nie da. Musiałbyś wszystkie obliczenia robić w poczwórnej precyzji, co tanie wcale nie jest. No tylko zmniejsza problem, a nie go rozwiązuje.

Jeśli już upierasz się przy tej metodzie, to lepiej już wziąć za duże h, niż za małe. Duże h zdegeneruje metodę do czegoś w rodzaju wielowymiarowej metody siecznych. Za małe spowoduje utratę zbieżności bo pewne pochodne będą Ci się zerować lub zmieniać znak.

BTW: Ja nie twierdzę, że metoda ilorazów różnicowych nie działa. Ona będzie pewnie działać w większości sytuacji. Niemniej, wcześniej czy później natrafisz na bardzo wredne funkcje, które nie będą chciały się poddać takiemu numerycznemu rózniczkowaniu, nawet na poczwórnej precyzji. Np. f(x) = 1 + e^(-x). Dla dużych x pochodna będzie dążyła do zera, ale z ilorazu różnicowego będzie Ci wychodziło albo 0, albo totalne bzdury. Natomiast policzenie wprost ze wzoru na pochodną f'(x) = -e^(-x) da bardzo dokładne wartości.

OK. Chodziło mi o metody proste,
zatem ostatecznie mogę użyć metodę siecznych, która i tak będzie chyba szybsza dla złożonych funkcji,
gdzie koszt wyliczania pochodnych jest duży.

2^4 = 16,
1.618^6 = 18

4 iteracje Newtona to tyle co 6 iteracji siecznych.