2011-01-11 28 views
8

Esto debería ser muy simple. Tengo una función f(x), y quiero evaluar f'(x) para un x dado en MATLAB.cómo evaluar la derivada de la función en matlab?

Todas mis búsquedas han surgido con matemática simbólica, que no es lo que necesito, necesito diferenciación numérica.

E.g. Si defino: fx = inline('x.^2')

Quiero encontrar f'(3) decir, lo que sería 6, no quiero encontrar 2x

Respuesta

7

Para obtener una diferencia numérica (diferencia simétrica), se calcula (f(x+dx)-f(x-dx))/(2*dx)

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

Como alternativa, puede crear un vector de valores de la función y aplicar DIFF, es decir

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

Tenga en cuenta que diff toma la diferencia directa, y que asume que dx es igual a 1.

Sin embargo, en su En este caso, es mejor que defina fx como polynomial, y evalúe la derivada de la función, en lugar de los valores de la función.

+0

+1 por una buena respuesta :) – posdef

+0

Estoy en lo cierto al pensar que cuanto más pequeño elijo dx, más precisa será mi respuesta. Si es así, ¿hay una constante MATLAB para un número real MUY pequeño? Algo como pi para 3.14 ... o i para sqrt (-1)? – lms

+2

@codenoob: 'eps' te da un número muy pequeño. Sin embargo, para la mayoría de los propósitos prácticos, 0.0001 debería ser suficiente. Además, si trabajas con polinomios y usas 'polyder', no tienes que preocuparte por el tamaño de' dx'. – Jonas

4

¿trató diff (calcula las diferencias y se aproxima a un derivado), gradient o polyder (calcula la derivada de un polinomio) funciones?

Puede leer más sobre estas funciones usando help <commandname> en la consola MATLAB, o usar el navegador de funciones en el menú Ayuda.

+0

+1 por ser un poco más rápido :) – Jonas

0

Para una función dada en forma analítica, se puede evaluar la derivada en un punto deseado con el siguiente código:

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

Para los derivados numéricos puros utilizar las soluciones ya dadas por Jonas y posdef.

6

Al carecer de la caja de herramientas simbólica, nada le impide usar Derivest, una herramienta para la diferenciación numérica adaptativa automática.

derivest(@sin,pi) 
ans = 
      -1 

Para su ejemplo lo hace muy bien. De hecho, incluso proporciona una estimación del error en la aproximación resultante.

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

fp = 
      6 
errest = 
    3.6308e-14 
+0

Esto es interesante, gracias. Sin embargo, lo he implementado usando el método de diferencia simétrica. – lms

9

Si su función es conocida por ser dos veces diferenciable, utilice

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

que es de segundo orden preciso en h. Si solo es diferenciable una vez, use

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

que es de primer orden en h.

Esto es teoría. En la práctica, las cosas son bastante complicadas. Tomaré la segunda fórmula (primer orden) ya que el análisis es más simple. Haz el segundo orden uno como ejercicio.

La primera observación es que debe asegúrese de que (x + h) - x = h, de lo contrario obtendrá grandes errores. De hecho, f (x + h) yf (x) están cerca uno del otro (digamos 2.0456 y 2.0467), y cuando los restas, pierdes muchas cifras significativas (aquí es 0.0011, que tiene 3 cifras significativas menos que x). Entonces, cualquier error en h es probable que tenga un gran impacto en el resultado.

Así que, primer paso, arregle un candidato h (le mostraré en un minuto cómo elegirlo), y tome como h para su cálculo la cantidad h '= (x + h) - x. Si está utilizando un lenguaje como C, debe tener cuidado de definir h o x como volátil para que ese cálculo no se optimice.

A continuación, la elección de h. El error en (*) tiene dos partes: el error de truncamiento y el error de redondeo. El error de truncamiento es porque la fórmula no es exacta:

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

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

El error de redondeo proviene del hecho de que f (x + h) yf (x) están cerca uno del otro. Se puede estimar aproximadamente como

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

donde epsilon_f es la relativa precisión en el cálculo de f (x) (o f (x + h), que está cerca). Esto debe evaluarse a partir de su problema. Para funciones simples, epsilon_f se puede tomar como la máquina épsilon. Para los más complicados, puede ser peor que eso en órdenes de magnitud.

Así que usted quiere h que minimice e1(h) + e2(h). Conectar todo juntos y la optimización de los rendimientos en h

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

que ha de estimarse a partir de su función. Puede tomar estimaciones aproximadas. En caso de duda, tome h ~ sqrt (epsilon) donde epsilon = precisión de la máquina. Para la elección óptima de h, la precisión relativa a la que se conoce la derivada es sqrt (epsilon_f), es decir. la mitad de las cifras significativas son correctas.

En resumen: demasiado pequeño a h => error de redondeo, demasiado grande a h => error de truncamiento.

Para la segunda fórmula orden, mismos rendimientos de cálculo

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

y una precisión fraccionaria de (epsilon_f)^(2/3) para el derivado (que es típicamente una o dos cifras significativas mejor que el primero fórmula de orden, suponiendo doble precisión).

Si esto es demasiado impreciso, no dude en pedir más métodos, hay muchos trucos para obtener una mayor precisión. La extrapolación de Richardson es un buen comienzo para funciones suaves. Pero esos métodos generalmente se calculan bastantes veces, esto puede ser o no lo que usted desea si su función es compleja.

Si va a utilizar derivadas numéricas muchas veces en diferentes puntos, resulta interesante construir una aproximación de Chebyshev.

+0

¡Respuesta muy interesante y detallada, gracias! ¿Puedo preguntar cuál es tu experiencia? – lms

+0

@codenoob: matemática, principalmente. –

Cuestiones relacionadas