Página web del curso

Derivación numérica

In [1]:
import numpy as np
import matplotlib.pyplot as plt

Derivada primera de funciones de una variable

Buscamos obtener una aproximación de $f'$ en un punto $x$ en función de valores de $f$ en un número finito de puntos cercanos a $x$.

¿Por qué?

  • Porque sólo conocemos los valores de $f$ en un número finito de puntos $x$.
  • Porque es menos costoso calcular la aproximación que el valor exacto.

A partir de la definición de derivada

$$ f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}{h} $$

podemos obtener tres aproximaciones: para $h>0$

Fórmula progresiva

$$ f'(x)\approx \frac{f(x+h)-f(x)}{h}$$

Fórmula regresiva

$$ f'(x)\approx \frac{f(x)-f( x -h)}{h}$$

Fórmula centrada

$$ f'(x)\approx \frac{f(x+h)-f(x-h)}{2h}$$

Tener en cuenta que si solo usamos una de estas fórmulas no podemos calcular la derivada al menos en un extremo del intervalo. Por ejemplo, si usamos la fórmula progresiva no podemos calcular la derivada aproximada en $b$ porque necesitamos $f(b+h)$ que, en general, es desconocido.

Ejercicio 1

Dada la función $f(x)=\ln(x)$ en $[1,2]$ y para h = 0.1 y h = 0.01.

  • Calcular y dibujar su derivada exacta.
  • Calcular y dibujar, sobre la figura anterior, su derivada aproximada usando diferencias finitas progresivas, df_p, regresivas, df_r, y centradas, df_c en puntos equiespaciados, con espacio entre ellos h, del intervalo $[1,2]$. En los cálculos no se pueden usar puntos que estén fuera de este intervalo.
  • Calcular y dibujar, en una figura nueva, el error en cada punto para cada una de las aproximaciones, $\left|f'(x_i)-df_a(x_i)\right|$, donde $x_i$ son los puntos y $a=p,~r$ o $c$.
  • Calcular el error global relativo de cada aproximación usando la fórmula
$$ Error_{Global\,rel}=\frac{\left\Vert f'(x_i)-df_a(x_i)\right\Vert_2}{\left\Vert f'(x_i)\right\Vert_2}=\sqrt{\frac{\sum_{i=1}^n (f'(x_i)-df_a(x_i))^2}{\sum_{i=1}^n (f'(x_i))^2}}.$$

Notas

  • Los puntos equiespaciados se pueden obtener, por ejemplo:
    • Usando un bucle.
    • Usando el comando np.arange.
    • Usando el comando np.linspace.
  • Los resultados numéricos pueden variar ligeramente en los últimos 3 o 4 decimales del error dependiendo del método usado.
  • Puede ser útil la orden np.diff(v) que calcula $(v_2-v_1,v_3-v_2,\ldots,v_n-v_{n-1})$
  • Para calcular la norma $2$ (euclidea) de un vector v se puede usar el comando norm. Por ejemplo
    from numpy.linalg import norm
vn = norm(v)
    $$ \mathrm{norm}(v)=\left\Vert v\right\Vert_2= \sqrt{\sum_{i=1}^n v_i^2}$$
In [6]:
%run Ejercicio1.py

h = 0.1 
ErrorG_p = 0.03772359200890384 
ErrorG_r = 0.03854271699344694 
ErrorG_c = 0.001942510922541181

h = 0.01 
ErrorG_p = 0.0038146660730015015 
ErrorG_r = 0.0038228493284284126 
ErrorG_c = 2.0608592494300922e-05

Vemos que la fórmula centrada da lugar a un error menor. Esto se debe a que las fórmulas progresiva y regresiva son de orden 1 y la fórmula centrada es de orden 2.

