3.2 Pivoteo parcial por filas con Python

Referencia: Chapra 9.4.2 p268, pdf 292. Burden 9Ed 6.2 p372. Rodriguez 4.0 p105

Los métodos de solución a sistema de ecuaciones, tienen en los primeros pasos en común usar la matriz aumentada y pivoteada por filas. Para compartir estos pasos y simplificar la presentación de los métodos, se presentan como una de las primeros algoritmos a implementar.

Para mostrar todo el desarrollo del proceso se usa como referencia un ejercicio.


Ejercicio

Referencia: Rodriguez cap4.0 Ejemplo 1 pdf.105

Ejemplo 1. Un comerciante compra tres productos A, B, C, pero en las facturas únicamente consta la cantidad comprada en Kg y el valor total de la compra. Se necesita determinar el precio unitario de cada producto.  Dispone de solo tres facturas con los siguientes datos:

Ejemplo:
Cantidad Valor ($)
Factura X1 X2 X3 Pagado
1 4 2 5 60.70
2 2 5 8 92.90
3 5 4 3 56.30

Los precios unitarios se pueden representar por las variables x1, x2, x3 para escribir el sistema de ecuaciones que muestran las relaciónes de cantidad, precio y valor pagado:

\begin{cases} 4x_1+2x_2+5x_3 = 60.70 \\ 2x_1+5x_2+8x_3 = 92.90 \\ 5x_1+4x_2+3x_3 = 56.30 \end{cases}

Sistema de ecuaciones como Matriz y Vector

El sistema de ecuaciones se escribe en la forma algebraica como matrices y vectoresl de la forma Ax=B

\begin{pmatrix} 4 & 2 & 5 \\ 2 & 5 & 8 \\5 & 4 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 60.70 \\ 92.90 \\ 56.30 \end{pmatrix}

Para el algoritmo la matriz A y el vector B se escriben como arreglos.

A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

B = np.array([[60.70],
              [92.90],
              [56.30]])

Observe que:

  • Las matrices y vectores se ingresan con arreglos de la libreria numpy
  • el vector B se escribe en forma de columna
  • No se usan listas, de ser el caso se convierten hacia arreglos con np.array()

Si el vector B está como fila, se aumenta una dimensión [B] y se aplica la transpuesta

Bfila = np.array([60.70,92.90,56.30])
Bcolumna = np.transpose([Bfila])
print(Bcolumna)
>>> Bcolumna
array([[60.7],
       [92.9],
       [56.3]])

En el desarrollo de las soluciones, para mantener sincronía entre las operaciones entre filas de la matriz A y el vector B se usa la matriz aumentada.


Matriz Aumentada

Se realiza al concatenar la matriz A con el vector B en forma de columnas (axis=1).

AB = np.concatenate((A,B),axis=1)

el resultado AB se muestra como:

>>> AB
array([[ 4. ,  2. ,  5. , 60.7],
       [ 2. ,  5. ,  8. , 92.9],
       [ 5. ,  4. ,  3. , 56.3]])

Pivoteo parcial por filas

Para el pivoteo por fila de la matriz aumentada AB, tiene como primer paso revisar la primera columna desde la diagonal en adelante.

columna = [|4|,
           |2|,
           |5|]
dondemax = 2

El procedimiento de pivoteo se realiza si la posición dónde se encuentra el valor de  mayor magnitud no corresponde a la diagonal de la matriz (posición 0 de la columna).

En el ejercicio se encuentra que la magnitud de mayor valor está en la última fila, por lo que en AB se realiza el intercambio entre la fila 3 y la fila 1

AB = [[ 5. , 4. , 3. , 56.3],
      [ 2. , 5. , 8. , 92.9],
      [ 4. , 2. , 5. , 60.7]]

Se repite al paso anterior, pero para la segunda columna formada desde la diagonal.

columna = [|5|,
           |2|]
dondemax = 0

como la posición dondemax es la primera, índice 0, se determina que ya está en la diagonal de AB y no es necesario realizar el intercambio de filas.

Se repite el proceso para la tercera columna desde la diagonal, que resulta tener solo una casilla (columna =[5]) y no ser requiere continuar.

