2.3.1 Método de Newton-Raphson Ejemplo con Python

Referencia: Burden ejemplo 1 p51

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.

f(x) = x^3 + 4x^2 -10 f'(x) = 3x^2 + 8x x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)}

Para el desarrollo se inicia la búsqueda desde un punto en el intervalo [1,2], por ejemplo el extremo derecho, x1=2.

iteración 1

f(2) = (2)^3 + 4(2)^2 -10 = 14 f'(2) = 3(2)^2 + 8(2) = 28 x_{2} = 2 -\frac{14}{28} = 1.5 tramo = |2 -1.5| = 0.5

iteración 2

f(1.5) = (1.5)^3 + 4(1.5)^2 -10 = 2.375 f'(1.5) = 3(1.5)^2 + 8(1.5) = 18.75 x_{3} = 1.5 -\frac{2.375}{18.75} = 1.3733 tramo = |1.5 -1.3733| = 0.1267

iteración 3

f(1.3733) = (1.3733)^3 + 4(1.3733)^2 -10 = 0.1337 f'(1.3733) = 3(1.3733)^2 + 8(1.3733) = 16.6442 x_{4} = 1.3733 -\frac{0.1337}{16.6442} =1.3652 tramo = |1.3733 -1.3652| = 0.0081

La tabla resume los valores de las iteraciones

Método de Newton-Raphson
iteración xi xnuevo tramo
1 2 1.5 0.5
2 1.5 1.3733 0.1267
3 1.3733 1.3653 0.0081
4

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.

 ['xi', 'xnuevo', 'tramo']
[[2.0000 1.5000 5.0000e-01]
 [1.5000 1.3733 1.2667e-01]
 [1.3733 1.3653 8.0713e-03]
 [1.3653 1.3652 3.2001e-05]]
raiz en:  1.3652300139161466
con error de:  3.200095847999407e-05

Al algoritmo básico se les añade lo necesario para mostrar la tabla con los valores de las iteraciones.

# Método de Newton-Raphson
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

import numpy as np

# INGRESO
fx  = lambda x: x**3 + 4*(x**2) - 10
dfx = lambda x: 3*(x**2) + 8*x

x0 = 2
tolera = 0.001

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

scipy.optimize.newton

El método de Newton-Raphson se encuentra implementado en Scipy, que también puede ser usado de la forma:

>>> import scipy.optimize as opt
>>> opt.newton(fx,x0, fprime=dfx, tol = tolera)
1.3652300139161466
>>> 

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


Tarea

Calcule la raíz de f(x) = e-x – x, empleando como valor inicial x0 = 0

  • convertir el algoritmo a una función
  • Revisar las modificaciones si se quiere usar la forma simbólica de la función.
  • incorpore la grafica básica 2D de la función f(x)

2.3 Método de Newton-Raphson – Concepto

Referencia: Burden 2.3 p66, Chapra 6.2 p148, Rodríguez 3.3 p52

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

2.2.1 Método de la Posición Falsa – Ejemplo con Python

Referencia: Burden 9Ed ejemplo 1 p51

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

posicionfalsa01_GIF


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:

c = b - f(b) \frac{a-b}{f(a)-f(b)}

iteración 1

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

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

# SALIDA
print(raiz)

Algoritmo aumentado para mostrar la tabla de cálculos

# Algoritmo Posicion Falsa para raices
# busca en intervalo [a,b]
# tolera = error

import 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)
while not(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 in range(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.

iteración:   0
[a,c,b]:     [1.     1.2632 2.    ]
[fa,fc,fb]:  [-5.     -1.6023 14.    ]
[tramo]:     1.0
iteración:   1
[a,c,b]:     [1.2632 1.3388 2.    ]
[fa,fc,fb]:  [-1.6023 -0.4304 14.    ]
[tramo]:     0.26315789473684204
iteración:   2
[a,c,b]:     [1.3388 1.3585 2.    ]
[fa,fc,fb]:  [-0.4304 -0.11   14.    ]
[tramo]:     0.0756699440909967
iteración:   3
[a,c,b]:     [1.3585 1.3635 2.    ]
[fa,fc,fb]:  [-0.11   -0.0278 14.    ]
[tramo]:     0.019718502996940224
iteración:   4
[a,c,b]:     [1.3635 1.3648 2.    ]
[fa,fc,fb]:  [-2.7762e-02 -6.9834e-03  1.4000e+01]
[tramo]:     0.005001098217311428
iteración:   5
[a,c,b]:     [1.3648 1.3651 2.    ]
[fa,fc,fb]:  [-6.9834e-03 -1.7552e-03  1.4000e+01]
[tramo]:     0.0012595917846898175
iteración:   6
[a,c,b]:     [1.3651 1.3652 2.    ]
[fa,fc,fb]:  [-1.7552e-03 -4.4106e-04  1.4000e+01]
[tramo]:     0.0003166860575976038
raiz:   1.3652033036626001
error:  7.958577822231305e-05

Ejemplo 2

Referencia: Burden ejemplo 3 p.74/pdf.84

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.

2.2 Método de la Posición Falsa – Concepto

Referencia: Burden p56, Chapra 5.3 p131

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:

\frac{f(a)}{c-a} = -\frac{f(b)}{b-c} \frac{f(a)}{c-a} = \frac{f(b)}{c-b}

que al despejar c, se obtiene:

c = b - f(b) \frac{(a-b)}{f(a)-f(b)}

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

 

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()