Raíces de ecuaciones no lineales

Contenidos

Raíz de una ecuación

Una ecuación escalar es una expresión de la forma

$$f(x)=0$$

donde $f:\Omega\subset\mathbb{R}\rightarrow\mathbb{R}$ es una función.

Cualquier número $\alpha$ tal que $f(\alpha)=0$ se dice que es una solución de la ecuación o una raíz de $f$.

Ejemplo 1. Calcular de forma aproximada las raíces de la ecuación:

$$x^3-2x^2+1=0$$

Insertamos las figuras en la notebook, en lugar de que cada figura abra una ventana nueva

In [4]:
%matplotlib inline

Con from __future__ import division hacemos que el resultado de la división con / dé siempre como resultado un float. Si no lo incluyéramos la división con / entre dos enteros tendría como resultado otro entero.

Importamos los módulos matplotlib.pyplot y numpy.

In [2]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
In [3]:
x = np.linspace(-1,2)              # definimos un vector con 50 elementos en (-1,2)
f = lambda x : x**3 - 2*x**2 + 1   # definimos la función
OX = 0*x                           # definimos un vector de ceros del tamaño de x
In [4]:
plt.plot(x,f(x))                   # dibujamos la función 
plt.plot(x,OX,'k-')                # dibujamos el eje X
plt.show()

Vemos que tiene tres raíces que son, aproximadamente, $x\simeq -0.6$, $x\simeq 0.9$ y $x\simeq 1.6$.

Separación de raíces

A partir de la gráfica anterior podemos separar las raíces en intervalos que cumplan las condiciones del Teorema de Bolzano:

Sea $f:[a,b]\rightarrow\mathbb{R}$ una función contínua con $f(a)f(b)<0$. Entonces existe al menos una raíz de $f$ en $[a,b]$.

Es decir, buscamos intervalos donde $f$ tome signos distintos en los extremos.

Y si además la función es monótona creciente o decreciente está garantizado que en el intervalo $[a,b]$ existe una única raíz.

En nuestro ejemplo, tres intervalos que cumplen estas condiciones son:

$$ [-1,-0.1]\quad [0.5,1.1] \quad [1.5,2] $$

Método de bisección

El método de bisección parte de una función $f$ contínua en un intervalo $[a,b]$ donde se cumplen las condiciones del teorema de Bolzano, es decir, la función tiene distinto signo en los extremos. Para que funcione no es necesario que las raíces estén separadas, pero sí es conveniente (si están las tres raíces en $[a,b]$ sólo aproximará una).

El método de bisección genera una sucesión de intervalos que:

  • Cumplen las condiciones del teorema de Bolzano.
  • Y por lo tanto contiene siempre al menos una raíz.
  • Cada intervalo tiene una longitud que es la mitad del anterior.

En estas condiciones:

  • La aproximación de la raíz es uno de los extremos del intervalo y por lo tanto
  • El error máximo cometido es la longitud del intervalo.

Algoritmo

  • Sea $a_1=a$, $b_1=b$.
  • Para $k=1,2,\ldots,MaxIteraciones$
    • Calcular el punto medio $x_k=\frac{a_k+b_k}{2}$ y evaluar $f(x_k)$.
    • Si $x_k$ satisface el criterio de parada, parar.
    • En el caso contrario,
      • si $f(a_k)f(x_k)<0$ entonces:

        $ a_{k+1}=a_{k},$ $ b_{k+1}=x_k,$

      • en el caso contrario:

        $ a_{k+1}=x_{k},$ $ b_{k+1}=b_k.$

In [2]:
%run biseccion_dibu.py

Criterios de parada

Como la sucesión ${x_k}$ tal que $x_k\rightarrow\alpha$, con $f(\alpha)=0$ es, en general, infinita, vamos a introducir un criterio para decidir cuando paramos.

Los criterios de parada se basan en:

  • El error absoluto $|x_k−x_{k−1}|\lt tol$.
  • El error relativo $\frac{|x_k−x_{k−1}|}{|x_k|}\lt tol$.
  • El residual $|f(x_k)| \lt tol$.

Ejercicio 1

Escribir una función llamada biseccion que tenga como argumentos de entrada la función $f$, los extremos del intervalo $a$, $b$, la cota del error absoluto $tol$ y el número máximo de iteraciones $maxiter$ y como argumentos de salida la solución aproximada y el número de iteraciones realizadas.

Usarla para aproximar las tres raíces de $f(x)= x^3 - 2 x^2 + 1$, con $tol=10^{-12}$ y $maxiter=200$.

Escribir las tres raíces y el número de iteraciones realizadas para aproximarlas.

In [5]:
%run Ejercicio1.py
-0.61803398875 40
1.0 40
1.61803398875 39

El método de Newton-Raphson

Es el método más usado cuando podemos calcular la derivada exacta. En cada paso, se sustituye la función por su recta tangente en el punto $x_k$ y se calcula su raíz, que es una aproximación de la raíz de la función.

Algoritmo

  • Sea $x_0$ un punto inicial.
  • Para $k=1,2,\ldots,MaxIter$:
    • Calcular $x_k=x_{k-1}-\frac{f(x_{k-1})}{f'(x_{k-1})}$
    • Si $x_k$ satisface el criterio de parada, parar.
    • En el caso contrario, hacer otra iteración.

En el dibujo vemos la sucesión de rectas tangentes (en rojo) para calcular la primera raíz tomando $x_0=-1$.

In [3]:
%run newton_dibu.py

Ejercicio 2

Escribir un programa llamado Newton que tenga como argumentos de entrada la función $f$, su derivada $df$, el punto inicial $x_0$, la cota de error absoluto $tol$ y el número máximo de iteraciones $maxiter$ y como argumentos de salida la solución aproximada y el número de iteraciones realizadas.

Utilizarlo para aproximar las tres raíces de la función $f(x)= x^3 - 2 x^2 + 1$, con $tol=10^{-12}$ y $maxiter=200$.

In [6]:
%run Ejercicio2.py
-0.61803398875 6
1.0 1
1.61803398875 6

Método de la secante

Es una variante del método de Netwon-Raphson. Sustituimos la tangente por una secante.

El problema del método de Newton-Raphson es que necesitamos tener la $f'(x)$.

En el método de la secante aproximamos $f'(x)$ con diferencias finitas: $$f'(x_{k-1})\approx\frac{f(x_{k-1})-f(x_{k-2})}{x_{k-1}-x_{k-2}}.$$

Necesitamos dos puntos iniciales, $x_0$ y $x_1$ para la primera iteración: $$ x_k= x_{k-1} -f( x_{k-1})\frac{x_{k-1}-x_{k-2}}{f(x_{k-1})-f(x_{k-2})}$$

En el dibujo vemos la sucesión de rectas secantes (en rojo) para calcular la primera raíz tomando $x_0=0$ y $x_1=-1$. En verde, la posición de las sucesivas raíces aproximadas.

In [4]:
%run secante_dibu.py

Ejercicio 3

Escribir una función llamada secante que tenga como argumentos de entrada la función $f$, los puntos iniciales $x_0$ y $x_1$, la cota de error absoluto $tol$ y el número máximo de iteraciones y como argumentos de salida la solución aproximada y el número de iteraciones realizadas.

Utilizarla para aproximar las tres raíces de la función $f(x)= x^3 - 2 x^2 + 1$, con $tol=10^{-12}$ y $maxiter=200$.

In [7]:
%run Ejercicio3.py
-0.61803398875 9
1.0 6
1.61803398875 8

Método de punto fijo

En el método de punto fijo sustituimos la ecuación $f(x)=0$, por la ecuación equivalente $g(x)=x$ de modo que si $\alpha$ es un punto fijo de $g$, es decir $g(\alpha)=\alpha$, entonces también será una raíz de $f$, o lo que es lo mismo $f(\alpha)=0$.

Se dice que la función $g:[a,b]\to\mathbb{R}$ tiene un punto fijo $\alpha\in [a,b]$ si $g(\alpha)=\alpha$. El método del punto fijo se basa en la iteración

$$ x_{k+1}=g(x_k),\quad k\geq 0, $$

donde $x_0$ es la aproximación inicial que debemos proporcionar.

No existe una manera única de expresar esta equivalencia. Veamos un ejemplo.

Sea $f(x)=x^3 - 2 x^2 + 1$ y buscamos las raíces de $f$ tal que $x^3 - 2 x^2 + 1 = 0$. Podemos reorganizar la ecuación como $x^3 + 1 = 2 x^2$, o lo que es lo mismo

$$ x = -\sqrt{\frac{x^3+1}{2}}$$

La solución de esta segunda ecuación, también llamada punto fijo de la función

$$ g(x) = -\sqrt{\frac{x^3+1}{2}}$$

será también una raíz de $f$.

Comprobémoslo gráficamente.

In [8]:
f = lambda x: x**3 - 2*x**2 +1
g = lambda x: - np.sqrt((x**3+1)/2)
In [9]:
a = -1.; b = 0;
x = np.linspace(a, b)              # Vector de 50 elementos equiespaciados en (a,b)
plt.plot(x, f(x),'g-',label='f')              
plt.plot(x, g(x),'r-',label='g')
plt.plot(x, 0*x,'k-')              # Eje OX
plt.plot(x, x,'b-',label='y = x')

raiz = -0.61803
plt.plot(raiz,0,'ro',label=u'raíz')     
plt.plot(raiz,raiz,'bo',label='punto fijo')

plt.legend(loc='best')
plt.show()

Y vemos que el punto donde $f(x)=0$ coincide con el punto donde $x=g(x)$

Por lo tanto, el punto fijo de $$ g(x) = -\sqrt{\frac{x^3+1}{2}}\quad (1)$$ es una raíz de $f$.

El Teorema de la aplicación contractiva dice: sea $g$ derivable definida en el intervalo $[a,b] \subset \mathbb{R} $ y $x_0 \in [a,b]$ un punto del intervalo. Supongamos que

  • $x \in [a,b] \quad \Rightarrow \quad g(x)\in [a,b]$
  • $\left|g'(x)\right|\le k\lt 1$ para todo $x\in [a,b]$

Entonces $g$ tiene un único punto fijo $\alpha \in [a,b]$, y la sucesión $x_n$ definida como $x_{i+1}=f(x_i)$ que tiene como punto inicial $x_0$ converge a $\alpha$ con orden al menos lineal.

In [10]:
a = -0.8; b = -0.3;
x = np.linspace(a, b) 

plt.plot(x, g(x),'r-', label='g')
plt.plot(x, x, 'b-',label='y = x')
plt.plot([a,b],[b,a],'k-') # Dibuja la otra diagonal

pf = -0.61803
plt.plot(pf,pf,'bo',label='punto fijo') 

plt.axis([a, b, a, b])     # Gráfica en [a,b]x[a,b]
plt.legend(loc='best')
plt.show()

En la gráfica vemos que se cumplen las condiciones del teorema de la aplicación contractiva para la función $g$ en el intervalo $[-0.8, -0.3]$. Si tomamos como punto inicial de la iteración de punto fijo un punto de este intervalo el algoritmo de punto fijo converge.


Ejercicio 4

Escribir un programa que contenga una función punto_fijo que tenga como argumentos de entrada una función de iteración $g$, el punto inicial $x_0$, la cota de error absoluto $tol$ y el número máximo de iteraciones $maxiter$ y como argumentos de salida la solución aproximada y el número de iteraciones realizadas.

Utilizarlo para aproximar la raíz de la función $f(x)= e^{-x}-x$ con la función de iteración $g(x)=e^{-x}$.

Seguir los siguientes pasos:

  • Dibujar $f$ en el intervalo $(-1,1)$. Dibujar sobre la misma gráfica $g$ y la recta $y = x$ para comprobar que el valor del punto fijo de $g$ coincide con el de la raíz de $f$.
  • Dibujar sobre la misma gráfica $g$ y la recta $y = x$ buscando un intervalo $[a,b]$ donde se cumplan las condiciones del teorema de la aplicación contractiva para estimar un punto inicial.
  • Usar la función punto_fijo con el valor inicial estimado ($tol=10^{-12}$ y $maxiter=200$).
  • Escribir la raiz y el número de iteraciones.
In [11]:
%run Ejercicio4.py
0.567143290409 46

Ejercicio 5

Sea la función $f(x)= x - \cos(x)$, con $x \in \left(0,1\right)$. Comprobar gráficamente que la ecuación $f(x)=0$ tiene la misma raíz que $g_i(x)=x$ con $i=1,2,3,4$ siendo

$$g_1(x)=\cos(x)\quad g_2(x)=\arccos(x) \quad g_3(x)=2x-\cos(x) \quad g_4(x)=x-\dfrac{x-\cos(x)}{1+\sin(x)}$$

Estudiar gráficamente qué funciones $g_i$ que pueden usarse para calcular el punto fijo por el método de iteración de punto fijo y, para estas funciones, encontrar gráficamente un intervalo donde se cumplan las condiciones del teorema de la aplicación contractiva. Utilizándolas, calcular el punto fijo de las funciones del apartado anterior con una tolerancia $tol = 10^{-12}$ y dar el número de iteraciones. Comprobar que es un punto fijo dibujándolo como un punto azul sobre la gráfica.

¿Cuál es la mejor función $g_i$ para aproximar la solución por el método de punto fijo? ¿Por qué?

In [12]:
%run Ejercicio5.py
Se pueden usar g1 y g4 porque, en un entorno del punto fijo,
su pendiente en valor absoluto es menor que 1

No se pueden usar g2 y g3 porque, en un entorno del punto fijo,
su pendiente en valor absoluto es mayor que 1
g1
0.739085133215 69

g4
0.739085133215 5

Funciones de python para calcular raíces

SciPy es una biblioteca open source de herramientas y algoritmos matemáticos para Python. SciPy contiene módulos para optimización, álgebra lineal, integración, interpolación, funciones especiales, FFT, procesamiento de señales y de imagen, resolución de ODEs y otras tareas para la ciencia e ingeniería.

Las funciónes que calculan raíces están incluídas dentro del módulo de optimización (scipy.optimize).

Si solo queremos usar las funciones de un módulo, pordemos cargar únicamente ese módulo.

In [13]:
import scipy.optimize as op

La función newton calcula las raíces utilizando el método de la secante o el de Newton-Raphson.


Ejercicio 6

Utilizando el método de newton y usando la función newton del módulo scipy.optimize calcular todas las soluciones positivas de la ecuación $\mathrm{sen}x-0.1x=0$. Empieza dibujando la función para estimar los valores iniciales a utilizar. Usar como parámetros $tol=10^{-12}$ y $maxiter=100$. Comprobar que son raíces dibujándolos como puntos rojos sobre la gráfica.

In [14]:
%run Ejercicio6.py
2.85234189445
7.0681743581
8.42320393236

Ejercicio 7

Utilizando el método de la secante y usando la función newton del módulo scipy.optimize calcular la solución positiva más cercana a cero de la ecuación $\cosh(x)\cos(x)-1=0$. Empieza dibujando la función para estimar los valores iniciales a utilizar. Usar como parámetros $tol=10^{-12}$ y $maxiter=100$. Comprobar que es una raiz dibujándola como un punto rojo sobre la gráfica.

In [15]:
%run Ejercicio7.py
4.73004074486

Ejercicio 8

Utilizando la función bisect del módulo scipy.optimize calcular $\sqrt[3]{75}$. Dibújala para estimar el intervalo inicial a utilizar. Usar como parámetros $tol=10^{-12}$ y $maxiter=100$. Comprobar que es una raiz dibujándola como un punto rojo sobre la gráfica.

In [16]:
%run Ejercicio8.py
4.21716332651

Cálculo de extremos y puntos de inflexión de una función

Si $f$ es una función derivable en un punto $x_0$ de su dominio y $x_0$ es un extremo relativo, entonces

$$ f'(x_0)=0 $$

Es decir, los extremos son raíces de la función derivada.

Si $f$ es una función derivable dos veces en un punto $x_0$ de su dominio y $x_0$ es un punto de inflexión, entonces

$$ f''(x_0)=0 $$

Es decir, los puntos de inflexión son raíces de la función derivada segunda.


Ejercicio 9

Dada la función: $$f(x)=x^3+\ln (x+7) \cos (4 x)-1$$

Calcular los extremos relativos (máximos y mínimos) de $f$, es decir,las raíces de $f'$:

  • Definir las funciones lambda de $f$, $f'$ y $f''$ utilizando el módulo sympy (mirar Nota al final de esta práctica).
  • Dibujar $f'$ con el eje OX y probar con distintos intervalos para el valor de x, hasta encontrar un intervalo que contenga todas sus raíces.
  • Obtener con un error máximo de $10^{-12}$ todos los ceros de la función $f'$ utilizando la función newton de scipy.optimize.
  • Dibujar $f$, de forma que contenga todos sus extremos relativos y, mirando la gráfica, decidir cuales de los puntos anteriores son máximos y cuales mínimos. Dibujar sobre la función los máximos relativos con puntos rojos y los mínimos con puntos verdes.
  • Comprobar que las derivadas segundas de $f$ en los máximos son negativas y en los mínimos positivas.
In [17]:
%run Ejercicio9.py
x1 =  -1.34140197854
x2 =  -0.856813731846
x3 =  0.00458534319858
x4 =  0.738973116799
d2f(x1) =  -26.0434065713
d2f(x2) =  22.3860101409
d2f(x3) =  -31.1536138346
d2f(x4) =  36.4367932146

Ejercicio 10

Dada la función: $$f(x)= e^{-x^2}\ln (x^2+1)$$

Calcular los puntos de inflexión de $f$, es decir, las raíces de $f''$:

  • Definir las funciones lambda de $f$, $f'$ y $f''$ utilizando el módulo sympy (mirar Nota al final de esta práctica).
  • Dibujar $f''$ con el eje OX y probar con distintos intervalos para el valor de x, hasta encontrar un intervalo que contenga todas sus raíces.
  • Obtener con un error máximo de $10^{-12}$ todos los ceros de la función $f''$ utilizando la función newton de scipy.optimize.
  • Dibujar $f$, de forma que contenga todos sus puntos de inflexión. Dibujar sobre la función los puntos de inflexión con puntos azules.
In [18]:
%run Ejercicio10.py
x1 =  -1.34109575434
x2 =  -0.393009307011
x3 =  0.393009307011
x4 =  1.34109575434

NOTA: Calcular la derivada de una función

Para calcular derivadas de forma simbólica podemos utilizar el módulo sympy.

In [19]:
import sympy as sym

Comenzamos declarando una variable simbólica y la función

In [20]:
x = sym.Symbol('x', real=True)

Declaramos la función simbólica y podemos calcular sus derivadas

In [21]:
f_sim   = sym.cos(x)*(x**3+1)/(x**2+1)
df_sim  = sym.diff(f_sim,x)
d2f_sim = sym.diff(df_sim,x)

Para usarlas como funciones numéricas las lambificamos.

In [22]:
f   = sym.lambdify([x], f_sim,'numpy') 
df  = sym.lambdify([x], df_sim,'numpy') 
d2f = sym.lambdify([x], d2f_sim,'numpy') 

Y ya podemos, por ejemplo, dibujarlas con plot

In [23]:
x = np.linspace(-3,3,100)
plt.plot(x,f(x),x,df(x),x,d2f(x))
plt.legend(['f','df','d2f'])
plt.show()