El resultado del pivoteo por fila se muestra como:

matriz pivoteada por fila:
AB = [[ 5. ,  4. ,  3. , 56.3],
      [ 2. ,  5. ,  8. , 92.9],
      [ 4. ,  2. ,  5. , 60.7]]

Algoritmo en Python

Para realizar el algoritmo, es de recordar que para realizar operaciones en una matriz sin alterar la original, se usa una copia de la matriz (np.copy). Se puede comparar y observar los cambios entre la matriz original y la copia a la que se aplicaron cambios

Si no es necesaria la comparación entre el antes y despues, no se realiza la copia y se ahorra el espacio de memoria, detalle importante para matrices de «gran tamaño» y una computadora con «limitada» memoria.

# Pivoteo parcial por filas
# Solución a Sistemas de Ecuaciones

import numpy as np

# INGRESO
A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

B = np.array([[60.70],
              [92.90],
              [56.30]])

# PROCEDIMIENTO
# Matriz aumentada
AB  = np.concatenate((A,B),axis=1)
AB0 = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = abs(AB[i:,i])
    dondemax = np.argmax(columna)
    
    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB)

Función pivoteafila(M)

Los bloques de cada procedimiento que se repiten en otros métodos se convierten a funciones def-return, empaquetando las soluciones algoritmicas a problemas resueltos.

Se usa la matriz M para generalizar y diferenciar de A que es usada en los ejercicios en realizados en adelante.