¿Es posible conseguir aproximaciones de orden 2 para $f'(a)$ y $f'(b)$? Sí, si usamos un polinomio de interpolación de $f$ de $2^{o}$ grado en $x_{0}$, $x_{1}$, $x_{2}$ y $h \gt0 $ con $x_1= x_0+h$ y $x_2 = x_1+h$, siendo $y_{j}=f(x_{j})$. Una vez obtenido el polinomio de derivación, lo derivamos y sustituyendo los valores $x_{0}$ obtenemos la fórmula progresiva, sustituyendo $x_{1}$, la fórmula centrada y sustituyendo $x_{2}$ obtendremos la fórmula regresiva. Todas ellas de orden dos.

$$ \begin{array}{l} \mathbf{Progresiva}\\ f'(x_{0})\approx\dfrac{-3y_{0}+4y_{1}-y_{2}}{2h}\\ \mathbf{Centrada}\\ f'(x_{01})\approx\dfrac{-y_{0}+y_{2}}{2h}\\ \mathbf{Regresiva}\\ f'(x_{2})\approx\dfrac{y_{0}-4y_{1}+3y_{2}}{2h} \end{array} $$

Que también podemos escribir de la forma:

$$ \begin{array}{l} \mathbf{Progresiva}\\ f'(x)\approx \dfrac{1}{2h}\left(-3f(x)+4f(x+h)-f(x+2h)\right)\\ \mathbf{Centrada}\\ f'(x)\approx \dfrac{1}{2h}\left(f(x+h)-f(x-h)\right)\\ \mathbf{Regresiva}\\ f'(x)\approx \dfrac{1}{2h}\left(f(x-2h)-4f(x-h)+3f(x)\right) \end{array} $$

Ejercicio 2

Dada la función $f(x)= x^3 + x ^2 + x$ en $[0,1]$, y h = 0.02

  • Calcular su derivada aproximada utilizando la fórmula centrada para sus puntos interiores y progresiva y regresiva de orden 1 (dos puntos) para $f'(0)$ y $f'(1)$ respectivamente.
  • Calcular su derivada aproximada utilizando la fórmula centrada para sus puntos interiores y progresiva y regresiva de orden 2 (tres puntos) para $f'(0)$ y $f'(1)$ respectivamente.
  • Calcular la derivada exacta en $[0,1]$ y el error global relativo para las dos aproximaciones.
  • Dibujar la derivada exacta y la aproximada

Nota

Para unir vectores se puede, por ejemplo

  • Inicializar un vector y rellenarlo con un bucle.
  • Inicializar un vector y rellenarlo con slicing. Por ejemplo
    v = np.zeros(5)
v[0] = 3.
v[1:4] = np.array([5.,4,3])
v[-1] = 2.
  • Usar el comando np.concatenate
In [8]:
%run Ejercicio2.py

Error con derivación de orden 1 = 0.003427865337711315

Error con derivación de orden 2 = 0.00012590491500409013

Derivada segunda de funciones de una variable

Dados $x_{0}$, $x_{1}$, $x_{2}$ con $x_1= x_0+h$ y $x_2 = x_1+h$ ($h \gt0 $), siendo $y_{j}=f(x_{j})$, una fórmula de aproximación de la derivada segunda es

$$f''(x_1)\approx \frac{1}{h^2}(y_0-2y_1+y_2)$$

o también

$$f''(x)\approx \frac{f(x-h)-2f(x)+f(x+h)}{h^2}$$

Ejercicio 3

Dada la función $f(x)=\cos(2\pi x)$ en $[-1,0]$ y h = 0.001

Calcular y dibujar su derivada segunda exacta en $[-1,0]$ y su derivada segunda aproximada utilizando la fórmula anterior. Calcular el error relativo de la aproximación en el intervalo $[-1+h,-h]$.

In [9]:
%run Ejercicio3.py

Error relativo = 3.2898637920427984e-06

Ejercicios propuestos

Derivadas de funciones de dos variables

Sea $\Omega =[a,b]\times[c,d]\subset \mathbb{R}^2$. Introducimos particiones uniformes

