Pochodna funkcji wielu zmiennych

Witam!

Czy ktoś wie jak napisać program obliczający pochodną cząstkowa funkcji (dla przykładu) dwóch zmiennych w konkretnym punkcie? Tak, że nie potrzebna jest nam znajomość wzoru tej pochodnej, a jedynie wynik.

Obliczenie pochodnej funkcji jednej zmiennej jest banalne, dzięki skorzystaniu z definicji, lecz przy dwóch zmiennych już zaczynają się schody.

Zależy mi żeby opracować algorytm liczący gradient funkcji dla określonego punktu. Np. Podajemy x0=[1,1] i program podaje nam wartości pochodnych cząstkowych dla tych punktów.

Czy jest to możliwe bez pisania skomplikowanego programu liczącego wzór tej pochodnej?

Podejrzewam, że trzeba się oprzeć na liczeniu pochodnej z definicji, i jakoś to zmodyfikować, tak by liczyła wspomniane pochodne cząstkowe, lecz niestety nie potrafię tego zrobić.

Pozdrawiam!

Problem sprowadza się do pochodnej funkcji jednej zmiennej.

Weźmy dwa wymiary. Liczymy gradient funkcji f(x,y) w punkcie P = (PX, PY).

Najpierw liczysz w tym punkcie pochodną funkcji jednowymiarowej powstałej z przecięcia f płaszczyzną y=PY.
Nazwijmy ją g, jej wzór to po prostu
g(x) = f(x, PY)

Podobnie liczysz dla kierunku y, czyli masz funkcje
h(y) = f(PX, y)

Te dwie pochodne dają gradient.

ujmujac obraozwo..

zalozmy ze masz z = f(x,y)

z def. wartosc pochodnej w x0 = lim_e->0 ( F(x0+e) - F(x0) ) / e
czyli, badasz jak wzrasta wartosc w momencie gdy argument 'wahnie' sie delikatnie

gradient z kolei to zestaw czastkowych pochodnych po kolejnych zmiennych, czyli tutaj: [df/dx, df/dy]
nie ma w nim pochodnych po dwoch zmiennych
badajac wartosc gradientu w punkcie, defacto szukasz wartosci pochodnych df/dx i df/dy w tym punkcie - czyli patrzysz jak *f wzrasta przy wahnieciach x (y sie nie rusza), *jak f wzrasta przy wahnieciach y (x sie nie rusza).

stad, wartosc pochodnych w punkcie (x0,y0):
(df/dx)(x0,y0) = lim_e->0 ( f(x0+e, y0) - f(x0, y0) ) / e
(df/dy)(x0,y0) = lim_e->0 ( f(x0, y0+e) - f(x0, y0) ) / e

niestety, nie obliczysz lim nie przeksztalcajac algebraicznie wzorow, wiec musisz sie ograniczyc do numerycznego: wybierasz dostatecznie male, konkretne 'E' i zakladasz ze
g(x0,y0) = ( f(x0+E, y0) - f(x0, y0) ) / E ~= (df/dx)(x0,y0)
h(x0,y0) = ( f(x0, y0+E) - f(x0, y0) ) / E ~= (df/dy)(x0,y0)
pozostaje podstawic i wyliczyc.

np.
z = f(x,y) = x*y^{2 (oczywiscie, df/dx = y}2, df/dy = 2xy)
x0,y0 = (1,2)
E = 0.01
g(x0,y0) = f(1.01,2)-f(1,2))/0.01 = (4.04-4)*100 = 4 (algebraicznie: (df/dx)(1,2) = 4)
h(x0,y0) = f(1,2.01)-f(1,2))/0.01 = (0.0401-4)*100 = 4.01 (algebraicznie: (df/dy)(1,2) = 4)

czyli obliczylismy ze grad w (1,2) to (4,4.01). algebraicznie wytyczylismy ze (4,4). no coz.. przyjelismy sobie E=0.01. E jest > 0, to nie lim_e->0, wiec rozbieznosc sila rzeczy bedzie tym wieksza im bardziej E>>0

Wzory podali przedmówcy, za to większy problem przy numerycznym liczeniu pochodnych będziesz mieć z przyjęciem odpowiedniego E. Dasz za małe, będziesz mieć duży błąd zaokrągleń. Przyjmiesz za duże, będziesz mieć duży błąd obcięcia.

Ok a co jeśli pochodna na podstawie 2 punktów jest zbyt mało dokładna?? Czy można zwiększyć dokładność obliczając pochodną na podstawie np. 4 punktów?? Staram sie znaleźć gradient funkcji stress Kruskala i algorytm dwupunktowy daje mało ciekawe rezultaty. Metoda BFGS minimalizacji funkcji zbiega sie zbyt wolno. Podobno jest jakiś wzór przybliżający gradient tej funkcji. Zna ktoś może??

Z moich doświadczeń wynika, że dla BFGS i funkcji z małą liczbą parametrów (poniżej 100) w miarę dobrze sprawuje się pochodna centralna:

Możesz też spróbować pochodnej 5-punktowej, chociaż jest ona dość kosztowna:

Jednak osobiście radziłbym zastosować gotowe rozwiązanie: http://www.chokkan.org/software/liblbfgs/. Biblioteka napisana w C.