Aproximación de raíces, puntos fijos y optimización de funciones en una variable

Vamos a practicar diversas técnicas para encontrar ceros (raíces) de funciones no lineales de una variable real.

Comenzamos por un ejemplo sencillo: $$ f(x) = x + \frac{1}{2}\sin(x) -1 $$

  • Comenzamos por dibujar la función cuyas raíces queremos calcular.
In [43]:
%format long es necesario para poder ver los números con más decimales
format long;
In [44]:
f = @(x)(x + sin(x)/2 - 1);

xmin = -5;
xmax = 5;
xs = linspace(xmin, xmax,200);
ys = f(xs);
plot(xs, ys);
title('Grafica de f')
xlabel('x')
ylabel('y')

Observar las funciones nos puede dar pistas sobre el número y ubicación aproximada de las raíces, pero necesitamos demostrar la existencia y unicidad con alguna de estas técnicas:

  • Si la función es continua, y en un intervalo cambia de signo, tiene que haber una raíz en ese intervalo (teorema de Bolzano).
  • Si en un intervalo la derivada no cambia de signo, no puede haber más de una raíz (consecuencia del teorema de Rolle).

En este caso concreto, tenemos un cambio de signo en el intervalo (-1,1), luego hay una raíz de $f$:

In [45]:
f(-1)
f(1)

ans = -2.42073549240395
ans =  0.420735492403948
In [46]:
hold on;
plot(xs,ys)
xleft = -1;
xright = 1;

plot([xmin, xmax], [0,0], 'k-')
plot([xleft, xright], [0,0], 'o')
axis xy
hold off;

La derivada de $f$ es fácil de calcular: $$ f'(x) = 1 + \frac{1}{2}\cos(x) $$ y podemos acotarla fácilmente por debajo si escribimos: $$ 1 = f'(x) - \frac{1}{2}\cos(x) $$ y usando la desigualdad triangular: $$ 1 \leq |f'(x)| + \frac{1}{2}|\cos(x)| $$ luego $$ |f'(x)| \geq 1 - \frac{1}{2}|\cos(x)| \geq 1 - \frac{1}{2} = \frac{1}{2} > 0 $$ Es decir, $f'$ no se anula nunca. Como $f'(0)>0$ y $f'$ no cambia de signo, $f'$ siempre es positiva, luego $f$ es creciente en todo $\mathbb{R}$, luego $f$ tiene sólo una raíz.

Ahora bien, encontrar esa raíz de forma exacta no es posible, y nos tenemos que contentar con aproximarla.

Implementamos el método de bisección:

  • Partimos de un intervalo [xleft, xright] donde $f$ cambia de signo.
  • Elegimos la tolerancia xtol para el valor de $x$.
  • En cada iteración:
    • Evaluamos $f$ en xmed = (xleft + xright)/2.
    • Nos quedamos con un nuevo intervalo donde hay un cambio de signo y de longitug la mitad que el anterior.
    • Terminamos cuando el intervalo obtenido tenga longitud menor que xtol.
In [47]:
hold on;
plot(xs,ys)

xleft  = -4;
xright = 4;
xtol = 1e-2;
plot([xleft], [0], 'x')
plot([xright], [0], 'x')

j = 1;
sleft  = sign(f(xleft));
sright = sign(f(xright));
while xright - xleft > xtol
    xmed = (xleft + xright)/2;
    smed = sign(f(xmed));
    if smed==sleft
        xleft = xmed;
    else
        xright = xmed;
    end
    plot([xmed], [0], 'o');
    j = j + 1;
end
legends = cell(1, j+2);
legends{1} = 'f';
legends{2} = 'left end';
legends{3} = 'right end';
for k=1:j-1
    legends{k+3} = sprintf('step %d', k);
end

legend(legends)
hold off;

Ejercicio

Empaqueta el método de bisección como una función con entradas:

  • Una función f continua
  • Los extremos de un intervalo a y b.
  • Una tolerancia para el error xtol