\begin{eqnarray*} a=x_1<x_2<\ldots<x_{M-1}<x_M=b \\ c=y_1<y_2<\ldots<y_{N-1}<y_N=d \end{eqnarray*}

Si llamamos a un punto $(x_m,y_n)$ sus vecinos son

Podemos construir las fórmulas de forma análoga a como hicimos con las fórmulas para una variable. Por ejemplo

  • Derivada parcial (progresiva):

$\partial_x f(x_m,y_n)$

$$ (\delta_x^+f)(x_m,y_n)= \frac{f(x_{m+1},y_n)-f(x_{m},y_n)}{h_x} $$
  • Gradiente (centrado):

$\nabla f(x_m,y_n)=(\partial_x f(x_m,y_n),\partial_y f(x_m,y_n))$

$$ \left(\frac{f(x_{m+1},y_n)-f(x_{m-1},y_n)}{2h_x}, \frac{f(x_{m},y_{n+1})-f(x_{m},y_{n-1})}{2h_y}\right)\quad (2) $$
  • Divergencia (centrada):

$\mathrm{div} (\mathbb{f}(x_m,y_n))=\partial_x f_1(x_m,y_n)+\partial_y f_2(x_m,y_n)$

$$ \frac{f_1(x_{m+1},y_n)-f_1(x_{m-1},y_n)}{2h_x}+ \frac{f_2(x_{m},y_{n+1})-f_2(x_{m},y_{n-1})}{2h_y} \quad (3) $$
  • Laplaciana (centrada):

$\Delta f(x_m,y_n)=\partial_{xx} f(x_m,y_n)+ \partial_{yy} f(x_m,y_n)$

$$ \frac{f(x_{m+1},y_n)-2f(x_m,y_n)+f(x_{m-1},y_n)}{h_x^2}+ \frac{f(x_{m},y_{n+1})-2f(x_m,y_n)+f(x_{m},y_{n-1})}{h_y^2} \quad (4) $$

Pero antes de usar estas fórmulas, recordemos como se representan las funciones de dos variables, tanto escalares como vectoriales, en python.

Se representan las funciones de dos variables en coordenadas cartesianas en dominios rectangulares. Como en el caso de las funciones de una variable nosotros decidimos el número de puntos que utilizaremos en la representación.

Si queremos representar la función escalar $f(x,y)=x^2+y^2$ en $x \in [0,5]$ con paso $h_x=1$ e $y \in [-2,2]$ con paso $h_y=1$, primero creamos la malla

In [2]:
h = 1.
x = np.arange(0., 5. + h, h)
y = np.arange(-2., 2. + h, h)
x, y = np.meshgrid(x, y)

y hemos creado dos matrices

In [3]:
print(x)

[[0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]
 [0. 1. 2. 3. 4. 5.]]
In [4]:
print(y)

[[-2. -2. -2. -2. -2. -2.]
 [-1. -1. -1. -1. -1. -1.]
 [ 0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.]
 [ 2.  2.  2.  2.  2.  2.]]

Calculamos ahora los valores de $z=f(x,y)=x^2+y^2$ para cada par de $(x_{ij},y_{ij}) \quad i=1,\ldots,m,\,j=1,\ldots,n$ siendo $m$ y $n$ el número de filas y columnas tanto de x como de y. Operando elemento a elemento obtenemos una matriz z de dimensiones $m\times n$

In [5]:
z = -(x**2 + y**2)
print(z)

[[ -4.  -5.  -8. -13. -20. -29.]
 [ -1.  -2.  -5. -10. -17. -26.]
 [ -0.  -1.  -4.  -9. -16. -25.]
 [ -1.  -2.  -5. -10. -17. -26.]
 [ -4.  -5.  -8. -13. -20. -29.]]

Con estas tres matrices x, y, z ya podemos representar la función $f$. Por ejemplo, si representamos la función usando las curvas de nivel