def pivoteafila(M):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    # Pivoteo por filas AB
    tamano = np.shape(M)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(M[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(M[i,:])
            M[i,:] = M[dondemax+i,:]
            M[dondemax+i,:] = temporal
    return(M)

2.6 Sistemas de Ecuaciones no lineales – Newton-Raphson

Referencia: Chapra 6.5 p162, Chapra Ejercicio 6.11 p166

Con el método de Newton-Raphson para múltiples ecuaciones, determine las raíces para:

x^2+xy =10 y + 3xy^2 = 57

Observe que un par correcto de raíces es x=2 y y=3.
Use como valore iniciales x=1.5, y=3.5

Planteamiento

Las ecuaciones se expresan de la forma f(x,y) = 0

x^2+xy -10 = 0 y + 3xy^2 -57 = 0

Se puede usar extensines de los métodos abiertos para resolver ecuacioens simples, por ejemplo Newton-Raphson.

u_{i+1} = u_i + (x_{i+1}-x_i)\frac{\partial u_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial u_i}{\partial y} v_{i+1} = v_i + (x_{i+1}-x_i)\frac{\partial v_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial v_i}{\partial y}

ecuaciones que se pueden reordenar y encontrar la solución a partir de la matriz Jacobiano.

Instrucciones en Python

Usando un algoritmo para resolver el Jacobiano y estimar los puntos luego de cada iteración se obtienen:

iteración:  1
Jacobiano con puntos iniciales: 
Matrix([[6.50000000000000, 1.50000000000000], [36.7500000000000, 32.5000000000000]])
determinante:  156.12499999999994
puntos xi,yi: 2.03602882305845 2.84387510008006
error: 0.656124899919936
iteración:  2
Jacobiano con puntos iniciales: 
Matrix([[6.91593274619696, 2.03602882305845], [24.2628767545662, 35.7412700376474]])
determinante:  197.78430344142245
puntos xi,yi: 1.99870060905582 3.00228856292451
error: 0.158413462844444
iteración:  3
Jacobiano con puntos iniciales: 
Matrix([[6.99968978103616, 1.99870060905582], [27.0412098452019, 37.0040558756713]])
determinante:  204.96962918261596
puntos xi,yi: 1.99999998387626 2.99999941338891
error: 0.00228914953559523
iteración:  4
Jacobiano con puntos iniciales: 
Matrix([[6.99999938114143, 1.99999998387626], [26.9999894410015, 36.9999926704397]])
determinante:  204.9999473486533
puntos xi,yi: 1.99999999999998 3.00000000000008
error: 5.86611161867978e-7
Resultado: 
1.99999999999998 3.00000000000008
>>> 

Algoritmo presentado para dos ecuaciones y dos incógnitas, en la unidad 3 se puede ampliar la propuesta. Revisar el método de Gauss-Seidel y Jacobi.

# Ejercicio Chapra Ej:6.11
# Sistemas de ecuaciones no lineales
# con método de Newton Raphson para xy

import numpy as np
import sympy as sym

def matrizJacobiano(variables, funciones):
    n = len(funciones)
    m = len(variables)
    # matriz Jacobiano inicia con ceros
    Jcb = sym.zeros(n,m)
    for i in range(0,n,1):
        unafi = sym.sympify(funciones[i])
        for j in range(0,m,1):
            unavariable = variables[j]
            Jcb[i,j] = sym.diff(unafi, unavariable)
    return Jcb

# PROGRAMA ----------
# INGRESO
x = sym.Symbol('x')
y = sym.Symbol('y')

f1 = x**2 + x*y - 10
f2 = y + 3*x*(y**2)-57

x0 = 1.5
y0 = 3.5

tolera = 0.0001

# PROCEDIMIENTO
funciones = [f1,f2]
variables = [x,y]
n = len(funciones)
m = len(variables)

Jxy = matrizJacobiano(variables, funciones)

# valores iniciales
xi = x0
yi = y0

# tramo inicial, mayor que tolerancia
itera = 0
tramo = tolera*2

while (tramo>tolera):
    J = Jxy.subs([(x,xi),(y,yi)])

    # determinante de J
    Jn = np.array(J,dtype=float)
    determinante =  np.linalg.det(Jn)

    # iteraciones
    f1i = f1.subs([(x,xi),(y,yi)])
    f2i = f2.subs([(x,xi),(y,yi)])

    numerador1 = f1i*Jn[n-1,m-1]-f2i*Jn[0,m-1]
    xi1 = xi - numerador1/determinante
    numerador2 = f2i*Jn[0,0]-f1i*Jn[n-1,0]
    yi1 = yi -numerador2/determinante
    
    tramo = np.max(np.abs([xi1-xi,yi1-yi]))
    xi = xi1
    yi = yi1

    itera = itera +1
    print('iteración: ',itera)
    print('Jacobiano con puntos iniciales: ')
    print(J)
    print('determinante: ', determinante)
    print('puntos xi,yi:',xi,yi)
    print('error:',tramo)
    
# SALIDA
print('Resultado: ')
print(xi,yi)

2.5.1 Método de la Secante – Ejemplo con Python

Referencia: Burden ejemplo 1 p51

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

[xa ,	 xb , 	 xc , 	 tramo]
[ 1.5     1.504   1.3736  0.1264]
[ 1.3736  1.5     1.3658  0.0078]
[  1.3658e+00   1.3736e+00   1.3652e+00   5.2085e-04]
raiz en:  1.36523214292

El algoritmo a implementar es:

# Método de la secante
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

import numpy as np

def secante_tabla(fx,xa,tolera):
    dx = 4*tolera
    xb = xa + dx
    tramo = dx
    tabla = []
    while (tramo>=tolera):
        fa = fx(xa)
        fb = fx(xb)
        xc = xa - fa*(xb-xa)/(fb-fa)
        tramo = abs(xc-xa)
        
        tabla.append([xa,xb,xc,tramo])
        xb = xa
        xa = xc

    tabla = np.array(tabla)
    return(tabla)

# PROGRAMA ---------------------
# INGRESO
fx = lambda x: x**3 + 4*x**2 - 10

a  = 1
b  = 2
xa = 1.5
tolera = 0.001
tramos = 100

# PROCEDIMIENTO
tabla = secante_tabla(fx,xa,tolera)
n = len(tabla)
raiz = tabla[n-1,2]

# SALIDA
np.set_printoptions(precision=4)
print('[xa ,\t xb , \t xc , \t tramo]')
for i in range(0,n,1):
    print(tabla[i])
print('raiz en: ', raiz)

En el caso de añadir la gráfica para la primera iteración:

# GRAFICA
import matplotlib.pyplot as plt

# Calcula los puntos a graficar
xi = np.linspace(a,b,tramos+1)
fi = fx(xi)
dx = (b-xa)/2
pendiente = (fx(xa+dx)-fx(xa))/(xa+dx-xa)
b0 = fx(xa) - pendiente*xa
tangentei = pendiente*xi+b0

fxa = fx(xa)
xb = xa + dx
fxb = fx(xb)

plt.plot(xi,fi, label='f(x)')

plt.plot(xi,tangentei, label='secante')
plt.plot(xa,fx(xa),'go', label='xa')
plt.plot(xa+dx,fx(xa+dx),'ro', label='xb')
plt.plot((-b0/pendiente),0,'yo', label='xc')

plt.plot([xa,xa],[0,fxa],'m')
plt.plot([xb,xb],[0,fxb],'m')

plt.axhline(0, color='k')
plt.title('Método de la Secante')
plt.legend()
plt.grid()
plt.show()

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


Tarea

convertir el algoritmo a una función de python con respuesta simple

2.4.1 Método del Punto fijo – Ejemplo con Python

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

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

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

Desarrollo Analítico

Al igual que los métodos anteriores, es conveniente determinar el intervalo [a,b] donde es posible evaluar f(x). Se revisa si hay cambio de signo en el intervalo para buscar una raíz.

Para el ejemplo, el rango de observación será [0,1], pues f(0)=1 es positivo y f(1)=-0.63 es negativo.

Para el punto fijo, se reordena la ecuación para para tener una ecuación con la variable independiente separada.

Se obtiene por un lado la recta identidad y=x, por otro se tiene la función g(x).

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

x = e^{-x} g(x) = e^{-x}

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

  • iniciando desde el extremo izquierdo x=a,
  • se determina el valor de b = g(x) ,
  • se determina la diferencia de la aproximación o error = |b-a|, tolera=0.001
  • se proyecta en la recta identidad, como el nuevo punto de evaluación
    x=b.
  • se repite el proceso para el nuevo valor de x, hasta que el error sea menor al tolerado.
  • En caso que el proceso no converge, se utiliza un contador de iteraciones máximo, para evitar tener un lazo infinito.

iteración 1

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

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.


Algoritmo en Python

El algoritmo en python se presenta como una función para usarla fácilmente como un bloque en otros ejercicios. Se requiere la función gx, el punto inicial y la tolerancia, la variable de número de iteraciones máxima, iteramax, permite controlar la convergéncia de la función con el método.

Se presenta el algoritmo mejorado, convirtiendo el procedimiento del método a una función Python.

# Algoritmo de punto fijo
# [a,b] intervalo de búsqueda
# error = tolera

import numpy as np

def puntofijo(gx,a,tolera, iteramax = 15):
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    while(tramo>=tolera and i<=iteramax ):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        i = i + 1
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax ):
        respuesta = np.nan
    return(respuesta)

# PROGRAMA ---------

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

a = 0       # intervalo
b = 1
tolera = 0.001
iteramax = 15  # itera máximo
muestras = 51  # gráfico
tramos = 50

# PROCEDIMIENTO
respuesta = puntofijo(gx,a,tolera)

# SALIDA
print(respuesta)

la respuesta obtenida del problema es

0.566908911921

Para obtener la gráfica básica se determinan los puntos para cada función fx y gx. Se añaden las siguiente instrucciones al algoritmo anterior:

# GRAFICA
# calcula los puntos para fx y gx
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
yi = xi

import matplotlib.pyplot as plt

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

if (respuesta!= np.nan):
    plt.axvline(respuesta)
plt.axhline(0, color='k')
plt.title('Punto Fijo')
plt.legend()
plt.show()

Scipy.optimize.fixed_point

El método del punto fijo se encuentra implementado en Scipy, que también puede ser usado de la forma:

>>> import scipy.optimize as opt
>>> opt.fixed_point(gx,a,xtol=0.001,maxiter=15)
array(0.5671432948307147)

el 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 raiz.
  • 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.
  • Realizar el mismo ejercicio para
f(x) = x^{2}-x-2 = 0

en el rango [-1,2]
Realice sus observaciones y recomendaciones.

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

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