y salidas

  • Si f(a)*f(b)>0, se rompe y no devuelve nada (imprime un mensaje de error, lanza una "excepción", etc)
  • Si f(a)*f(b)<=0, devuelve un intervalo que contiene la raíz de longitud inferior a xtol
In [ ]:
In [ ]:

Ejercicio

  • Busca ahora las raíces de
$$ f(x) = x^3 - 4x + \log(1+x^2) $$

Atención: ¿hay más de una raíz? ¿cuántas son positivas?

In [ ]:
In [ ]:
In [ ]:

Ejercicio

  • Busca una función en internet que implemente el algoritmo de bisección, compara sus argumentos, la forma de reportar el resultado y su comportamiento en caso de argumentos incorrectos con la función que has creado un poco más arriba.
  • Compara el resultado de esta función con la tuya en uno de los ejemplos anteriores.
In [ ]:
In [ ]:
In [ ]:

Resumen: raíces y puntos fijos en una variable

  • Si f(a) y f(b) tienen signos opuestos $\Rightarrow$ f tiene una raíz en el intervalo [a,b]
    • Método de bisección, que produce cotas del error cometido.
    • Tras $n$ iteraciones, la raíz está acotada en un intervalo de longitud $(b-a)/2^n$.
  • Si $g:[a,b]\rightarrow [a,b]$ $\Rightarrow$ g tiene un punto fijo en el intervalo [a,b]
  • Si $g:[a,b]\rightarrow [a,b]$ y $|g'(x)|\leq k < 1$, $\Rightarrow$
    • $g$ tiene un único punto fijo $x^*$ en el intervalo [a,b]
    • la sucesión $x_{n+1}=g(x_n)$ converge a $x^*$, para cualquier punto inicial $x_0$.
    • $|x_n - x^*| < \frac{k^n}{1-k}|x_1-x_0|$

