El método del punto fijo es un método abierto, también llamado de iteración de un punto o sustitución sucesiva, que reordena la ecuación
f(x)=0
de la forma en que x esté del lado izquierdo de la ecuación, para buscar la intersección entre la recta identidad y la curva g(x),
y = x
x=g(x)
Observe que la raíz de f(x) se encuentra en el mismo valor de x donde ocurre la intersección entre la recta identidad en color verde y la función g(x) en color naranja. Se usa la linea vertical en color morado en x=raíz como referencia de lo indicado.
El método consiste en establecer un punto inicia x0 para la búsqueda, que se usa para calcular el valor g(x0).
En la siguiente iteración el nuevo valor para x es g(x0), que se refleja en la recta identidad y nuevamente se usa para calcular g(x).
El resultado iterativo se muestra en la figura animada, donde se observa que el resultado es convergente.
La ecuación mostrada tiene una raíz en [1,2], ya que f(1)=-5 y f(2)=14.
Muestre los resultados parciales del algoritmo de Newton-Raphson con una tolerancia de 0.0001
Observe que el error representado por el tramo se va reduciendo entre cada iteración. Se debe repetir las iteraciones hasta que el error sea menor al valor tolerado.
El método de Newton-Raphson se implementa como algoritmo básico en Python
De la sección anterior, luego de realizar tres iteraciones para la tabla, notamos la necesidad de usar un algoritmo para que realice los cálculos repetitivos y muestre la tabla o directamente el resultado.
Se deduce a partir de la interpretación gráfica o por medio del uso de la serie de Taylor.
De la gráfica, se usa el triángulo formado por la recta tangente que pasa por f(xi), con pendiente f'(xi) y el eje x.
f'(x_i) = \frac{f(x_i) - 0}{x_i - x_{i+1}}
El punto xi+1 es la intersección de la recta tangente con el eje x, que es más cercano a la raíz de f(x), valor que es usado para la próxima iteración.
Reordenando la ecuación de determina la fórmula para el siguiente punto:
x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)}
El error se determina como la diferencia entre los valores sucesivos encontrados |xi+1 – xi|
La gráfica animada muestra el proceso aplicado varias veces sobre f(x) para encontrar la raíz.
Tarea
Use la serie de Taylor hasta la primera derivada para encontrar el siguiente punto de aproximación xi+1
La ecuación mostrada tiene una raíz en [1,2], ya que f(1)=-5 y f(2)=14. Muestre los resultados parciales del algoritmo de la posición falsa con una tolerancia de 0.0001
# Algoritmo Posicion Falsa para raices# busca en intervalo [a,b]import numpy as np
# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.001
# PROCEDIMIENTO
tramo = abs(b-a)
whilenot(tramo<=tolera):
fa = fx(a)
fb = fx(b)
c = b - fb*(a-b)/(fa-fb)
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if (cambia > 0):
tramo = abs(c-a)
a = c
else:
tramo = abs(b-c)
b = c
raiz = c
# SALIDAprint(raiz)
Algoritmo aumentado para mostrar la tabla de cálculos
# Algoritmo Posicion Falsa para raices# busca en intervalo [a,b]# tolera = errorimport numpy as np
# INGRESO
fx = lambda x: x**3 + 4*(x**2) -10
a = 1
b = 2
tolera = 0.0001
# PROCEDIMIENTO
tabla = []
tramo = abs(b-a)
fa = fx(a)
fb = fx(b)
whilenot(tramo<=tolera):
c = b - fb*(a-b)/(fa-fb)
fc = fx(c)
unafila = [a,c,b,fa,fc,fb,tramo]
cambio = np.sign(fa)*np.sign(fc)
if cambio>0:
tramo = abs(c-a)
a = c
fa = fc
else:
tramo = abs(b-c)
b = c
fb = fc
unafila[6]=tramo
tabla.append(unafila)
tabla = np.array(tabla)
ntabla = len(tabla)
# SALIDA
np.set_printoptions(precision=6)
print('método de la Posición Falsa ')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print(' ','tramo')
for i inrange(0,ntabla,1):
print(i, tabla[i,0:3],tabla[i,3:6])
print(' ', tabla[i,6])
print('raiz: ',c)
Observe el número de iteraciones realizadas, hasta presentar el valor de la raíz en 1.3652 con un error de 0.000079585 en la última fila de la tabla. Sin embargo, observe que la tabla solo muestra cálculos de filas completas, el último valor de c y error no se ingresó a la tabla, que se muestra como c y tramo, y es el más actualizado en los cálculos.
# Algoritmo de falsa posicion para raices# Los valores de [a,b] son seleccionados# desde la gráfica de la función# error = toleraimport numpy as np
defposicionfalsa(fx,a,b,tolera,iteramax = 20,
vertabla=False, precision=6):
'''fx en forma numérica lambda
Los valores de [a,b] son seleccionados
desde la gráfica de la función
error = tolera
'''
fa = fx(a)
fb = fx(b)
tramo = np.abs(b-a)
itera = 0
cambia = np.sign(fa)*np.sign(fb)
if cambia<0: # existe cambio de signo f(a) vs f(b)if vertabla==True:
print('método de la Posición Falsa ')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print(' ','tramo')
np.set_printoptions(precision)
while (tramo >= tolera and itera<=iteramax):
c = b - fb*(a-b)/(fa-fb)
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if vertabla==True:
print(itera,np.array([a,c,b]),
np.array([fa,fc,fb]))
if (cambia > 0):
tramo = np.abs(c-a)
a = c
fa = fc
else:
tramo = np.abs(b-c)
b = c
fb = fc
if vertabla==True:
print(' ',tramo)
itera = itera + 1
respuesta = c
# Valida respuestaif (itera>=iteramax):
respuesta = np.nan
else:
print(' No existe cambio de signo entre f(a) y f(b)')
print(' f(a) =',fa,', f(b) =',fb)
respuesta=np.nan
return(respuesta)
# PROGRAMA ----------------------# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.0001
# PROCEDIMIENTO
respuesta = posicionfalsa(fx,a,b,tolera,vertabla=True)
# SALIDAprint('raíz en: ', respuesta)
La gráfica se puede obtener añadiendo las siguientes instrucciones:
# GRAFICAimport matplotlib.pyplot as plt
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()
El método de la posición falsa, falsa posición, regla falsa o regula falsi considera dividir el intervalo cerrado [a,b] donde se encontraría una raíz de la función f(x) basado en la cercanía a cero que tenga f(a) o f(b).
El método une f(a) con f(b) con una línea recta, la intersección de la recta con el eje x representaría una mejor aproximación hacia la raíz.
Al reemplazar la curva de f(x) por una línea recta, se genera el nombre de «posición falsa» de la raíz. El método también se conoce como interpolación lineal.
A partir de la gráfica, usando triángulos semejantes, considerando que f(a) es negativo en el ejemplo, se estima que:
Calculado el valor de c, éste reemplaza a uno de los valores iniciales [a,b], cuyo valor evaluado tenga el mismo signo que f(c)
Nota: La forma de la expresión presentada para c, se usa para comparar con el método de la secante. Se obtiene sumando y restando b y reagrupando.
Control de iteraciones
Las correcciones del intervalo que se realizan en cada iteración tienen a ser más pequeñas, por lo que el control de iteraciones se realizan sobre la porción o tramo que se redujo el intervalo.
Si la reducción del intervalo es por la izquierda, tramo = c – a
Si la reducción del intervalo es por la derecha, tramo = b – c
La ecuación mostrada tiene una raíz en [1,2], ya que f(1)=-5 y f(2)=14 y existe cambio de signo. Muestre los resultados parciales del algoritmo de la bisección con una tolerancia de 0.0001
El desarrollo del ejercicio tradicionalmente realizado con lápiz, papel y calculadora, muestra el orden y detalle de las operaciones que se pueden traducir a un algoritmo en Python. El objetivo además de desarrollar la comprensión del método, permite en una evaluación observar si el estudiante conoce el método y usa apropiadamente los valores en cada iteración.
Se realiza la gráfica los puntos [c,f(c)] de la tabla para observar el resultado, resaltando que los puntos al final se aglomeran alrededor de la solución o raíz de la ecuación.
El video presenta el desarrollo básico conceptual del algoritmo en Python para una comprensión del proceso paso a paso.
Instrucciones en Python del Algoritmo básico del video
# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = toleraimport numpy as np
import matplotlib.pyplot as plt
# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.001
# PROCEDIMIENTO
tramo = b-a
whilenot(tramo<tolera):
c = (a+b)/2
fa = fx(a)
fb = fx(b)
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if cambia < 0:
a = a
b = c
if cambia > 0:
a = c
b = b
tramo = b-a
# SALIDAprint(' raiz en: ', c)
print('error en tramo: ', tramo)
El algoritmo presentado en el video se puede mejorar, por ejemplo simplificando las dos condicionales en uno.
Considere que en cada iteración se evalúa la función en tres puntos y se puede optimizar sustituyendo los valores de los extremos y solo evaluando el centro.
Finalmente se puede convertir el procedimiento en una función de Python.
Instrucciones en Python
# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = toleraimport numpy as np
defbiseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
'''
Algoritmo de Bisección
Los valores de [a,b] son seleccionados
desde la gráfica de la función
error = tolera
'''
fa = fx(a)
fb = fx(b)
tramo = np.abs(b-a)
itera = 0
cambia = np.sign(fa)*np.sign(fb)
if cambia<0: # existe cambio de signo f(a) vs f(b)if vertabla==True:
print('método de Bisección')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print(' ','tramo')
np.set_printoptions(precision)
while (tramo>=tolera and itera<=iteramax):
c = (a+b)/2
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if vertabla==True:
print(itera,[a,c,b],np.array([fa,fc,fb]))
if (cambia<0):
b = c
fb = fc
else:
a = c
fa = fc
tramo = np.abs(b-a)
if vertabla==True:
print(' ',tramo)
itera = itera + 1
respuesta = c
# Valida respuestaif (itera>=iteramax):
respuesta = np.nan
else:
print(' No existe cambio de signo entre f(a) y f(b)')
print(' f(a) =',fa,', f(b) =',fb)
respuesta=np.nan
return(respuesta)
# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.001
# PROCEDIMIENTO
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDAprint('raíz en: ', respuesta)
La gráfica se puede obtener añadiendo las siguientes instrucciones:
# GRAFICAimport matplotlib.pyplot as plt
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()
Scipy.optimize.bisect
El método de la bisección se encuentra también implementado en las librería Scipy, que también puede ser usado de la forma:
>>> import scipy.optimize as opt
>>> fx = lambda x: x**3 + 4*x**2 - 10
>>> opt.bisect(fx,1,2,xtol=0.001)
1.3642578125
que es el valor de la variable ‘a’ de la tabla para la última iteración del ejercicio. Lo que muestra que el algoritmo realizado tiene un valor más aproximado.
Sin embargo por didáctica y mejor comprensión de los métodos y su implementación en algoritmos que es parte del objetivo de aprendizaje, se continuará desarrollando la forma básica en Python.
El método más simple para buscar raíces de una ecuación, se basa en el teorema del valor intermedio, búsqueda binaria, partición de intervalos o de Bolzano.
En el intervalo donde existe un cruce por cero de la función f(x), el algoritmo busca la raíz al reducir el intervalo en la mitad (bisección), seleccionando el sub- intervalo donde se mantenga el cambio de signo de la función f(x).
Los pasos a seguir son los siguientes:
el intervalo [a,b] se divide siempre en la mitad c.
Si la función f(x) cambia de signo sobre un intervalo, se evalúa el valor de la función en el punto medio f(c).
La posición de la raíz se determina en el punto medio del sub-intervalo, izquierdo o derecho, dentro del cual ocurre un «cambio de signo».
el proceso se repite hasta obtener una mejor aproximación
La gráfica muestra una animación del proceso, observe la forma en que progresivamente se acercan los puntos [a,b], donde se mantienen valores con signo diferente entre f(a) y f(b).
Para describir mejor el método, observamos la gráfica en una sola iteración.
Para la primera iteración se tiene que la función tiene un cambio de signo dentro del intervalo [a,b].
El intervalo se divide en la mitad, representado por el punto c, obteniendo el sub-intervalo izquierdo [a,c] o sub-intervalo derecho [c,b].
El sub-intervalo que contiene la función con un cambio de signo, se convierte en el nuevo intervalo por analizar en la siguiente iteración.
El error del método de la bisección se estima como el ancho o tamaño del intervalo [a,b] de la última iteración realizada. Si el error es menor que la tolerancia del ejercicio, el algoritmo se detiene y se considera encontrada la raíz.
Suponga que f ∈ C[a,b] y f(a)*f(b)<0, f es una función en el intervalo [a,b] y que presenta un cambio de signo.
|p_n - p| \leq \frac{b-a}{2^n} \text{donde } n \geq 1
la desigualdad implica que pn converge a p con una razón de convergencia de orden:
O \Big(\frac{1}{2^n}\Big)
es decir:
p_n =p+O \Big( \frac{1}{2^n} \Big)
Con lo que se puede determinar el número de iteraciones necesarias para encontrar la raíz, tal como se muestra en el siguiente ejercicio.
El ejercicio se presenta como un complemento para la sección 1.4 que permite obtener una gráfica animada.
Esta sección es complementaria y usada solo como referencia para exponer el tema. Normalmente se da una explicación breve en el laboratorio de computadoras.
Ejercicio 1. Taylor para cos(x), tabla y gráfica en Python
Referencia: Burden 7Ed Capítulo 1.1 Ejemplo 3. p11, 10Ed p8. Chapra, 4.1 p80. Taylor Series (Wikipedia)
Continuando con el Ejemplo01, se generaliza el algoritmo para crear una tabla de polinomios de Taylor de diferente grado.
Se complementa el ejercicio con el gráfico de cada polinomio para interpretar los resultados, alrededor de x0 = 0
Sección complementaria, no obligatoria para la parte algorítmica
En el ejercicio presentado requiere resolver con varios grados de polinomio, por lo que se generaliza convirtiendo el procedimiento del algoritmo anterior al formato de función def-return. Cada polinomio intermedio se añade a una tabla de resultados:
# Aproximación Polinomio de Taylor alrededor de x0# f(x) en forma simbólica con sympyimport numpy as np
import sympy as sym
defpolitaylor(fx,x0,n):
k = 0
polinomio = 0
while (k <= n):
derivada = fx.diff(x,k)
derivadax0 = derivada.subs(x,x0)
divisor = np.math.factorial(k)
terminok = (derivadax0/divisor)*(x-x0)**k
polinomio = polinomio + terminok
k = k + 1
return(polinomio)
# PROGRAMA -------------# Capitulo 1 Ejemplo 2, Burden p11, pdf 21# INGRESO
x = sym.Symbol('x')
fx = sym.cos(x)
x0 = 0
n = 10 # Grado polinomio Taylor
a = -5 # intervalo [a,b]
b = 5
muestras = 51
# PROCEDIMIENTO# tabla polinomios
px_tabla = []
for grado inrange(0,n,1):
polinomio = politaylor(fx,x0,grado)
px_tabla.append(polinomio)
# SALIDAprint('grado : polinomio')
for grado inrange(0,n,1):
px = px_tabla[grado]
print(str(grado)+ ' : '+str(px))
# print('polinomio: ')# sym.pprint(px)# print()
Con lo que se obtiene los polinomios para cada grado calculado.
La forma gráfica de cada polinomio se obtiene evaluando cada polinomio para obtener las líneas en el intervalo [a, b] para cada punto del vector xi .
Se utiliza un cierto número de muestrasen cada intervalo [a,b].
El resultado es una matriz, px_lineas, cuya fila representa el gradodel polinomio, y la columna contiene los valores del polinomio de cada grado evaluado en cada punto xi
# GRAFICA - TABLA polinomios ------
xi = np.linspace(a,b,muestras)
# Forma lambda, simplifica evaluación
fxn = sym.utilities.lambdify(x,fx,'numpy')
fi = fxn(xi)
# lineas de cada grado de polinomio
px_lineas = np.zeros(shape=(n,muestras), dtype =float)
for grado inrange(0,n,1):
polinomio = px_tabla[grado]
px = sym.utilities.lambdify(x,polinomio,'numpy')
px_lineas[grado] = px(xi)
los valores se pueden graficar añadiendo las siguientes instrucciones:
# SALIDA - GRAFICAimport matplotlib.pyplot as plt
plt.plot(xi,fi,'r',label=str(fx))
for grado inrange(0,n,2):
etiqueta = 'grado: '+str(grado)
plt.plot(xi, px_lineas[grado],
'-.',label = etiqueta)
ymax = 2*np.max(fi)
plt.xlim([a,b])
plt.ylim([-ymax,ymax])
plt.xlabel('x')
plt.ylabel('y')
plt.title('Aproximación con Polinomios de Taylor')
plt.legend()
plt.show()