In [6]:
import matplotlib.cm as cm # colormap

plt.figure()
plt.contour(x, y, z, 20, cmap=cm.bwr)
plt.title(r'$f(x,y)=-x^2-y^2$', fontsize=16)
plt.show()
In [7]:
plt.figure()
plt.contourf(x, y, z, 20, cmap=cm.bwr)
plt.colorbar()
plt.title(r'$f(x,y)=-x^2-y^2$', fontsize=16)
plt.show()

También podemos representar la función como una superficie

In [9]:
from mpl_toolkits import mplot3d

fig = plt.figure(figsize=(10,5))
ax = plt.axes(projection ='3d') 

ax.plot_surface(x, y, z, cmap=cm.bwr)
plt.title(r'$f(x,y)=-x^2-y^2$', fontsize=16)
plt.show()
In [10]:
fig = plt.figure(figsize=(10,5))
ax = plt.axes(projection ='3d') 

ax.contour3D(x, y, z, 50, cmap=cm.bwr)
plt.title(r'$f(x,y)=-x^2-y^2$', fontsize=16)
plt.show()

Si queremos representar una función vectorial $\mathbb{f}=(f_1(x,y),f_2(x,y))$ podemos usar el comando quiver

In [13]:
f1 = lambda x, y: (x**2 + y**2)/3
f2 = lambda x, y: (x**2 - y**2)/2
plt.figure()
ax = plt.gca()
ax.quiver(x,y,f1(x,y),f2(x,y))
ax.set_xlim([-1,6])
ax.set_ylim([-3,3])
plt.show()

La orden numpy

fy, fx = np.gradient(f,hx,hy)

devuelve el gradiente (calculado numéricamente) de la función escalar f, usando los pasos hx y hy. Lo calcula usando una fórmula centrada de orden dos para los puntos interiores y diferencias de orden uno o dos, progresivas o regresivas en los bordes.

Por lo tanto, la matriz que devuelve tiene las mismas dimensiones que las matrices de los argumentos de entrada.

Ejercicio 4

Calcular el gradiente de una función usando la orden numpy gradient:

  • Dibujar a función $f(x,y)=x e^{-x^2-y^2}$ en $[-2,2]\times[-2,2]$ con hx = 0.1 y hy = 0.1.
  • Calcular el gradiente de $f$ y dibujarlo en la misma figura con la función $f$.
  • Calcular el error global relativo de la aproximación usando np.linalg.norm.
In [12]:
%run Ejercicio4

Error relativo = 0.012265

Ejercicio 5

  • Calcular la divergencia de $\mathbb{f}(x,y)=(f_1(x,y),f_2(x,y))(x^2+y^2,x^2-y^2)$ en $[-2,2]\times[-2,2]$ utilizando las salidas de la función gradient para las funciones $f_1$ y $f_2$, con hx = 0.1 y hy = 0.1. Dibujar en la misma figura $\mathbb{f}$ y $\text{div} (\mathbb{f})$
  • Se aprecia un efecto en los bordes debido a una mala aproximación de la divergencia en estos puntos. Rehacer los dibujos restringiendo la representación de la divergencia a $[-2+h_x,2-h_x]\times[-2+h_y,2-h_y].$
  • Calcular el error global relativo en ambos casos usando np.linalg.norm.
In [11]:
%run Ejercicio5

Error relativo 1 = 9.333167e-03
Error relativo 2 = 9.847095e-16

Ejercicio 6

  • Calcular, usando la fórmula (4), la laplaciana de $f(x,y)=x e^{-x^2-y^2}$ en $[-2+h_x,2-h_x]\times[-2+h_y,2-h_y]$ con hx = 0.1 y hy = 0.1. Dibuja la función y su Laplaciana.
  • Calcular la solución exacta y el error global relativo usando np.linalg.norm.
In [12]:
%run Ejercicio6

Error relativo = 5.734265e-03