2.5.1 Método de la Secante – Ejemplo con Python

[ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


1. Ejercicio

Referencia: Burden 2.1 ejemplo 1 p38

La ecuación mostrada tiene una raíz en el intervalo [1,2], ya que f(1) = -5 y f(2) = 14
Muestre los resultados parciales del algoritmo de la secante con una tolerancia de 0.0001

f(x) = x^3 + 4x^2 -10 =0

Secante Metodo animado

[ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [función]
..


2. Desarrollo Analítico

siendo:

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

Se probará con valores iniciales que tienen el mismo signo f(x[i-1]), f(x[i]), que no se considera en el método cerrado de la posición falsa

itera = 0

X-1 = a =2.5,  X0 =b=4, tolera = 0.001

f(2.5) = (2.5)^3 + 4(2.5)^2 -10 = 30.625 f(4) = (4)^3 + 4(4)^2 -10 =118 x_{1}= 4 - 118\frac{(2.5 - 4)}{(30.625 - 118)} = 1.9742 errado = |1.9742 - 4| = 2.02575

itera = 1

X0 = 4,  X1 = 1.9742

f(4) =118 f(1.9742) = ( 1.9742)^3 + 4( 1.9742)^2 -10 =13.2855 x_{2}= 1.9742 - (13.2855)\frac{(4 - 1.9742)}{(118 -13.2855)} =1.7172 errado = |1.7172 - 1.9742| = 0.2570

itera = 2

X1 =1.9742, X2 =1.7172

f(1.9742) = 13.2855 f(1.7172) = (1.7172)^3 + 4(1.7172)^2 -10 =6.8594 x_{3}= 1.7172 - (6.8594)\frac{(1.9742 - 1.7172)}{(6.8594 -13.2855)} = 1.4428 errado = | 1.4428- 1.7172| =0.2743

desarrollando una tabla con el algoritmo se tiene:

método de la Secante
i [ x[i-1], xi, x[i+1], f[i-1], fi ] tramo
0 [  2.5    4.       1.974249  30.625   118.      ] 2.0257510729613735
1 [  4.     1.974249 1.717233 118.       13.285584] 0.2570160557153698
2 [1.974249 1.717233 1.442883  13.285584  6.859484] 0.27434949602777015
3 [1.717233 1.442883 1.376796   6.859484  1.331607] 0.06608786561380797
4 [1.442883 1.376796 1.365656   1.331607  0.19207 ] 0.011139180030541596
5 [1.376796 1.365656 1.365232   0.19207   0.007041] 0.00042390961991034537
raiz en : 1.3652324200312267

[ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


3. Algoritmo en Python

El algoritmo básico a usar es:

# Método de secante
import numpy as np

# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 4
tolera = 0.001
iteramax = 20
precision = 5

c = (a+b)/2 # probar dos puntos en un mismo lado de f(x)
# PROCEDIMIENTO
xi_1 = c
xi = b
itera = 0
tramo = np.abs(xi-xi_1)
print('i','[ x[i-1], xi, x[i+1], f[i-1], fi ]','tramo')
np.set_printoptions(precision)
while not(tramo<tolera or itera>iteramax):
    fi_1 = fx(xi_1)
    fi = fx(xi)
    xi1 = xi-fi*(xi_1 - xi)/(fi_1-fi)
    tramo = np.abs(xi1-xi)
    print(itera,np.array([xi_1,xi, xi1, fi_1,fi]),tramo)
    xi_1 = xi
    xi = xi1
    itera = itera + 1

respuesta = xi

# SALIDA
print('raiz: ',respuesta)

 

[ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


4. Función en Python

Una función resumida, se recomienda realizar la validación de cambio de signo entre [a,b]

# Método de secante
import numpy as np

def secante_raiz(fx,a,b,tolera, iteramax=20,
                 vertabla=True, precision=6):
    '''fx en forma numérica lambda
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    '''
    xi_1 = a
    xi = b
    itera = 0
    tramo = np.abs(xi-xi_1)
    if vertabla==True:
        print('método de la Secante')
        print('i','[ x[i-1], xi, x[i+1], f[i-1], fi ]','tramo')
        np.set_printoptions(precision)
    while not(tramo<tolera or itera>iteramax):
        fi_1 = fx(xi_1)
        fi = fx(xi)
        xi1 = xi-fi*(xi_1 - xi)/(fi_1-fi)
        tramo = np.abs(xi1-xi)
        if vertabla==True:
            print(itera,np.array([xi_1,xi,xi1,fi_1,fi]),tramo)
        xi_1 = xi
        xi = xi1
        itera = itera + 1
    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')
    respuesta = xi
    
    return(respuesta)

# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10
a = 1
b = 4
tolera = 0.001
c = (a+b)/2 # probar dos puntos en un mismo lado de f(x)

# PROCEDIMIENTO
respuesta = secante_raiz(fx,c,b,tolera,
                         vertabla=True,precision = 5)
# SALIDA
print('raíz en: ', respuesta)

Scipy.optimize.newton – Secante

El método de la secante se encuentra implementado en Scipy en la forma de algoritmo de newton, que al no proporcionar la función para la derivada de f(x), usa el método de la secante:

>>> import scipy.optimize as opt
>>> opt.newton(fx,xa, tol=tolera)
1.3652320383201266

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


[ Secante ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

2.4.1 Método del Punto fijo – Ejemplo con Python

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


1. Ejercicio

Referencia: Burden 2.2 p41, Chapra 6.1 p143, Rodríguez 3.2 p44

Encontrar la solución a la ecuación, usando el método del punto fijo, considerando tolerancia de 0.001

f(x): e^{-x} - x = 0

Punto Fijo 01a_GIF

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

..


2. 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).

f(x): e^{-x} - x = 0 x = e^{-x}

g(x) = e^{-x}punto fijo 01b

Se buscará la intersección entre las dos expresiones .

Se puede iniciar la búsqueda por uno de los extremos del rango [a,b].

iteración 1

– iniciando desde el extremo izquierdo x=a,

x_0 = 0

– se determina el valor de b = g(x) ,

x_1 = g(0) = e^{-0} = 1

– se determina la diferencia de la aproximación o error = |b-a|

tramo =|1-0|= 1

– 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.

iteración 2

x_1 = 1 x_2 = g(1) = e^{-1} = 0.3678 tramo =|0.3678 - 1|= 0.6322

iteración 3

x_2 = 0.3678 x_3 = e^{-0.3678} = 0.6922 tramo =|0.6922-0.3678|= 0.3244

La tabla resume los valores de las iteraciones

Método del punto fijo
iteración x xnuevo = g(x) |error|
1 0 1.0 1
2 1.0 0.3678 0.6322
3 0.3678 0.6922 0.3244
4 0.6922
5

El proceso realizado en la tabla se muestra en la gráfica, para una función que converge.

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


3. Algoritmo en Python

El algoritmo en Python del video se adjunta para seguimiento. 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 convergencia de la función con el método.

El resultado del algoritmo se presenta como:

raíz en:  0.5669089119214953
errado :  0.0006477254067881466
itera  :  14

Instrucciones en Python

# Algoritmo de punto fijo
# [a,b] son seleccionados desde la gráfica
# error = tolera

import numpy as np

# INGRESO
fx = lambda x: np.exp(-x) - x
gx = lambda x: np.exp(-x)

a = 0
b = 1
tolera = 0.001
iteramax = 15

# PROCEDIMIENTO
i = 1 # iteración
b = gx(a)
tramo = abs(b-a)
while not(tramo<=tolera or i>iteramax):
    a = b
    b = gx(a)
    tramo = abs(b-a)
    i = i + 1
respuesta = b
# Valida convergencia
if (i>=iteramax):
    respuesta = np.nan

# SALIDA
print('raíz en: ', respuesta)
print('errado : ', tramo)
print('itera  : ', i)

Para obtener la gráfica básica se determinan los puntos para cada función fx y gx. Como los valore de a y b cambiaron durante el algoritmo, para la siguiente sección es necesario volver a escribirlos u obtener una copia de los valores en otra variable para reutilizarlos en ésta sección.

Se añaden las siguiente instrucciones al algoritmo anterior:

# GRAFICA
import matplotlib.pyplot as plt
# calcula los puntos para fx y gx
a = 0
b = 1
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
yi = xi

plt.plot(xi,fi, label='f(x)',
         linestyle='dashed')
plt.plot(xi,gi, label='g(x)')
plt.plot(xi,yi, label='y=x')

plt.axvline(respuesta, color='magenta',
            linestyle='dotted')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Punto Fijo')
plt.grid()
plt.legend()
plt.show()

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


4. Función en Python

El resultado del algoritmo como una función y mostrando la tabla de resultados por iteraciones es:

método del Punto Fijo
i ['xi', 'gi', 'tramo']
0 [0. 1. 1.]
1 [1.       0.367879 0.632121]
2 [0.367879 0.692201 0.324321]
3 [0.692201 0.500474 0.191727]
4 [0.500474 0.606244 0.10577 ]
5 [0.606244 0.545396 0.060848]
6 [0.545396 0.579612 0.034217]
7 [0.579612 0.560115 0.019497]
8 [0.560115 0.571143 0.011028]
9 [0.571143 0.564879 0.006264]
10 [0.564879 0.568429 0.003549]
11 [0.568429 0.566415 0.002014]
12 [0.566415 0.567557 0.001142]
13 [0.567557 0.566909 0.000648]
raíz en:  0.5669089119214953

Instrucciones en Python

# Algoritmo de punto fijo
# [a,b] son seleccionados desde
# la gráfica de la función
# error = tolera

import numpy as np

def puntofijo(gx,a,tolera, iteramax=20,vertabla=True, precision=6):
    """
    g(x) se obtiene al despejar una x de f(x)
    máximo de iteraciones predeterminado: iteramax
    si no converge hasta iteramax iteraciones
    la respuesta es NaN (Not a Number)
    """
    itera = 0
    b = gx(a)
    tramo = abs(b-a)
    if vertabla==True:
        print('método del Punto Fijo')
        print('i', ['xi','gi','tramo'])
        np.set_printoptions(precision)
        print(itera,np.array([a,b,tramo]))
    while(tramo>=tolera and itera<=iteramax):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        itera = itera + 1
        if vertabla==True:
            print(itera,np.array([a,b,tramo]))
    respuesta = b
    # Valida respuesta
    if itera>=iteramax:
        respuesta = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')
    return(respuesta)

# PROGRAMA ----------------------
# INGRESO
fx = lambda x: np.exp(-x) - x
gx = lambda x: np.exp(-x)

a = 0
b = 1
tolera = 0.001
iteramax = 15

# PROCEDIMIENTO
respuesta = puntofijo(gx,a,tolera,iteramax,vertabla=True)

# SALIDA
print('raíz en: ', respuesta)

De añadirse la parte gráfica, el bloque no requiere volver a escribir los valores de a y b, al ser usados como variable internas a la función y se mantienen en la parte principal del programa.

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 valor predeterminado de iteraciones si no se escribe es 500 y la tolerancia predeterminada es xtol=1e-08
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fixed_point.html


Tarea

  • Revisar lo que sucede cuando el valor inicial a esta a la derecha de la raíz.
  • Validar que la función converge, revisando que |g'(x)|<1, o que la función g(x) tenga pendiente menor a pendiente de la recta identidad.

[ Punto fijo ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

2.3.1 Método de Newton-Raphson – Ejemplo con Python

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


1. Ejercicio

ReferenciaBurden 2.1 ejemplo 1 p38

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


[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


2. 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

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


3. 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)

La gráfica presentada para revisar f(x) se realiza con las instrucciones:

# GRAFICA
import matplotlib.pyplot as plt
a = 1
b = 4
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()

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


4. Función en Python

Se convierte el algoritmo a una función, con partes para ver la tabla, y se obtienen los siguientes resultados:

i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 2.  14.  28.   1.5  0.5]
1 [ 1.5     2.375  18.75    1.3733  0.1267]
2 [1.3733e+00 1.3435e-01 1.6645e+01 1.3653e+00 8.0713e-03]
3 [1.3653e+00 5.2846e-04 1.6514e+01 1.3652e+00 3.2001e-05]
raíz en:  1.3652300139161466

Instrucciones en Python

# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx,dfx,xi, tolera, iteramax=100,
                   vertabla=False, precision=4):
    '''fx y dfx en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('método de Newton-Raphson')
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1

    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')

    return(xi)

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

x0 = 2
tolera = 0.001

# PROCEDIMIENTO
respuesta = newton_raphson(fx,dfx,x0, tolera, vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

La gráfica se la puede añadir usando las instrucciones dadas en el algoritmo básico realizado al inicio par ala comprensión del método.

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 gráfica básica 2D de la función f(x)

[ Newton-Raphson ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

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

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


1. Ejercicio

Referencia: Burden 2.1 ejemplo 1 p38

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

f(x) = x^3 + 4x^2 -10 =0

posicionfalsa01_GIF

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


2. 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 raíz.

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 continuar con las iteraciones como tarea

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


3. 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)
    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 in range(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.

método de la Posición Falsa 
i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [1.       1.263158 2. ] [-5.       -1.602274 14. ]
   0.26315789473684204
1 [1.263158 1.338828 2. ] [-1.602274 -0.430365 14. ]
   0.0756699440909967
2 [1.338828 1.358546 2. ] [-0.430365 -0.110009 14. ]
   0.019718502996940224
3 [1.358546 1.363547 2. ] [-0.110009 -0.027762 14. ]
   0.005001098217311428
4 [1.363547 1.364807 2. ] [-2.776209e-02 -6.983415e-03  1.400000e+01]
   0.0012595917846898175
5 [1.364807 1.365124 2. ] [-6.983415e-03 -1.755209e-03  1.400000e+01]
   0.0003166860575976038
6 [1.365124 1.365203 2. ] [-1.755209e-03 -4.410630e-04  1.400000e+01]
    7.958577822231305e-05
raíz en:  1.3652033036626001

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


4. Función en Python

# Algoritmo de falsa posicion para raices
# Los valores de [a,b] son seleccionados
# desde la gráfica de la función
# error = tolera

import numpy as np

def posicionfalsa(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 respuesta
        if (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)
# SALIDA
print('raíz en: ', respuesta)

La gráfica se puede obtener añadiendo las siguientes instrucciones:

# GRAFICA
import 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()

[ Falsa Posición ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]

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

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ función ]
..


1. Ejercicio

Referencia: Burden 2.1 ejemplo 1 p38

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

f(x) = x^3 + 4x^2 -10 =0

Método de la Bisección animado

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [función]
..


2. Desarrollo Analítico

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.

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

cambio de signo a la izquierda

a = 1, b= c = 1.5 tramo = |1.5-1| =0.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

cambio de signo a la derecha

a = c = 1.25, b=1.5 tramo = |1.5-1.25| = 0.25

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 0.5
2 1 1.25 1.5 -5 -1.79 2.37 0.25
3 1.25 1.5

La misma tabla se puede realizar 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 iteración i=9 con tramo=0.00097 que representa el error de estimación del valor vs tolerancia.

i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [1, 1.5, 2] [-5.     2.375 14.   ]
   0.5
1 [1, 1.25, 1.5] [-5.     -1.7969  2.375 ]
   0.25
2 [1.25, 1.375, 1.5] [-1.7969  0.1621  2.375 ]
   0.125
3 [1.25, 1.3125, 1.375] [-1.7969 -0.8484  0.1621]
   0.0625
4 [1.3125, 1.34375, 1.375] [-0.8484 -0.351   0.1621]
   0.03125
5 [1.34375, 1.359375, 1.375] [-0.351  -0.0964  0.1621]
   0.015625
6 [1.359375, 1.3671875, 1.375] [-0.0964  0.0324  0.1621]
   0.0078125
7 [1.359375, 1.36328125, 1.3671875] [-0.0964 -0.0321  0.0324]
   0.00390625
8 [1.36328125, 1.365234375, 1.3671875] [-3.2150e-02  7.2025e-05  3.2356e-02]
   0.001953125
9 [1.36328125, 1.3642578125, 1.365234375] [-3.2150e-02 -1.6047e-02  7.2025e-05]
   0.0009765625
raíz en:  1.3642578125
>>> 

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.

método de la bisección

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [función]
..


3. Algoritmo en Python

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

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [función]

..


4. Función en Python

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 = tolera
import numpy as np

def biseccion(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 respuesta
        if (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)
# SALIDA
print('raíz en: ', respuesta)

La gráfica se puede obtener añadiendo las siguientes instrucciones:

# GRAFICA
import 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.

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

[ Bisección ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [función]

1.4.2 Polinomio de Taylor – Gráfica animada en Python

Referencia: Burden 7Ed, Cap 1.1 Ejemplo 3.  p11, 10Ed p8. 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.4.1 Polinomio de Taylor – Función, Tabla y Gráfica con Python

Polinomio de Taylor: [ Ejercicio ] [ función politaylor() ] [ pprint() ]
[ gráfica ] [ función en Sympy ]
..


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

f(x) = \cos (x)
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

Polinomio de Taylor: [ Ejercicio ] [ función politaylor() ] [ pprint() ]
[ gráfica ] [ función en Sympy ]
..


2. Función politaylor(fx,x0,n) en Python

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

Polinomio de Taylor: [ Ejercicio ] [ función politaylor() ] [ pprint() ]
[ gráfica ] [ función en Sympy ]
..


3. 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     

Polinomio de Taylor: [ Ejercicio ] [ función politaylor() ] [ pprint() ]
[ gráfica ] [ función en Sympy ]
..


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

Polinomio de Taylor: [ Ejercicio ] [ función politaylor() ] [ pprint() ]
[ gráfica ] [ función en Sympy ]
..


5. Función de 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)
>>> 

Polinomio de Taylor: [ Ejercicio ] [ función politaylor() ] [ pprint() ]
[ gráfica ] [ función en Sympy ]

1.4 Polinomio de Taylor – Ejemplos con Sympy-Python

Ejemplo Taylor: [ cos(x) ] [ analítico ] [ algoritmo ] [ raíz() ] [ función ]
..


Ejemplo 1. Polinomio de Taylor para cos(x)

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

Para la siguiente función trigonométrica f(x), alrededor de x0=0, encontrar :

f(x) = \cos (x)

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

Ejemplo Taylor: [ cos(x) ] [ analítico ] [ algoritmo ] [ raíz() ] [ función ]
..


1.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 evalúa 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.

Ejemplo Taylor: [ cos(x) ] [ analítico ] [ algoritmo ] [ raíz() ] [ función ]
..


1.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:

Sympy – Fórmulas y funciones simbólicas

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.

Ejemplo Taylor: [ cos(x) ] [ analítico ] [ algoritmo ] [ raíz() ] [ función ]
..


Ejemplo 2. Polinomio de Taylor y el error de aproximación

Referencia: Burden 7Ed cap 1.1 Ejercicio 8. Burden 10Ed p8

Obtenga el tercer polinomio de Taylor P3(x) para la función f(x), alrededor de x0=0.

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

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


2.1 Desarrollo analítico

Las instrucciones piden calcular los errores reales como la diferencia entre f(x) y el polinomio de Taylor p(x).

Usando los pasos del ejemplo01, se determina el polinomio de Taylor con n=3. Verifique su respuesta con el polinomio mostrado y calcula los valores de la tabla:

P_3(x) = 1 + \frac{1}{2}x - \frac{1}{8} x^2 +\frac{1}{16} x^3
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 puede observar que al alejarse x del punto de referencia x0, el error aumenta. El error se marca en amarillo entre las curvas f(x) y el polinomio p(x).


Ejemplo Taylor: [ cos(x) ] [ analítico ] [ algoritmo ] [ raíz() ] [ función ]
..


2.2 Desarrollo algorítmico como función, instrucciones en Python

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

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.3 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 evalúa 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()

Ejemplo Taylor: [ cos(x) ] [ analítico ] [ algoritmo ] [ raíz() ] [ función ]

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

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


Ejemplo 1 . Raíces en intervalo

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

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 ]  [ 2. Raíces en intervalo ]
..


Ejemplo 2 . Raíces en intervalo

Referencia: Burden 7Ed cap1.1 Ejercicio 1

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 ]  [ 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 ]  [ 2. Raíces en intervalo ]

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

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


Ejemplo 1 . Máximo en intervalo

Referencia: Burden 7Ed capítulo 1.1-ejemplo 1 p6; Burden 10Ed  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 requiere 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 librería Scipy dispone de varios algoritmos de optimización que se desarrollarán durante el curso.

valor inicial de x0

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]
>>>

 

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 ]  [ 2. Máximo en intervalo ]
..


Ejemplo 2. Máximo en un intervalo

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

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 ]  [ 2. Máximo en intervalo ]