Un posible inconveniente en el método de Newton Raphson es implementar la evaluación de la derivada. La derivada también se puede aproximar mediante una diferencia finita dividida hacia atrás:
f'(x_i) = \frac{f(x_{i-1})-f(x_i)}{x_{i-1}-x_i}
la que se sustituye en la ecuación del método de Newton-Raphson para obtener:
Encontrar la solución a la ecuación, usando el método del punto fijo
f(x):e^{-x} - x = 0
Desarrollo Analítico
Al igual que los métodos anteriores, es conveniente determinar el intervalo [a,b] donde es posible evaluar f(x). Se revisa si hay cambio de signo en el intervalo para buscar una raíz.
Para el ejemplo, el rango de observación será [0,1], pues f(0)=1 es positivo y f(1)=-0.63 es negativo.
Para el punto fijo, se reordena la ecuación para para tener una ecuación con la variable independiente separada.
Se obtiene por un lado la recta identidad y=x, por otro se tiene la función g(x).
Se buscará la intersección entre las dos expresiones .
x = e^{-x} g(x) = e^{-x}
Se puede iniciar la búsqueda por uno de los extremos del rango [a,b]. Por ejemplo:
iniciando desde el extremo izquierdo x=a,
se determina el valor de b = g(x) ,
se determina la diferencia de la aproximación o error = |b-a|, tolera=0.001
se proyecta en la recta identidad, como el nuevo punto de evaluación
x=b.
se repite el proceso para el nuevo valor de x, hasta que el error sea menor al tolerado.
En caso que el proceso no converge, se utiliza un contador de iteraciones máximo, para evitar tener un lazo infinito.
El proceso realizado en la tabla se muestra en la gráfica, para una función que converge.
Algoritmo en Python
El algoritmo en python se presenta como una función para usarla fácilmente como un bloque en otros ejercicios. Se requiere la función gx, el punto inicial y la tolerancia, la variable de número de iteraciones máxima, iteramax, permite controlar la convergéncia de la función con el método.
Se presenta el algoritmo mejorado, convirtiendo el procedimiento del método a una función Python.
# Algoritmo de punto fijo# [a,b] intervalo de búsqueda# error = toleraimport numpy as np
defpuntofijo(gx,a,tolera, iteramax = 15):
i = 1 # iteración
b = gx(a)
tramo = abs(b-a)
while(tramo>=tolera and i<=iteramax ):
a = b
b = gx(a)
tramo = abs(b-a)
i = i + 1
respuesta = b
# Validar respuestaif (i>=iteramax ):
respuesta = np.nan
return(respuesta)
# PROGRAMA ---------# INGRESO
fx = lambda x: np.exp(-x) - x
gx = lambda x: np.exp(-x)
a = 0 # intervalo
b = 1
tolera = 0.001
iteramax = 15 # itera máximo
muestras = 51 # gráfico
tramos = 50
# PROCEDIMIENTO
respuesta = puntofijo(gx,a,tolera)
# SALIDAprint(respuesta)
la respuesta obtenida del problema es
0.566908911921
Para obtener la gráfica básica se determinan los puntos para cada función fx y gx. Se añaden las siguiente instrucciones al algoritmo anterior:
# GRAFICA# calcula los puntos para fx y gx
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
yi = xi
import matplotlib.pyplot as plt
plt.plot(xi,fi, label='f(x)')
plt.plot(xi,gi, label='g(x)')
plt.plot(xi,yi, label='y=x')
if (respuesta!= np.nan):
plt.axvline(respuesta)
plt.axhline(0, color='k')
plt.title('Punto Fijo')
plt.legend()
plt.show()
Scipy.optimize.fixed_point
El método del punto fijo se encuentra implementado en Scipy, que también puede ser usado de la forma:
>>> import scipy.optimize as opt
>>> opt.fixed_point(gx,a,xtol=0.001,maxiter=15)
array(0.5671432948307147)
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 buscarla intersección entre la recta identidad y la curva g(x), como se muestran en los siguientes ejemplos.
y = x
x=g(x)
Observe que la raiz 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=raiz 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.
Ejemplo 1
f(x):e^{-x} - x = 0
se reordena para tener:
x = e^{-x} g(x) = e^{-x}
Ejemplo 2
f(x): x^2 - 2x -3 = 0
se reordena para tener:
x = \frac {x^2 - 3}{2} g(x) = \frac {x^2 - 3}{2}
Ejemplo 3
f(x): \sin (x) = 0
puede ser complicado despejar x, por lo que se simplifica el proceso sumando x en ambos lados.
x = \sin (x) + x g(x) = \sin (x) + x
El método proporciona una fórmula para predecir un valor nuevo de x en función del valor anterior:
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
f(x) = x^3 + 4x^2 -10 =0
Desarrollo Analítico
El método requiere obtener la derivada f'(x) de la ecuación para el factor del denominador.
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.
Las demás iteraciones se dejan como tarea
Algoritmo con Python
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 raiz.
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 raiz 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
f(x) = x^3 + 4x^2 -10 =0
Desarrollo Analítico
Semejante a los métodos anteriores, el método posición falsa, falsa posición, regla falsa o regula falsi, usa un intervalo [a,b] para buscar la raiz.
Se divide el intervalo en dos partes al calcular el punto c que divide al intervalo siguiendo la ecuación:
el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo
tramo = |c-a| = | 1.2631 - 1| = 0.2631 a = c = 1.2631
iteración 2
a = 1.2631 , b = 2 f(1.2631) = -1.6031 f(2) = 14 c = 2 - 14 \frac{1.2631-2}{-1.6031-14} = 1.3388 f(1.3388) = (1.3388)^3 + 4(1.3388)^2 -10 = -0.4308
el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo
tramo = |c-a| = |1.3388 - 1.2631| = 0.0757 a = c = 1.3388
iteración 3
a = 1.3388 , b = 2 f(1.3388) = -0.4308 f(2) = 14 c = 2 - 14 \frac{1.3388-2}{-0.4308-14} = 1.3585 f(1.3585) = (1.3585)^3 + 4(1.3585)^2 -10 = -0.1107
el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo
tramo = |c-a| = |1.3585 - 1.3388| = 0.0197 a = c = 1.3585
valores que se resumen en la tabla
Método de posición falsa
a
c
b
f(a)
f(c)
f(b)
tramo
1
1.2631
2
-5
-1.6031
14
0.2631
1.2631
1.3388
2
-1.6031
-0.4308
14
0.0757
1.3388
1.3585
2
-0.4308
-0.1107
14
0.0197
1.3585
…
2
se puede contrinuar con las iteraciones como tarea
Algoritmo en Python
Algoritmo básico del video
# 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)
tabla.append([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
tabla = np.array(tabla)
ntabla = len(tabla)
# SALIDA
np.set_printoptions(precision=4)
for i inrange(0,ntabla,1):
print('iteración: ',i)
print('[a,c,b]: ', tabla[i,0:3])
print('[fa,fc,fb]: ', tabla[i,3:6])
print('[tramo]: ', tabla[i,6])
print('raiz: ',c)
print('error: ',tramo)
Observe el número de iteraciones realizadas, hasta presentar el valor de la raiz en 1.3652 con un error de 0.0003166860575976038 en la útlima 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.
Con el método de la posición falsa, determine la cantidad de iteraciones necesarias para resolver:
f(x) = \cos (x) - x
Realice una tabla que muestre: [ i, a, b, c, f(c) ].
Seleccióne el rango [a,b] usando una gráfica y compruebe los resultados.
Compare con otros métodos para encontrar raíces.
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 raiz.
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 redución del intervalo es por la izquierda, tramo = c – a
Si la redución del intervalo es por la derecha, tramo = b – c
La ecuación mostrada tiene una raiz 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
f(x) = x^3 + 4x^2 -10 =0
Desarrollo Analítico
Como parte del desarrollo del ejercicio se presenta las iteraciones para el algoritmo, tradicionalmente realizadas con una calculadora.
Se grafican los puntos [c,f(c)] de la tabla para obsevar el resultado de forma gráfica, resaltando que los puntos al final se aglomeran alrededor de la solución o raiz de la ecuación.
Escriba sus observaciones y preguntas sobre los resultados.
Algoritmo en Python
El video presenta el desarrollo básico conceptual del algoritmo en Python para una comprensión del proceso de creación algoritmica.
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)
Mejorando el algoritmo
El algoritmo presentado en el video se puede mejorar, por ejemplo simplificando los dos condicionales en uno.
Considere que siempre se evalúa la función en tres puntos y se puede optimizar sustituyendo los valores de los extremos y solo evaluando el centro.
El resultado se mejorar usando los valores finales del último intervalo [a,b] y obteniendo una nueva midad ‘c’ al final.
Para mostrar la tabla de la sección anterior se agregan los resultados parciales a una tabla que permita mostrar el resultado al final
# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = toleraimport numpy as np
# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 2
tolera = 0.001
# PROCEDIMIENTO
tabla = []
tramo = b-a
fa = fx(a)
fb = fx(b)
i = 1
while (tramo>tolera):
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
i = i + 1
cambia = np.sign(fa)*np.sign(fc)
if (cambia<0):
b = c
fb = fc
else:
a=c
fa = fc
tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)
raiz = c
# SALIDA
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), tramo]')
# print(tabla)# Tabla con formato
n=len(tabla)
for i inrange(0,n,1):
unafila = tabla[i]
formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
unafila = formato.format(*unafila)
print(unafila)
print('raiz: ',raiz)
El ultimo complemento al algoritmo consiste en realizar la gráfica, seleccionando solo las columnas correspondientes a [c,f(c)].
Se ordenan los datos en forma ascendente antes de graficarlos usando solo sus índices con np.argsort(xi)
# Algoritmo de Bisección# GRAFICAimport matplotlib.pyplot as plt
xi = tabla[:,2]
yi = tabla[:,5]
# ordena los puntos para la grafica
orden = np.argsort(xi)
xi = xi[orden]
yi = yi[orden]
plt.plot(xi,yi)
plt.plot(xi,yi,'o')
plt.axhline(0, color="black")
plt.xlabel('x')
plt.ylabel('y')
plt.title('Bisección en f(x)')
plt.grid()
plt.show()
Scipy.optimize.bisect
El método de la bisección se encuentra también implementado en las libreria Scipy, que también puede ser usado de la forma:
>>> import scipy.optimize as opt
>>> 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 se basa en el teorema del valor intermedio, conocido como método de la bisección, búsqueda binaria, partición de intervalos o de Bolzano.
Es un tipo de búsqueda incremental en el que:
el intervalo se divide siempre en la mitad.
Si la función cambia de signo sobre un intervalo, se evalúa el valor de la función en el punto medio.
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 el proceso en forma animada, 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).
Observamos la gráfica para una sola iteración y asi describir mejor el método.
Para la primera iteración se tiene como procedimiento que la función tiene un cambio de signo en el intervalo [a,b].
En 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 a ser analizado en la siguiente iteración
Cota de Error
Referencia: Teorema 2.1 Burden 9Ed. p,51.
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)
cantidad de iteraciones
Referencia: Ejemplo 2 Burden 9Ed. p.52/pdf.62
Determine la cantidad de iteraciones necesarias para resolver
f(x) = x^3 + 4x^2 -10 =0
con exactitud de 10 – 3 en el intervalo [1,2].
Desarrollo: Se busca encontrar un entero n que satisface la ecuación: