2.1.1 Método de la Bisección – Ejemplo con Python

Referencia: Burden 9Ed. ejemplo 1 p50

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.

iteración 1

a = 1, b=2 c = \frac{a+b}{2} = \frac{1+2}{2} = 1.5 f(1) = (1)^3 + 4(1)^2 -10 = -5 f(1.5) = (1.5)^3 + 4(1.5)^2 -10= 2.37 f(2) = (2)^3 + 4(2)^2 -10 =14 tramo = 2-1 =1

cambio de signo a la izquierda

a = 1, b= c = 1.5

iteración 2

a = 1, b=1.5 c = \frac{1+1.5}{2} = 1.25 f(1) = -5 f(1.25) = (1.25)^3 + 4(1.25)^2 -10 = -1.794 f(1.5) = 2.37 tramo = 1.5-1 = 0.5

cambio de signo a la derecha

a = c = 1.25, b=1.5

iteración 3

continuar como tarea.

La tabla resume los valores de las iteraciones

tabla para Bisección
i a c b f(a) f(c) f(b) tramo
1 1 1.5 2 -5 2.37 14 1
2 1 1.25 1.5 -5 -1.794 2.37 0.5
3 1.25 1.5

La misma tabla se puede relizar con un algoritmo para tener los resultados más rápido y observar el comportamiento del método.

Observe los resultados de f(c), principalmente en la iteracion i=9 con tramo=0.004, respecto a la obtenida en la última iteración.

[i, a,   c,    b,    f(a),   f(c), f(b),  tramo]
 1 1.000 1.500 2.000 -5.000  2.375 14.000 1.000 
 2 1.000 1.250 1.500 -5.000 -1.797  2.375 0.500 
 3 1.250 1.375 1.500 -1.797  0.162  2.375 0.250 
 4 1.250 1.312 1.375 -1.797 -0.848  0.162 0.125 
 5 1.312 1.344 1.375 -0.848 -0.351  0.162 0.062 
 6 1.344 1.359 1.375 -0.351 -0.096  0.162 0.031 
 7 1.359 1.367 1.375 -0.096  0.032  0.162 0.016 
 8 1.359 1.363 1.367 -0.096 -0.032  0.032 0.008 
 9 1.363 1.365 1.367 -0.032  0.000  0.032 0.004 
10 1.363 1.364 1.365 -0.032 -0.016  0.000 0.002 
11 1.364 1.365 1.365 -0.016 -0.008  0.000 0.001 
raiz:  1.36474609375
>>> 

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 = tolera

import 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
while not(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

# SALIDA
print('       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 = tolera
import 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 in range(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
# GRAFICA
import 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.

Referencia: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html

2.1 Método de la Bisección – Concepto

Referencia: Burden 2.1 p48, Chapra 5.2 p124, Rodríguez 3.1 p36

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:

|p_n -p| \leq \frac{b-a}{2^{n}} 2^{-n}< 10^{-3}

usando logaritmos:

-n \log _{10}( 2) < -3 n > \frac{3}{\log _{10}( 2)} = 9.96

En consecuencia se requieren unas diez iteraciones para lograr la aproximación de 10-3. Verifique los resultados con los valores calculados.

Unidad 2 Raíces de ecuaciones en una variable

Métodos Cerrados

Método de Bisección: Concepto / / Ejemplo01

Método de la Posicíón Falsa: Concepto / / Ejemplo01

Métodos Abiertos

Método del Punto Fijo: Concepto / / Ejemplo01

Método de Newton-Raphson: Concepto / / Ejemplo01

Método de la Secante: Concepto / / Ejemplo01

Tema de introducción

Un asteroide recién descubierto pasará este jueves muy cerca de la Tierra. 23 de septiembre, 2020. https://www.eluniverso.com/noticias/2020/09/23/nota/7987777/asteroide-recien-descubierto-pasara-este-jueves-muy-cerca-tierra

3Eva_2020PAOI_T1 Distancia mínima en trayectoria

 

1.3.3 Polinomio de Taylor – Gráfica animada en Python

Referencia: Burden 7Ed, Cap 1.1 Ejemplo 3.  p11, 9Ed p11. Chapra, 4.1 p80, Taylor Series (Wikipedia)

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.

la parte adicional se muestra a partir de:

# SALIDA
# GRAFICA CON ANIMACION ------------


Algoritmo completo en Python

El tema en detalle se desarrolla en Movimiento circular – Una partícula, animación con matplotlib-Python del curso de Fundamentos de Programación

El algoritmo requiere disponer de «imagemagick», complemento que se puede descargar en:

https://imagemagick.org/script/download.php

# Aproximación Polinomio de Taylor alrededor de x0
# f(x) en forma simbólica con sympy

import numpy as np
import sympy as sym

def politaylor(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 in range(0,n,1):
    polinomio = politaylor(fx,x0,grado)
    px_tabla.append(polinomio)

# SALIDA
print('grado :  polinomio')
for grado in range(0,n,1):
    px = px_tabla[grado]
    print(str(grado)+ ' : '+str(px))
    
    # print('polinomio: ')
    # sym.pprint(px)
    # print()

# 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 in range(0,n,1):
    polinomio = px_tabla[grado]
    px = sym.utilities.lambdify(x,polinomio,'numpy')
    px_lineas[grado] = px(xi)

# SALIDA
# GRAFICA CON ANIMACION ------------
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Parametros de trama/foto
narchivo = 'Taylor01' # nombre archivo
retardo = 700   # milisegundos entre tramas
tramas = len(px_lineas)
ymax = 2*np.max(np.abs(fi))

# GRAFICA figura
figura, ejes = plt.subplots()

# Función Base
fx_linea, = ejes.plot(xi,fi,'r')

# Polinomios de tablapoli grado = 0
px_unalinea, = ejes.plot(xi, px_lineas[0],
                         '-.', label='grado: 0')

# Configura gráfica
plt.xlim([a,b])
plt.ylim([-ymax,ymax])
plt.axhline(0, color='k')  # horizontal en cero
plt.title('Polinomio Taylor: '+'f(x) = ' + str(fx))
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# cuadros de texto en gráfico
txt_x = (b+a)/2
txt_y = ymax*(1-0.1)
texto_poli = ejes.text(txt_x, txt_y*(1),
                      'p(x):',
                      horizontalalignment='center')
texto_grado = ejes.text(txt_x, txt_y*(1-0.1),
                        'grado:',
                        horizontalalignment='center')

# Nueva Trama
def unatrama(i,xi,pxi):
    
    # actualiza cada linea
    px_unalinea.set_xdata(xi)
    px_unalinea.set_ydata(pxi[i])
    etiquetap = 'p'+str(i)+'(x) = '+str(px_tabla[i])
    px_unalinea.set_label(etiquetap)
    
    # actualiza texto
    texto_poli.set_text(etiquetap)
    texto_grado.set_text('Grado: '+str(i))
    
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numerocolor = i%10
        lineacolor = 'C'+str(numerocolor)
    px_unalinea.set_color(lineacolor)
    
    return (px_unalinea, texto_poli, texto_grado)

# Limpia Trama anterior
def limpiatrama():
    
    px_unalinea.set_ydata(np.ma.array(xi, mask=True))
    px_unalinea.set_label('')
    
    texto_poli.set_text('')
    texto_grado.set_text('')
    
    return (px_unalinea,texto_poli, texto_grado)

# Trama contador
i  = np.arange(0,tramas,1)
ani = animation.FuncAnimation(figura,
                              unatrama,
                              i ,
                              fargs = (xi,px_lineas),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=True)

# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif',writer='imagemagick')

# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

1.3.2 Polinomio de Taylor – Tabla y Gráfica

Referencia: Burden 7Ed Capítulo 1.1 Ejemplo 3. p11, 9Ed p11. Chapra, 4.1 p80. Taylor Series (Wikipedia)

Continuando con el ejercicio01, 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.

f(x) = \cos (x)

alrededor de x0 = 0

grado :  polinomio
0 : 1
1 : 1
2 : -x**2/2 + 1
3 : -x**2/2 + 1
4 : x**4/24 - x**2/2 + 1
5 : x**4/24 - x**2/2 + 1
6 : -x**6/720 + x**4/24 - x**2/2 + 1
7 : -x**6/720 + x**4/24 - x**2/2 + 1
8 : x**8/40320 - x**6/720 + x**4/24 - x**2/2 + 1
9 : x**8/40320 - x**6/720 + x**4/24 - x**2/2 + 1


3 Generalizando a una función politaylor(fx,x0,n)

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 sympy

import numpy as np
import sympy as sym

def politaylor(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 in range(0,n,1):
    polinomio = politaylor(fx,x0,grado)
    px_tabla.append(polinomio)

# SALIDA
print('grado :  polinomio')
for grado in range(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.

grado :  polinomio
0 : 1
1 : 1
2 : -x**2/2 + 1
3 : -x**2/2 + 1
4 : x**4/24 - x**2/2 + 1
5 : x**4/24 - x**2/2 + 1
6 : -x**6/720 + x**4/24 - x**2/2 + 1
7 : -x**6/720 + x**4/24 - x**2/2 + 1
8 : x**8/40320 - x**6/720 + x**4/24 - x**2/2 + 1
9 : x**8/40320 - x**6/720 + x**4/24 - x**2/2 + 1

3.1 Otras formas de presentación del polinomio

Otra forma de presentar la salida es ¨pretty print¨ con sym.pprint(), borre los # de las tres últimas líneas de código para obtener:

grado :  polinomio
0 : 1
polinomio: 
1

1 : 1
polinomio: 
1

2 : -x**2/2 + 1
polinomio: 
   2    
  x     
- -- + 1
  2     

3 : -x**2/2 + 1
polinomio: 
   2    
  x     
- -- + 1
  2     

4 : x**4/24 - x**2/2 + 1
polinomio: 
 4    2    
x    x     
-- - -- + 1
24   2     

5 : x**4/24 - x**2/2 + 1
polinomio: 
 4    2    
x    x     
-- - -- + 1
24   2     

6 : -x**6/720 + x**4/24 - x**2/2 + 1
polinomio: 
    6    4    2    
   x    x    x     
- --- + -- - -- + 1
  720   24   2     


3.2 Gráfica de resultados

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 muestras en cada intervalo [a,b].

El resultado es una matriz, px_lineas, cuya fila representa el grado del polinomio, y la columna contiene los valores del polinomio de cado 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 in range(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 - GRAFICA
import matplotlib.pyplot as plt

plt.plot(xi,fi,'r',label=str(fx))

for grado in range(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()


Función Sympy para polinomio de Taylor

Existe la función desarrollada en Sympy para generar el polinomio. Se muestra un ejemplo del uso de la función como referencia.

>>> import sympy as sym
>>> x = sym.Symbol('x')
>>> fx = sym.cos(x)
>>> polinomio = sym.series(fx,x,x0=0, n=3)
>>> polinomio
1 - x**2/2 + O(x**3)
>>> 

1.3.1 Polinomio de Taylor – Error de aproximación. Ejemplo02

Referencia: Burden 7Ed cap 1.1 Ejercicio 8. Burden p15, 9Ed p15

Obtenga el tercer polinomio de Taylor P3(x) para la función:

f(x) = \sqrt{x+1}

alrededor de x0=0.

Aproxime el resultado para x=0.5, 0.75, 1.25 y 1.75 usando P3(x) y calcule los errores reales.


1. Desarrollo analítico

Se indica calcular los errores reales como la diferencia entre f(x) y el polinomio de Taylor p(x).

Siguiendo los pasos del ejercicio01, Determine el polinomio de Taylor, realice los pasos y verifique su respuesta con:

P_3(x) = 1 + \frac{1}{2}x - \frac{1}{8} x^2 +\frac{1}{16} x^3

luego calcule los valores de la tabla:

x P3(x) \sqrt{x+1} |diferencia ó error|
0.5  1.22656250000000  1.22474487139159  0.00181762860841106
0.75
1.25
1.5

Realice las observaciones a los resultados obtenidos.

Al graficar los valores de la tabla, se tiene puede observar que al alejarse x del punto de referencia x0, el error aumenta. Se representa como lo marcado en amarillo entre las curvas f(x) y el polinomio p(x).


2. Desarrollo algorítmico, instrucciones en Python

Puede reutilizar la función del polinomio de Taylor con la fórmula simbólica usada en el enlace del Ejemplo 1.4

A partir del algoritmo básico, se convierte el procedimiento a una función def-return. Con la función politaylor() se crea el polinomio y se evalúa para calcular el error respecto al valor real de la expresión.

# Aproximación Polinomio de Taylor alrededor de x0
# función en forma simbólica con sympy

import numpy as np
import sympy as sym

# Calcula n términos del polinomio de Taylor
# funcionx es simbólica
def politaylor(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.1 Ejecicio 8, Burden p15, pdf 25
# Calcule el error con polinomio Taylor grado 3

# INGRESO

# variable x es simbólica
x = sym.Symbol('x')
fx = sym.sqrt(x+1)

x0 = 0 
xi = 0.5 # donde se evalua el polinomio
n  = 3

# PROCEDIMIENTO

# Referencia, f(xi) real
fxi = fx.subs(x,xi)

# Aproximado con Taylor
polinomio = politaylor(fx,x0,n)
pxi = polinomio.subs(x,xi)

error_real = fxi - pxi

# SALIDA
print(' Taylor:     ', polinomio)
print(' xi:         ', xi)
print(' estimado  : ', pxi)
print(' real:       ', fxi)
print(' error real: ', error_real)

cuyo resultado para xi=0.5 es:

 Taylor:      x**3/16 - x**2/8 + x/2 + 1
 xi:          0.5
 estimado  :  1.22656250000000
 real:        1.22474487139159
 error real:  -0.00181762860841106

complete la tabla usando también el algoritmo en Python


2.1 Gráfica del Error entre f(x) y p(x)

Esta es una sección complementaria para realizar la gráfica mostrada en el ejemplo. Requiere el uso de la librería matplotlib. Puede revisar la sección de Recursos/Resumen Python/Gráficas 2D de línea para más detalles.

  • Para evaluar en el intervalo se requiere convertir las expresiones simbólicas a la forma numérica lambda: fxn, pxn
  • Para la gráfica, se usa el intervalo [a,b] con las muestras necesarias para una buena resolución de imagen. Se obtiene el vector xin
  • Se evalua fxn y pxn en el intervalo, obteniendo los valores en los vectores: fxni y pxni.
  • Se realiza la gráfica entre xin vs fxni y pxni
  • Para destacar el error de truncamiento, se rellena el espacio en color amarillo entre fxni y pxni, usando plt.fill_between() .
  • Para resaltar x0, se traza una línea vertical.
# cambia a forma lambda
fxn = sym.lambdify(x,fx,'numpy')
pxn = sym.lambdify(x,polinomio,'numpy')

# intervalo usando xi como referencia
a = x0        # izquierda
b = x0 + 3*xi # derecha
muestras = 51

# evaluar en intervalo
xin = np.linspace(a,b,muestras)
fxni = fxn(xin)
pxni = pxn(xin)

# Gráfica
plt.plot(xin,fxni,label='f(x)')
plt.plot(xin,pxni,label='p(x)')

plt.fill_between(xin,pxni,fxni,color='yellow')
plt.axvline(x0,color='green')

plt.title('Polinomio Taylor: f(x) vs p(x)')
plt.legend()
plt.xlabel('xi')
plt.show()

1.3 Polinomio de Taylor – Ejemplo01 con Sympy-Python

Referencia: Burden 7Ed Capítulo 1.1 Ejemplo 3.  p11, 9Ed p11. Chapra, 4.1 p80.   Taylor Series (Wikipedia)

Para la siguiente función trigonométrica:

f(x) = \cos (x)

alrededor de  x0 = 0, encontrar:
a) el segundo polinomio de Taylor (n=2),
b) el tercer polinomio de Taylor (n=3), para aproximar cos(0.01)
c) con el resultado anterior y su término residuo aproximar

\int_{0}^{0.1} \cos(x) dx

1 Desarrollo analítico

Para el polinomio de Taylor se tiene que:

P_{n}(x) = \sum_{k=0}^{n} \frac{f^{(k)}(x_0)}{k!} (x-x_0)^k P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 + + \frac{f'''(x_0)}{3!}(x-x_0)^3 + \text{...}

Para la expresión completa, se desarrollan las derivadas y se se evalua cada expresión en x0 = 0, como se muestra a continuación

f(x) = cos(x) f(0) = 1
f'(x) = -sen(x) f'(0) = 0
f”(x) = -cos(x) f”(0) = -1
f’”(x) = sen(x)  f’”(0) = 0
f4(x) = cos(x)  f4(0) = 1

En el literal a) para n=2 y x0=0:

\cos (x) = 1 + \frac{0}{1} (x-0) + \frac{-1}{2}(x-0)^2 + +\frac{\sin(\xi(x))}{6}(x-0)^3

A la expresión se añade un término más para estimar el error, como residuo o error de truncamiento, evaluado en ξ(x).

\cos (x) = 1 - \frac{1}{2}x^2 + \frac{\sin(\xi(x))}{6}x^3

con lo que si x=0.01

\cos (0.01) = 1 - \frac{1}{2}(0.01)^2 + \frac{1}{6}(0.01)^3 \sin(\xi(x)) = 0.99995 + 0.16 \text{x} 10^{-6} \sin(\xi(x))

El término del error es es del orden 10-6, la aproximación coincide por lo menos con los cinco primeros dígitos.

El residuo o error de truncamiento ξ(x) está entre 0 y x,

0<ξ(x) <0.01

Observe que los términos impares evaluados en x0=0 se anulan, por lo que el polinomio solo cambia con términos pares.

Tarea: revisar y continuar con los siguientes literales.


2 Desarrollo Algorítmico con Sympy-Python

Una forma de obtener los polinomios de Taylor es crear una función que resuelva el polinomio. Para el algoritmo, usar la forma simbólica es una opción para crear la expresión.

Por facilidad se usan funciones matemáticas expresadas de forma simbólica con Sympy, con lo que se obtiene las derivadas y se crea el polinomio para el grado requerido.

De ser necesario, revisar los conceptos para Sympy en:

Fórmulas y funciones simbólicas con Python – Sympy

2.1 Algoritmo para construir el polinomio de Taylor de grado n

El algoritmo usa la forma simbólica de la expresión para crear el polinomio.

El procedimiento consiste en crear cada término k-ésimo y añadirlo a la expresión del polinomio. Al final se presenta solo la expresión del polinomio

# Aproximación Polinomio de Taylor alrededor de x0
# f(x) en forma simbólica con sympy
# Burden 7Ed Capítulo 1.1 Ejemplo 3.p11,pdf21;9Ed p11.

import numpy as np
import sympy as sym

# INGRESO
x  = sym.Symbol('x')
fx = sym.cos(x) 
muestras = 51
x0 = 0
grado = 2       # grado>0
n  = grado + 1  # Términos de polinomio

# PROCEDIMIENTO

k = 0 # contador de términos
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

# SALIDA
print(polinomio)

un ejemplo de ejecución del algoritmo con n=3:

1 - x**2/2

Para una interpretación gráfica del resultado, luego el polinomio se evalúa en el intervalo [a, b] que incluya x0.

En las siguientes páginas 1.4.2 y 1.4.3 se propone como tema complementario, no obligatorio para la parte algorítmica, se generaliza el algoritmo usando el formato def-return y se añade la gráfica que incluye el polinomio con diferentes grados.

1.2.1 Aproximación numérica – Raíces en intervalo

[Ejemplo 1. Raíces en intervalo ]  [ Ejemplo 2. Raíces en intervalo ] ..


Ejemplo 1 . Raíces en intervalo

Referencia:  Burden 7Ed Capítulo 1.1 Ejemplo 2 p11; Burden 9Ed p8

Para la expresión mostrada encuentre una solución en el intervalo [0,1],

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

considere nombrar a la parte izquierda de la ecuación como f(x), para encontrar la solucion como f(x) = 0.

f(x) = x^5 -2x^3 + 3x^2 -1

dado que al evaluar en los extremos del intervalo [0,1]:

f(0) = (0)^5 -2(0)3 + 3(0)^2 -1 = -1 f(1) = 1^5 -2(1)^3 + 3(1)^2 -1 = 1

Los valores de la función en los extremos son -1 y 1, existe cambio de signo, y dado que f es contínua, por el teorema del valor intermedio existe un valor de x en el intervalo, tal que se satisface que el valor de la función es cero.

La gráfica de la ecuación muestra el punto o raíz a buscar:

Se pueden dividir en muchas muestras el intervalo x=[a,b] y buscar los puntos xi donde la función f(xi) cambia de signo.

Usando un algoritmo con muestras = 1001 se encuentra que existen dos puntos donde se encuentra la raíz.

raiz entre posiciones i:  [ 479.  480.]
entre los valores: 
 [ xi, fi]
[[ 0.7185     -0.00162876]
 [ 0.72        0.00219576]]

Algoritmo en Python

En el intervalo [a,b] se crean nuevos puntos de muestras para realizar la gráfica. Las muestras se usan para buscar un nuevo intervalo entre x[i] y x[i+1] donde ocurre un cambio de signo en f(x[i]).

# Burden Capítulo 1.1 Ejemplo 2 p11, pdf 21
# raices del polinomio en [a,b]

import numpy as np
import matplotlib.pyplot as plt

funcionx = lambda x: x**5 - 2*(x**3) + 3*(x**2) -1

# INGRESO
a = 0
b = 1.5
muestras = 1001

# PROCEDIMIENTO
xi = np.linspace(a,b,muestras)
fi = funcionx(xi)

# Busca cambios de signo
donde=[] # donde[i,xi,fi]
for i in range(0,muestras-1,1):
    antes = fi[i]
    despues = fi[i+1]
    signo = (np.sign(antes))*(np.sign(despues))
    if (signo<0):
        donde.append([i,xi[i],antes])
        donde.append([i+1,xi[i+1],despues])
        
donde = np.array(donde)
  
# SALIDA
print('raiz entre posiciones i: ', donde[:,0])
print('entre los valores: ')
print(' [ xi, fi]')
print(donde[:,1:])
      
# GRAFICA
plt.plot(xi,fi)
plt.xlabel('x')
plt.ylabel('fx')
plt.axhline(y=0, color='g')
plt.show()

Al observar los resultados, realice sus comentarios y recomendaciones relacionadas con la respuesta obtenida para mejorar la respuesta

Usando Scipy.Optimize

Continuando con el uso de funciones de scipy se obtiene una de las raíces, empezando la búsqueda desde 0.4.

raiz en :  [ 0.71913933]
>>>

para la otra raíz usar un nuevo punto de partida. Compare respuestas con el método anterior.

las instrucciones usadas son:

# Burden Capítulo 1.1 Ejemplo 2 p11, pdf 21
# raices del polinomio en [a,b]

import numpy as np
import scipy.optimize as opt

fx = lambda x: x**5 - 2*(x**3) + 3*(x**2) -1

# INGRESO
a  = 0
b  = 1.5
x0 = 0.4

# PROCEDIMIENTO
# fx pasa por cero, cerca de x0
donde = opt.fsolve(fx,x0)
  
# SALIDA
print('raiz en : ', donde)

[Ejemplo 1. Raíces en intervalo ]  [ Ejemplo 2. Raíces en intervalo ] ..


Ejemplo 2 . Raíces en intervalo

Referencia: Burden 7Ed cap1.1 Ejercicio 1.  p15, Burden 9Ed p15

Demuestre que las siguientes ecuaciones tienen al menos una solución en los intervalos dados

x \cos (x) - 2 x^{2} + 3 x -1 = 0

en el intervalo [0.2, 0.3] y [1.2, 1.3]

2.1 Desarrollo analítico

Del «teorema del valor intermedio«, si hay cambio de signo al evaluar la función en los puntos x=0.2 y x=0.3, debe existir un punto c donde se cumple la expresión.

f(x) = x \cos (x) - 2 x^{2} + 3 x -1 f(0.2) = 0.2 \cos (0.2) - 2 (0.2)^2 +3(0.2) -1 = -0.2839 f(0.3) = 0.2 \cos (0.3) - 2 (0.3)^2 +3(0.3) -1 = 0.00660094

Hay cambio de signo de la función en el intervalo, por lo que la ecuación debe pasar por cero, y se cumple la igualdad.

2.2 Desarrollo numérico con Python

Por simplicidad se usa la ventana iterativa. Se evalúa la función en los puntos extremos del intervalo y con los resultados se continúa de forma semejante a la sección de desarrollo analítico.

>>> import numpy as np
>>> fx = lambda x: x*np.cos(x) - 2*(x**2) + 3*x -1
>>> fx(0.2)
-0.28398668443175157
>>> fx(0.3)
0.0066009467376817454

es decir, por cambio de signo, debe haber un cruce por cero de la función en el intervalo.

2.3 Desarrollo con gráfica

Como existen varios intervalos, [0.2, 0.3] y [1.2, 1.3] se unifican los intervalos entre los extremos x=[0.2, 1.3].

Para la gráfica se crean 100 tramos (pasos o divisiones) en el intervalo que equivalen a 101 muestras en el intervalo

>>> muestras=101
>>> xi = np.linspace(0.2,1.3,muestras)
>>> fi = fx(xi)
>>> xi
array([ 0.2 , 0.211, 0.222, 0.233, 0.244, 0.255, 0.266, 0.277, ... 1.256, 1.267, 1.278, 1.289, 1.3 ])
>>> fi
array([-0.28398668, -0.24972157, -0.21601609, -0.18287411, -0.15029943, ...,  -0.13225152])

Una gráfica permite observar mejor la función en el intervalo.
Se necesita llamar a la librería matplotlib.pyplot que se resume como plt.

>>> import matplotlib.pyplot as plt
>>> plt.plot(xi,fi)
[<matplotlib.lines.Line2D object at 0x0000020C67820E80>]
>>> plt.axhline(0,color='g')
<matplotlib.lines.Line2D object at 0x0000020C678204E0>
>>> plt.show()

plt.plot() crea la gráfica con los vectores xi y fi , añadiendo como referencia en este caso una linea horizontal que pasa por cero, usando plt.axhline(). Finalmente se muestra la gráfica con plt.show().

De la gráfica, fácilmente se puede observar que existen dos puntos «c» que cumplen con la igualdad y que se encuentran en los intervalos de evaluación, con lo que se comprueba que existe solución en los intervalos presentados en el problema.

Resumen de instrucciones Python

# Raices en intervalo - Ejemplo02
import numpy as np
import matplotlib.pyplot as plt

fx = lambda x: x*np.cos(x) - 2*(x**2) + 3*x -1

# INGRESO
a = 0.2
b = 1.3
muestras = 101

# PROCEDIMIENTO
xi = np.linspace(a,b,muestras)
fi = fx(xi)

# SALIDA - GRAFICA
plt.plot(xi,fi)
plt.axhline(0,color='g')
plt.show()

Tarea: continúe con el ejercicio usando la función de scipy.optimize.fsolve() y compare resultados.

[Ejemplo 1. Raíces en intervalo ]  [ Ejemplo 2. Raíces en intervalo ]


3. Tarea

Continuar con el ejercicio observando que si el dominio es [-2,2] se tiene que:

gráfica obtenida con:

# Burden Capítulo 1.1 Ejemplo 2 p11, pdf 21
# raices del polinomio en [a,b]

import numpy as np
import matplotlib.pyplot as plt

funcionx = lambda x: x**5 - 2*(x**3) + 3*(x**2) -1

# INGRESO
a = -2
b = 2
muestras = 1001

# PROCEDIMIENTO
xi = np.linspace(a,b,muestras)
yi = funcionx(xi)

# SALIDA
plt.plot(xi,yi)
plt.axhline(0, color='green')
plt.axvline(0, color='green')
plt.show()

[Ejemplo 1. Raíces en intervalo ]  [ Ejemplo 2. Raíces en intervalo ]

1.2 Aproximación numérica – Máximo en intervalo

[Ejemplo 1. Máximo en intervalo ]  [ Ejemplo 2. Máximo en intervalo ] ..


Ejemplo 1 . Máximo en intervalo

Referencia: Burden 7Ed capítulo 1.1-ejemplo 1 p6; Burden 9Ed  p5

Determine el valor máximo de |f(x)|  en los intervalos: [1, 2] y [0.5, 1]. Siendo la función:

f(x) = 5 \cos(2x) - 2x \sin(2x)

Se puede usar dos opciones para el desarrollo: la analítica y la numérica.

1.1 Solución analítica

Se determina la derivada de f(x) y se determina el valor de x cuando f'(x) toma el valor de cero.

f(x) = 5 \cos(2x) - 2x \sin(2x) f'(x) = 5 (- 2 \sin(2x)) - [2x (2 \cos(2x)) + 2 \sin(2x) ] f'(x) = - 12 \sin(2x) - 4x \cos(2x)

f'(x) en el rango [1,2] toma el valor de cero en:

0 = - 12 \sin(2x) - 4x \cos(2x)

Situación que require un poco de trabajo adicional para encontrar el punto buscado…

1.2 Solución numérica

Otra forma es determinar el valor usando un método numérico, cuya precisión dependerá de la cantidad de muestras discreta, o tamaño de paso, que se utilicen para la evaluación.

Una gráfica permite estimar las intersecciones con los ejes y extremos de las funciones.

el máximo se encuentra en: 1.358
con el valor f(x) de: 5.67530054527

El valor máximo de |fx| en magnitud se cumple cuando la derivada es cero en un punto del intervalo.

Algoritmo en Python

  • Para observar la función, se realiza la gráfica en el rango [0.5, 2].
  • El algoritmo base corresponde al usado para una gráfica 2D, si no dispones de información previa, consulte el enlace: Gráficas 2D de línea
  • La función fx se escribe en formato lambda por simplicidad. Si no tiene  información previa sobre funciones numéricas en formato lambda revise el enlace: Funciones def-return vs lambda.
  • La precisión a usar es de mil tramos, o mil uno muestras en el intervalo [a,b], que es (2-0.5)/1000 = 0.0015‬
  • Se usa el algoritmo de búsqueda de posición del valor mayor en la función valor absoluto «fxabs».

Las instrucciones usadas en Python son:

# Burden capítulo 1.1-ejemplo 1 p6, pdf16
# Determine el maximo entre [a,b] para fx
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: 5*np.cos(2*x)-2*x*np.sin(2*x)
a = 0.5
b = 2
muestras = 1001

# PROCEDIMIENTO
xi = np.linspace(a,b,muestras)
fi = fx(xi)

fiabs = np.abs(fi)
donde = np.argmax(fiabs)

# SALIDA
print('el máximo se encuentra en: ', xi[donde])
print('con el valor f(x): ', fiabs[donde])

# GRAFICA
plt.plot(xi,fi, label='f(x)')
plt.plot(xi,fiabs, label='|f(x)|')
plt.axhline(y=0, color='g')
plt.xlabel('x')
plt.ylabel('fx')
plt.legend()
plt.title('f(x) y |f(x)|')
plt.show()

1.3 Usando Scipy.Optimize

La libreria Scipy dispone de varios algoritmos de optimización que se desarrollarán durante el curso. La comprensión de cada uno de ellos permite una aplicación efectiva de los algoritmos para obtener el resultado buscado.

Por ejemplo, usando la derivada de la función y un punto de partida x0 donde se supone, intuye o cercano donde se pretende obtener, se busca cuándo su valor es mínimo con la instrucción fsolve() se obtiene:

[ 1.35822987]
>>>

 

xalor inicial de x0

las instrucciones del algoritmo son:

# Burden capítulo 1.1-ejemplo 1 p6, pdf16
# Determine el maximo entre [a,b] para fx
# Encontrar el máximo cuando f'(x) pasa por cero

import numpy as np
import scipy.optimize as opt

# INGRESO
fx = lambda x: 5*np.cos(2*x)-2*x*np.sin(2*x)
dfx = lambda x: -12*np.sin(2*x)-4*x*np.cos(2*x)
a = 0.5
b = 2
muestras = 1001
x0 = 1 # punto inicial de búsqueda

# PROCEDIMIENTO
dondemax  = opt.fsolve(dfx,x0)

# SALIDA
print(dondemax)

compare con los resultados anteriores.

[Ejemplo 1. Máximo en intervalo ]  [ Ejemplo 2. Máximo en intervalo ] ..


Ejemplo 2. Máximo en un intervalo

Referencia: Burden 7Ed Capítulo 1.1 Ejercicio 3a p15, Burden 9Ed p15

Demuestre que f'(x) se anula al menos una vez en el  intervalo [0,1].

f(x) = 1 - e^{x} + (e-1)sen \Big( \frac{\pi}{2}x \Big)

2.1 Desarrollo analítico

Se usa el «teorema de Rolle«, si los extremos del intervalo son iguales, existe un punto intermedio c en el que la derivada es cero, en donde la función tiene un máximo.

f(0) = 1 - e^{0} + (e-1)sen(\frac{\pi}{2}0) = = 1 - 1 + (e-1)(0) = 0 f(1) = 1 - e^{1} + (e-1)sen(\frac{\pi}{2}1) = = 1 - e + (e-1)(1) = 0 f(0) = f(1)

por el teorema, debe existir un máximo, o existe un c tal que f'(c) = 0.

2.2 Desarrollo numérico y gráfico

Para encontrar el máximo, se evalúa en los extremos, se aplica Rolle y como comprobación se muestra la gráfica.

Puntos en extremos de intervalo
(xi,fi)
0 0.0
1 0.0

Algoritmo en Python

Semejante al ejercicio anterior, el punto de partida es el algoritmo para gráficas 2D.

Se plantea la función en formato lambda, usando el intervalo con 50 tramos o

# Burden Capítulo 1.1 Ejercicio 3a p15, pdf 25
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: 1-np.exp(x)+(np.exp(1)-1)*np.sin((np.pi/2)*x)
a = 0
b = 1
muestras = 51

# PROCEDIMIENTO
fa = fx(a)
fb = fx(b)

xi = np.linspace(a,b,muestras)
fi = fx(xi)

# SALIDA
print('Puntos en extremos de intervalo')
print('[xi,fi]')
print(a,fa)
print(b,fb)

# GRAFICA
plt.plot(xi,fi)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.axhline(0,color='g')
plt.show()

añada las instrucciones para encontrar el punto donde f'(x) pasa por cero, que es donde existe el máximo. use como referencia el ejemplo 1.

[Ejemplo 1. Máximo en intervalo ]  [ Ejemplo 2. Máximo en intervalo ]

1.1.1 Error-por tipos en computadoras

Referencia: Chapra 3.3 p56, Burden definición 1.15 p20

1. Error absoluto

Es la magnitud (sin signo) entre el valor «conocido» real X y el valor «estimado» Xk.

E = |X-X_k|

El valor dependerá de la magnitud de X, por ejemplo:

– Al contar monedas de 1 centavo, una persona cuenta Xk = 98 y una máquina contadora de monedas determina que el valor X = 100, el error absoluto es de 2 centavos ó 0.02 dólares.

¿Que pasaría si el conteo fuese con monedas de 1 dólar y se mantienen las mismas cantidades de monedas?


2. Error relativo

Continuando el tema del ejemplo anterior, se puede mejorar  dimensionando proporcionalmente los errores, es decir ponderarlos respecto a la magnitud usada.

e = \frac{|X-X_k|}{X}

Ponderar el error, calculando el error relativo para ambos ejemplos anteriores, se hacen comparables cuando las monedas son de un centavo o un dolar:

e = \frac{|100-98|}{100} = 0.02

3. Error de redondeo

Aparece cuando se usa una calculadora o computadora para los cálculos con números reales. La calculadora usa una cantidad finita de dígitos.
Por ejemplo:

  • el número π tiene un número infinito de dígitos,
  • si el número resultante de \sqrt{3} se eleva al cuadrado, se debería obtener 3

Sin embargo podemos comprobar que lo enunciado no se cumple al usar el computador,  así obtenemos el error de redondeo.

Usando Python se obtiene:

>>> import numpy as np
>>> numeropi=np.pi
>>> numeropi
3.141592653589793

>>> b=np.sqrt(3)
>>> b
1.7320508075688772
>>> b**2
2.9999999999999996
>>> 

4. Error de redondeo absoluto

si X_k es una aproximación de X, el error absoluto es

E = |X-X_k|

Este error se enfoca solo en la magnitud de las diferencias, no importa el signo.

5. Error de redondeo relativo

El error relativo es más significativo al usar la proporción del error en lugar del tamaño del valor.

e = \frac{|X-X_k|}{X}

6. Error de truncamiento

Resultan al usar una aproximación en lugar de un procedimiento matemático exacto. Es la diferencia entre una respuesta esperada y el valor calculado con una fórmula iterativa.

Ejemplo, al usar un polinomio de Taylor en lugar de la función original.