Resumen: métodos de Newton y Secante en una variable

  • El método de Newton busca raíces de f(x), iterando $$x_{n+1}=x_{n} - \frac{f(x_n)}{f'(x_{n})}$$
    • Es una iteración de punto fijo, con $g(x) = x - f(x)/f'(x)$.
    • No es fácil saber si convergerá... pero cuando lo hace es de orden 2
    • Hay que evaluar la derivada
  • El método de la secante es similar al de Newton, salvo porque: $$x_{n+1}=x_{n} - \frac{f(x_n)\left(x_{n}-x_{n-1}\right)}{f(x_{n})-f(x_{n-1})}$$
    • No hay que evaluar la derivada
    • El orden de convergencia es algo menor

Los métodos de Newton y Secante tienen diferencias cualitativas con el método de bisección:

  • no tenemos un intervalo que acota la raíz, sino un punto que la aproxima.
  • no tenemos una cota teórica del error entre el punto $x_n$ y la raíz, a menos que se verifiquen las condiciones de convergencia y podamos acotar la segunda derivada. Tenemos dos soluciones:
    • evaluar el error en la 'y': la diferencia entre $f(x)$ y $0$.
    • estimar el error en la 'x' comparando $x_n$ con $x_{n+1}$.
  • si no conseguimos verificar las hipótesis de los teoremas de convergencia (en particular acotando $f''$), no tenemos garantía de éxito.

Como no tenemos garantía de éxito, debemos programar nuestro método de forma defensiva. Si escribimos por ejemplo:

function x=newton(f, fp, x0, ytol)
    x = x0;
    errory = abs(f(x));
    while errory > ytol
        x = x - f(x)/fp(x);
        errory = np.abs(f(x));
    end
end

nos encontraremos en un bucle infinito siempre que el método no converja.

In [ ]:

Una forma de evitar este problema es imponer un máximo de iteraciones usando un contador:

function x=newton(f, fp, x0, ytol, maxiter)
    x = x0;
    errory = abs(f(x));
    k = 0; #contador
    while errory > ytol
        x = x - f(x)/fp(x);
        errory = abs(f(x));
        k = k + 1;
    end
end

En la definición anterior hemos usado el criterio de parada $|f(x_n)|<ytol$.

El otro criterio de parada habitual es $|x_{n+1}-x_n|<xtol$.

El código que corresponde a ese criterio de parada es:

function x=newton(f, fp, x0, xtol, maxiter)
    x = x0;
    %Comenzamos con un error infinito, simplemente para que
    %la ejecución del programa no se detenga antes de empezar
    errorx = inf;
    k = 0; %contador
    #el bucle se detiene cuando se consiga la precisión deseada,
    #  pero también termina prematuramente si
    #  se alcanza el máximo de iteraciones
    while errorx > xtol && k<maxiter
        x = x - f(x)/fp(x);
        errorx = abs(f(x));
        k = k + 1;
    end
end

Finalmente, ya que no hay garantías de éxito, es interesante devolver no sólo nuestra estimación de la raíz, sino nuestra estimación de la precisión obtenida:

function [x,errorx]=newton(f, fp, x0, xtol, maxiter)
    x = x0;
    %Comenzamos con un error infinito, simplemente para que
    %la ejecución del programa no se detenga antes de empezar
    errorx = inf;
    k = 0; %contador
    #el bucle se detiene cuando se consiga la precisión deseada,
    #  pero también termina prematuramente si
    #  se alcanza el máximo de iteraciones
    while errorx > xtol && k<maxiter
        x = x - f(x)/fp(x);
        errorx = abs(f(x));
        k = k + 1;
    end
end
In [56]:
function [x,errorx]=newton(f, fp, x0, xtol, maxiter)
    x = x0;
    %Comenzamos con un error infinito, simplemente para que
    %la ejecución del programa no se detenga antes de empezar
    errorx = inf;
    k = 0; %contador
    #el bucle se detiene cuando se consiga la precisión deseada,
    #  pero también termina prematuramente si
    #  se alcanza el máximo de iteraciones
    while errorx > xtol && k<maxiter
        x = x - f(x)/fp(x);
        errorx = abs(f(x));
        k = k + 1;
    end
end
In [58]:
f = @ (x)(x^2-2);
fp = @ (x)(2*x);

% 2 iteraciones son insuficientes para conseguir la precisión deseada
[x,e] = newton(f,fp, 1, 1e-8, 2)
% 10 iteraciones sí son suficientes para conseguir la precisión deseada
% La estimación del error es mejor que la precisión mínima que le pedimos
[x,e] = newton(f,fp, 1, 1e-8, 10)
[x,e] = newton(f,fp, 1, 1e-8, 100)

x =  1.41666666666667
e =  0.00694444444444464
x =  1.41421356237469
e =    4.51061410444709e-12
x =  1.41421356237469
e =    4.51061410444709e-12

Calcular derivadas mediante cálculo simbólico

Para el método de Newton, y para muchos otros, es necesario poder evaluar la derivada de una función. Nuestras opciones son:

  • Podemos calcular la derivada a mano, aunque:
    • es fácil equivocarse si el cálculo es largo.
    • no se puede automatizar: al cambiar la función hay que derivar de nuevo.
  • Podemos usar derivación numérica, de la que hablaremos más adelante.
  • Podemos usar cálculo simbólico, que siempre consigue calcular la derivada de una función. Es fácil encontrar ejemplos en internet: ¿conseguirás ser el primero en subir un ejemplo al foro? Los profesores comentarán después tu solución...
  • Podemos buscar otro método que no necesite el cálculo de la derivada.

Ejercicio

  • Compara el método de Newton con el de bisección en los ejemplos anteriores. Constata la importancia de una buena primera aproximación. Compara el número de iteraciones necesarias para conseguir la misma precisión con uno y otro método.
In [ ]:
In [ ]:
In [ ]:

Ejercicio

Implementa el método de la secante:

  • ¿Qué argumentos necesita?
  • ¿Qué criterio de parada vas a elegir?
  • Compara el método de la secante con el de Newton para calcular la raices cuadradas. Compara el número de iteraciones necesarias para conseguir la misma precisión con uno y otro método.
  • Busca otros métodos que implementen los métodos de newton y secante en internet. Compara las decisiones tomadas y prueba uno y otro método en algún ejemplo concreto.
In [ ]:
In [ ]:
In [ ]: