2011-01-11 16 views
8

To powinno być bardzo proste. Mam funkcję f(x) i chcę ocenić f'(x) dla danego x w MATLAB.jak ocenić pochodną funkcji w programie Matlab?

Wszystkie moje wyszukiwania wyszły z matematyki symbolicznej, która nie jest tym, czego potrzebuję, potrzebuję numerycznego zróżnicowania.

E.g. gdybym określić: fx = inline('x.^2')

Chcę znaleźć f'(3) powiedzieć, co byłoby 6, nie chcę, aby znaleźć 2x

Odpowiedz

7

Aby uzyskać różnicę numerycznej (różnica symetryczna), można obliczyć (f(x+dx)-f(x-dx))/(2*dx)

fx = @(x)x.^2; 
fPrimeAt3 = (fx(3.1)-fx(2.9))/0.2; 

Alternatywnie, można utworzyć wektor wartości funkcji i zastosować DIFF, tj

xValues = 2:0.1:4; 
fValues = fx(xValues); 
df = diff(fValues)./0.1; 

Zauważ, że diff trwa przednia różnica i zakłada, że ​​dx wynosi 1.

Jednak w twoim W takim przypadku lepiej jest zdefiniować fx jako wartość polynomial i oszacować pochodną funkcji, a nie wartości funkcji.

+0

+1 za fajną odpowiedź :) – posdef

+0

Mam rację myśląc, że im mniej wybieram dx, tym dokładniejsza będzie moja odpowiedź. Jeśli tak, to czy istnieje stała MATLAB dla BARDZO małej liczby rzeczywistej? Coś jak pi dla 3.14 ... lub i dla sqrt (-1)? – lms

+2

@codenoob: 'eps' daje bardzo małą liczbę. Jednak dla większości praktycznych celów wystarczy 0,0001. Ponadto, jeśli pracujesz z wielomianami i używasz 'polyder', nie musisz się martwić o rozmiar' dx'. – Jonas

4

Próbowałaś diff (oblicza różnice i przybliża pochodną), gradient lub polyder (oblicza pochodną wielomianu) funkcje?

Możesz przeczytać więcej na temat tych funkcji, korzystając z help <commandname> na konsoli MATLAB lub korzystając z przeglądarki funkcji w menu Pomoc.

+0

+1 za bycie trochę szybszym :) – Jonas

0

Dla danej funkcji w postaci analitycznej, można ocenić pochodną w żądanym punktem z kodu:

syms x 
df = diff(x^2); 
df3 = subs(df, 'x', 3); 
fprintf('f''(3)=%f\n', df3); 

do czystej pochodnych liczbowych wykorzystywać już podanych rozwiązań Jonas i posdef.

6

Z braku symbolicznej skrzynki narzędziowej nic nie stoi na przeszkodzie, aby użyć Derivest, narzędzia do automatycznej adaptacyjnej numeracji różnicowej.

derivest(@sin,pi) 
ans = 
      -1 

Na przykład bardzo ładnie. W rzeczywistości zapewnia nawet oszacowanie błędu w wynikowym przybliżeniu.

fx = inline('x.^2'); 
[fp,errest] = derivest(fx,3) 

fp = 
      6 
errest = 
    3.6308e-14 
+0

To jest interesujące, thx. Zaimplementowałem go za pomocą metody symetrycznej różnicy. – lms

9

Jeśli funkcja jest znany dwukrotnie różniczkowalne, użyj

f'(x) = (f(x + h) - f(x - h))/2h 

który jest drugi w kolejności dokładne godz. Jeśli jest tylko jednokrotnie różny, należy jednoznacznie ustawić wartość w wierszu na wartość 1.

To jest teoria. W praktyce sprawy są dość trudne. Przyjmę drugą formułę (pierwsze zamówienie), ponieważ analiza jest prostsza. Wykonaj drugie zamówienie jako ćwiczenie.

Pierwsza obserwacja jest taka, że ​​musisz musi upewnij się, że (x + h) - x = h, w przeciwnym razie dostaniesz ogromne błędy. Rzeczywiście, f (x + h) i f (x) są blisko siebie (powiedzmy 2.0456 i 2.0467), a kiedy je odejmujesz, tracisz wiele znaczących liczb (tutaj jest to 0,0011, co daje 3 znaczące liczby mniej niż x). Tak więc każdy błąd na h może mieć ogromny wpływ na wynik.

Tak więc, pierwszy krok, napraw kandydata h (pokażę ci po chwili jak go wybrać) i weź h do obliczenia ilość h '= (x + h) - x. Jeśli używasz języka takiego jak C, musisz zadbać o to, aby zdefiniować h lub x jako zmienne, aby nie można było zoptymalizować tych obliczeń.

Następnie wybór h. Błąd w (*) składa się z dwóch części: błędu obcięcia i błędu zaokrągleń. Błąd obcięcia jest fakt, że formuła nie jest dokładny:

(f(x + h) - f(x))/h = f'(x) + e1(h) 

gdzie e1(h) = h/2 * sup_{x in [0,h]} |f''(x)|.

Błąd zaokrąglania wynika z faktu, że f (x + h) i f (x) są blisko siebie. Można oszacować grubsza

e2(h) ~ epsilon_f |f(x)/h| 

gdzie epsilon_f jest względną dokładność obliczeń funkcji f (x) (lub F (x + h), który jest blisko). To musi być ocenione na podstawie twojego problemu. W przypadku prostych funkcji, epsilon_f może być traktowany jako epsilon maszyny. Dla bardziej skomplikowanych może być gorszy niż rząd wielkości.

Więc chcesz h, który minimalizuje e1(h) + e2(h). Podłączając wszystko razem i optymalizacji w h daje

h ~ sqrt(2 * epsilon_f * f/f'') 

która musi być określona na podstawie swojej funkcji. Możesz wziąć przybliżone szacunki. Jeśli masz wątpliwości, weź h ~ sqrt (epsilon) gdzie epsilon = dokładność maszyny. Dla optymalnego wyboru h względna dokładność, z jaką znana jest pochodna, to sqrt (epsilon_f), czyli. połowa znaczących liczb jest poprawna.

W skrócie: zbyt mały a h => błąd zaokrągleń, zbyt duży a h => błąd obcięcia.

Dla wzoru drugiego rzędu, takie same wydajności obliczeń

h ~ (6 * epsilon_f/f''')^(1/3) 

i częściowym dokładność (epsilon_f)^(2/3) dla pochodnej (zwykle jedna albo dwie cyfry znaczące lepsze niż pierwsza formuła zamówienia, zakładając podwójną precyzję).

Jeśli jest to zbyt nieprecyzyjne, poproś o więcej metod, istnieje wiele sztuczek, aby uzyskać lepszą dokładność. Ekstrapolacja Richardsona jest dobrym początkiem dla płynnych funkcji. Ale te metody zazwyczaj są obliczane dość często, może to być lub nie być to, co chcesz, jeśli twoja funkcja jest złożona.

Jeśli zamierzasz używać wielu pochodnych wiele razy w różnych punktach, interesujące jest skonstruowanie przybliżenia Czebyszewa.

+0

Bardzo ciekawa i szczegółowa odpowiedź, dziękuję! Czy mogę zapytać, jakie jest Twoje pochodzenie? – lms

+0

@codenoob: matematyka, głównie. –

Powiązane problemy