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
%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
.
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
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
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$.
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] $$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:
En estas condiciones:
- 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.$
%run biseccion_dibu.py
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:
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.
%run Ejercicio1.py
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.
- 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$.
%run newton_dibu.py
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$.
%run Ejercicio2.py
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.
%run secante_dibu.py
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$.
%run Ejercicio3.py
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.
f = lambda x: x**3 - 2*x**2 +1
g = lambda x: - np.sqrt((x**3+1)/2)
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
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.
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.
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:
punto_fijo
con el valor inicial estimado ($tol=10^{-12}$ y $maxiter=200$).%run Ejercicio4.py
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é?
%run Ejercicio5.py
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.
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.
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.
%run Ejercicio6.py
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.
%run Ejercicio7.py
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.
%run Ejercicio8.py
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.
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'$:
sympy
(mirar Nota al final de esta práctica).newton
de scipy.optimize
.%run Ejercicio9.py
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''$:
sympy
(mirar Nota al final de esta práctica).newton
de scipy.optimize
.%run Ejercicio10.py
Para calcular derivadas de forma simbólica podemos utilizar el módulo sympy.
import sympy as sym
Comenzamos declarando una variable simbólica y la función
x = sym.Symbol('x', real=True)
Declaramos la función simbólica y podemos calcular sus derivadas
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.
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
x = np.linspace(-3,3,100)
plt.plot(x,f(x),x,df(x),x,d2f(x))
plt.legend(['f','df','d2f'])
plt.show()