Edison Del Rosario – Página 51 – Métodos numéricos

2017-09-032024-09-22

3.3.1 Método de Gauss – determinante de matriz con Python

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

1. Determinante de matriz

Referencia: Rodríguez 4.3.9 p129, Burden 6.4 p296, Chapra 9.1.2 p250.

El determinante de una matriz cuadrada triangular superior también puede calcularse como el producto de los coeficientes de la diagonal principal, considerando el número de cambios de fila del pivoteo k.

det(A) = (-1)^k \prod_{i=1}^n a_{i,i}

Si observamos que en las secciones anteriores se tiene desarrollado los algoritmos para obtener la matriz triangular superior en el método de Gauss, se usan como punto de partida para obtener los resultados del cálculo del determinante.

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

2. Ejercicio

Referencia: Chapra Ejemplo 9.12 p277

Calcular el determinante de la matriz A.

A= \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix}

El algoritmo para ejercicio se convierte en una extensión de los algoritmos anteriores.

A = np.array([[3. , -0.1, -0.2],
              [0.1,  7. , -0.3],
              [0.3, -0.2, 10.  ]])

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

3. Algoritmo en Python

El algoritmo parte de lo realizado en método de Gauss, indicando que la matriz a procesar es solamente A. Se mantienen los procedimientos de «pivoteo parcial por filas» y » eliminación hacia adelante»

Para contar el número de cambios de filas, se añade un contador de pivoteado dentro del condicional para intercambio de filas.

Para el resultado del operador multiplicación, se usan todas las casillas de la diagonal al acumular las multiplicaciones.

# Determinante de una matriz A
# convirtiendo a diagonal superior 

import numpy as np

# INGRESO
A = np.array([[3. , -0.1, -0.2],
              [0.1,  7. , -0.3],
              [0.3, -0.2, 10.  ]], dtype=float)

# PROCEDIMIENTO

# Matriz aumentada
AB = np.copy(A) # Para usar algoritmo previo

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

# Para cada fila en AB
pivoteado = 0 # contador para cambio fila

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

        pivoteado = pivoteado + 1 # cuenta cambio fila
        
AB1 = np.copy(AB)

# eliminación hacia adelante
for i in range(0,n-1,1):
    pivote   = AB[i,i]
    adelante = i + 1
    for k in range(adelante,n,1):
        factor  = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor

# calcula determinante
multiplica = 1
for i in range(0,n,1):
    multiplica = multiplica*AB[i,i]
determinante = ((-1)**pivoteado)*multiplica

# SALIDA
print('Pivoteo parcial por filas')
print(AB1)
print('Cambios de fila, pivoteado: ',pivoteado)
print('eliminación hacia adelante')
print(AB)
print('determinante: ')
print(determinante)

Se aplica la operación de la fórmula planteada para el método, y se presenta el resultado:

Pivoteo parcial por filas
[[ 3.  -0.1 -0.2]
 [ 0.1  7.  -0.3]
 [ 0.3 -0.2 10. ]]
Cambios de fila, pivoteado:  0
eliminación hacia adelante
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]
determinante: 
210.35299999999995
>>>

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

4. Algoritmo como función en Numpy

El resultado se puede verificar usando la función de Numpy:

>>> np.linalg.det(A)
210.3529999999999

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]

2017-09-032024-10-01

3.3 Método de Gauss con Python

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]

1. Método de Gauss

Referencia: Chapra 9.2 p254, Burden 6.1 p273, Rodríguez 4.3 p119

El método de Gauss opera sobre la matriz aumentada y pivoteada por filas, añadiendo los procesos:.

eliminación hacia adelante:

A_{k} = A_{k} -A_{i}\frac{a_{k,i}}{a_{i,i}}

sustitución hacia atrás:

x_i = \frac{b_i^{(i-1)}-\sum_{j=i+1}^{n}a_{ij}^{(i-1)} x_{j}}{a_{ii}^{(i-1)}}

para i = n-1, n-2, …

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

2. Ejercicio

Referencia: Rodríguez 4.0 p105, 1Eva_IT2010_T3_MN Precio artículos

Se continúa a partir del resultado del tema de pivoteo parcial por filas para matrices:

\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}

Matriz aumentada y pivoteada por filas:

[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

3. Eliminación hacia adelante o eliminación Gaussiana

Consiste en simplificar la matriz a una triangular superior, con ceros debajo de la diagonal, usando operaciones entre filas, para obtener:

Elimina hacia adelante
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

Los índices de fila y columna en la matriz A[i,j] se usan de forma semejante a la nomenclatura de los textos de Álgebra Lineal. Progresivamente para cada fila, se toma como referencia o pivote el elemento de la diagonal (i=j). Luego, se realizan operaciones con las filas inferiores para convertir los elementos por debajo de la diagonal en cero. Las operaciones incluyen el vector B debido a que se trabaja sobre la matriz aumentada AB.

AB Matriz aumentada y pivoteada por filas:
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]

iteración fila 1, operación fila 1 y 2

Para la fila 1, con posición i=0, se usa el elemento a_i,i como pivote.

pivote = AB[i,i] = AB[0,0] = 5

Para las filas de que están después de la diagonal se referencian como k.Se obtiene el factor escalar de la operación entre filas de la formula

k = i+1 = 0+1 = 1

A_{k} = A_{k} -A_{i}\frac{a_{k,i}}{pivote}

factor = AB[1,0]/pivote = 2/5

y se realiza la operación entre fila k y la fila i para actualizar la fila k,

       [ 2. 5.  8.  92.9]
-(2/5)*[ 5. 4.  3.  56.3]
__________________________
     = [ 0. 3.4 6.8 70.38]

con lo que la matriz aumentada AB se actualiza a:

AB =
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 4.    2.    5.   60.7 ]]

iteración fila 1, operación fila 1 y 3

se continúa con la siguiente fila, quedando la matriz aumentada con la columna debajo de la primera diagonal en cero:

k = i+1 = 2
factor = 4/5

        [ 4.  2.  5.   60.7] 
- (4/5)*[ 5.  4.  3.   56.3]
_____________________________
      = [ 0. -1.2 2.6  15.66]

AB =
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.   -1.2   2.6  15.66]]

Como ya se terminaron las operaciones con la primera posición de la diagonal, el siguiente paso es usar la segunda posición, i =2.

iteración fila 2

Para la fila 2, con posición i=1, se toma el elemento de la diagonal a_i,i como pivote, la variable adelante indica la referencia de la tercera fila

pivote = A[i,i] = AB[1,1] = 3.4

Para las filas ubicadas adelante de la diagonal se referencian como k

adelante = k = i+1 = 1+1 = 2

Para aplicar la fórmula por filas, se obtiene el factor .

factor = AB[2,1]/pivote  = -1.2/3.4 = - 0,3529

            [ 0. -1.2 2.6 15.66]
-(-1.2/3.4)*[ 0.  3.4 6.8 70.38]
________________________________
         =  [ 0.  0.  5.  40.5 ]

AB =
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

Con lo que se completa el objetivo de tener ceros debajo de la diagonal.
Observe que no es necesario realizar operaciones para la última fila, por lo que k debe llegar solamente hasta la fila penúltima.

El resultado de la eliminación hacia adelante a ser usado en el próximo paso es:

Elimina hacia adelante
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

4. Sustitución hacia atrás

La fórmula se interpreta para facilitar el algoritmo

x_i = \frac{b_i^{(i-1)}-\sum_{j=i+1}^{n}a_{ij}^{(i-1)} x_{j}}{a_{ii}^{(i-1)}}

Para una fila i, el vector b[i] representa el valor de la constante en la fila i de la matriz aumentada, a[i] se refiere los valores de los coeficientes de la ecuación, de los que se usan los que se encuentran a la derecha de la diagonal.

Las operaciones se realizan de abajo hacia arriba desde la última fila. Para el ejercicio presentado se tiene que:

ultfila = n-1 = 3-1 = 2
ultcolumna = m-1 = 4-1 = 3

la matriz a procesar es:

Elimina hacia adelante
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

iteración 1, fila 3, i=2

Empieza desde la última fila de la matriz,

[ 0. 0. 5. 40.5 ]

0 x_1 + 0 x_2 + 5 x_3 = 40.5

El valor de la constante es b = 40.5 y no existen elementos hacia la derecha de la diagonal. No se usa la ultima columna que es de las constantes:

5 x_3 = 40.5

x_3 = 40.5/5 = 8.1

la respuesta se interpreta en el vector X como:

X= \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 8.1 \end{pmatrix}

iteración 2, fila 2, i = 1

De la penúltima fila se obtiene la ecuación para encontrar x₂

[ 0. 3.4 6.8 70.38]

0x_1 + 3.4 x_2 +6.8 x_3 = 70.38

se observa que b = 70.38 y a la derecha de a diagonal se tiene un solo valor de [6.8].

3.4 x_2 = 70.38 -6.8 x_3

usa el valor de x₃ encontrado en la iteración anterior

3.4 x_2 = 70.38 -6.8 (8.1)

Muestra la ecuación para la segunda fila.

x_2 = (70.38 -6.8 (8.1))/3.4 = 4.5

la respuesta se interpreta en el vector X como:

X= \begin{pmatrix} 0 \\ 4.5 \\ 8.1 \end{pmatrix}

iteración 3 fila 1, i=0

se sigue el mismo proceso para la siguiente incógnita X₁ que se interpreta como

[ 5. 4. 3. 56.3 ]

5x_1 + 4 x_2 + 3x_3 = 56.3

5x_1 = 56.3 - ( 4 x_2 + 3x_3)

x_1 = \frac{56.3 - ( 4 x_2 + 3x_3)}{5}

Se encuentra que la solución al sistema de ecuaciones es:

X= \begin{pmatrix} 2.8\\ 4.5 \\ 8.1 \end{pmatrix}

por sustitución hacia atrás
el vector solución X es:
[[2.8]
 [4.5]
 [8.1]]

Verificar respuesta

Para verificar que el resultado es correcto, se usa el producto punto entre la matriz a y el vector resultado X. La operación A.X = B debe dar el vector B.

verificar que A.X = B
[[60.7]
 [92.9]
 [56.3]]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

5. Algoritmo con Python

El algoritmo para el Método de Gauss, reutiliza las instrucciones para matriz aumentada y pivoteo parcial por filas.

Recordar: Asegurar que los arreglos sean de tipo Real (float), para que no se considere el vector como entero y realice operaciones entre enteros, generando errores por truncamiento.

La parte nueva a desarrollar corresponde al procedimiento de «eliminación hacia adelante» y el procedimiento de «sustitución hacia atrás».

# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B

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
casicero = 1e-15 # Considerar como 0

# Evitar truncamiento en operaciones
A = np.array(A,dtype=float) 

# 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
AB1 = np.copy(AB)

# eliminación hacia adelante
for i in range(0,n-1,1):
    pivote   = AB[i,i]
    adelante = i + 1
    for k in range(adelante,n,1):
        factor  = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor

# sustitución hacia atrás
ultfila = n-1
ultcolumna = m-1
X = np.zeros(n,dtype=float)

for i in range(ultfila,0-1,-1):
    suma = 0
    for j in range(i+1,ultcolumna,1):
        suma = suma + AB[i,j]*X[j]
    b = AB[i,ultcolumna]
    X[i] = (b-suma)/AB[i,i]

X = np.transpose([X])

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB1)
print('eliminación hacia adelante')
print(AB)
print('solución: ')
print(X)

Tarea

Revisar cuando la matriz pivoteada por filas tienen un elemento cero o muy cercano a cero pues la matriz sería singular. El valor considerado como casi cero podría ser 1×10^-15

A estas alturas, por la cantidad de líneas de instrucción es recomendable reutilizar bloques de algoritmos usando funciones def-return. Por ejemplo: pivoteo por filas, eliminación hacia adelante, sustitución hacia atrás.

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]

6. Algoritmo como función de Python

El resultado par el ejercicio anterior es:

Matriz aumentada
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Elimina hacia adelante:
 fila 0 pivote:  5.0
   factor:  0.4  para fila:  1
   factor:  0.8  para fila:  2
 fila 1 pivote:  3.4
   factor:  -0.3529411764705883  para fila:  2
 fila 2 pivote:  5.0
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
solución: 
[2.8 4.5 8.1]
>>>

Instrucciones en Python

# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print(AB)
    return(AB)

def gauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss elimina hacia adelante
    tarea: verificar términos cero
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    if vertabla==True:
        print('Elimina hacia adelante:')
    for i in range(0,n,1):
        pivote = AB[i,i]
        adelante = i+1
        if vertabla==True:
            print(' fila',i,'pivote: ', pivote)
        for k in range(adelante,n,1):
            if (np.abs(pivote)>=casicero):
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    if vertabla==True:
        print(AB)
    return(AB)

def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss sustituye hacia atras
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    # Sustitución hacia atras
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1):
        suma = 0
        for j in range(i+1,ultcolumna,1):
            suma = suma + AB[i,j]*X[j]
        X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
    return(X)

# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

B = [60.70,92.90,56.30]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_sustituyeAtras(AB,vertabla=True)

# SALIDA
print('solución: ')
print(X)

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]

2017-09-032024-09-27

3.2 Pivoteo parcial por filas con Python

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]

Referencia: Chapra 9.4.2 p268. Burden 6.2 p279. Rodríguez 4.0 p105

Los métodos de solución a sistemas de ecuaciones en los primeros pasos usan la matriz aumentada y pivoteada por filas. Como es un procedimiento usado en todos los métodos de la unidad, para simplificar se presenta como uno de los primeros algoritmos.

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

1. Ejercicio

Referencia: Rodríguez 4.0 p105, 1Eva_IT2010_T3_MN Precio artículos

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	X₁	X₂	X₃	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 x₁, x₂, x₃ para escribir el sistema de ecuaciones que muestran las relaciones 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}

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

2. Matriz, Vector y la Matriz Aumentada

El sistema de ecuaciones se escribe en la forma algebraica como matrices y vectores 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 como arreglos de la librería 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 la solución, se usa la matriz aumentada, para mantener sincronía entre las operaciones entre filas de la matriz A y el vector B.

La matriz aumentada AB se forma al concatenar por columnas la matriz A con el vector B, usando el parámetro 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]])

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

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

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

4. 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 después, 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)

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

5. 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 algorítmicas a problemas resueltos.

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

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

import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas, entrega matriz aumentada AB
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print(AB)
    return(AB)

# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

B = [60.70, 92.90, 56.30]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

# SALIDA
print('Resultado de Pivoteo parcial por filas')
print(AB)

El resultado del ejercicio es:

Matriz aumentada
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
Pivoteo parcial por filas
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
>>>

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]

2017-09-032024-09-28

3.1 Sistema de ecuaciones 3×3, planos 3D con Python

2017-09-032024-10-08

Unidad 3 Solución a sistemas de ecuaciones

Pivoteo parcial por filas con Python

Métodos Directos

Gauss

Gauss-Jordan

Métodos Iterativos

Jacobi

Gauss-Seidel

2017-09-022024-10-03

2.7 Métodos de raíces con gráficos animados en Python

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para raíces de ecuaciones, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.

La gráfica (graf_ani) se crea en una ventana (fig_ani), inicializando con la linea a partir f(x) y configurando los parámetros base para el gráfico.

Se usan procedimientos para crear unatrama() para marca de iteración y en cada cambio se limpia la trama manteniendo la base con limpiatrama().

En caso de requerir un archivo .gif animado se proporciona un nombre de archivo. Para crear el archivo se requiere de un complemento ‘imagemagick‘ a ser instalado.

otros ejemplos de animación en el curso de Fundamentos de Programación:

Movimiento circular – Una partícula, animación con matplotlib-Python

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]
..

Método de la Bisección con gráfico animado en Python

# Algoritmo de Bisección, Tabla
# Los valores de [a,b] son aceptables
# y seleccionados desde la gráfica de la función
# error = tolera

import numpy as np

def biseccion_tabla(fx,a,b,tolera,iteramax = 20,
                    vertabla=False, precision=6):
    '''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)
    tabla=[]
    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)
            unafila = [a,c,b,fa,fc,fb] 
            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)
            unafila.append(tramo)
            tabla.append(unafila)
            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
    tabla = np.array(tabla,dtype=float)
    return(respuesta,tabla)

# PROGRAMA #######################

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

a = 1
b = 2
tolera = 0.001

# PROCEDIMIENTO
[raiz,tabla] = biseccion_tabla(fx,a,b,tolera,vertabla=True)

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

# GRAFICA
import matplotlib.pyplot as plt
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
xc = tabla[:,1]
yc = tabla[:,4]

plt.plot(xi,fi, label='f(x)')
plt.plot([a,b],[fx(a),fx(b)],'o',
         color='red',label='[[a,b],[f(a),f(b)]]')
plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Bisección')
plt.grid()
plt.legend()
#plt.show()

# GRAFICA CON ANIMACION ------------
import matplotlib.animation as animation

xa = tabla[:,0]
ya = tabla[:,3]
xb = tabla[:,2]
yb = tabla[:,5]
xc = tabla[:,1]
yc = tabla[:,4]

# Inicializa parametros de trama/foto
narchivo = 'Biseccion' # nombre archivo
retardo = 700 # milisegundos entre tramas
tramas = len(xa)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
ymax = np.max(fi)
ymin = np.min(fi)
deltax = np.abs(b-a)
deltay = np.abs(ymax-ymin)
graf_ani.set_xlim([a-0.05*deltax,b+0.05*deltax])
graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.05*deltay])
# Lineas y puntos base
lineafx,= graf_ani.plot(xi,fi, label='f(x)')

puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                        color='red', label='a')
puntob, = graf_ani.plot(xb[0], yb[0],'o',
                        color='green', label='b')
puntoc, = graf_ani.plot(xc[0], yc[0],'o',
                        color='orange', label='c')
lineaa, = graf_ani.plot([xa[0],xa[0]],
                        [0,ya[0]],color='red',
                        linestyle='dashed')
lineab, = graf_ani.plot([xb[0],xb[0]],
                        [0,yb[0]],color='green',
                        linestyle='dashed')
lineac, = graf_ani.plot([xc[0],xc[0]],
                        [0,yc[0]],
                        color='orange',
                        linestyle='dashed')
# Configura gráfica
linea0 = graf_ani.axhline(0, color='k')
graf_ani.set_title('Bisección')
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Cada nueva trama
def unatrama(i,xa,ya,xb,yb,xc,yc):
    # actualiza cada punto
    puntoa.set_xdata(xa[i]) 
    puntoa.set_ydata(ya[i])
    puntob.set_xdata(xb[i])
    puntob.set_ydata(yb[i])
    puntoc.set_xdata(xc[i])
    puntoc.set_ydata(yc[i])
    # actualiza cada linea
    lineaa.set_ydata([ya[i], 0])
    lineaa.set_xdata([xa[i], xa[i]])
    lineab.set_ydata([yb[i], 0])
    lineab.set_xdata([xb[i], xb[i]])
    lineac.set_ydata([yc[i], 0])
    lineac.set_xdata([xc[i], xc[i]])
    return (puntoa, puntob, puntoc, lineaa, lineab, lineac,)

# Limpia trama anterior
def limpiatrama():
    puntoa.set_ydata(np.ma.array(xa, mask=True))
    puntob.set_ydata(np.ma.array(xb, mask=True))
    puntoc.set_ydata(np.ma.array(xc, mask=True))
    lineaa.set_ydata(np.ma.array([0,0], mask=True))
    lineab.set_ydata(np.ma.array([0,0], mask=True))
    lineac.set_ydata(np.ma.array([0,0], mask=True))
    return (puntoa, puntob, puntoc, lineaa, lineab, lineac,)

# contador de tramas
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(fig_ani,unatrama,
                              i ,
                              fargs=(xa, ya,
                                     xb, yb,
                                     xc, yc),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'_animado.mp4')
plt.show()

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

Método de la Posición Falsa con gráfico animado 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_tabla(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)
    tabla = []
    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)
            unafila = [a,c,b,fa, fc, fb]
            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
            unafila.append(tramo)
            tabla.append(unafila)
            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
    tabla = np.array(tabla,dtype=float)
    return(respuesta,tabla)

# PROGRAMA ----------------------

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

a = 1
b = 2
tolera = 0.001

# PROCEDIMIENTO
[raiz,tabla] = posicionfalsa_tabla(fx,a,b,tolera,
                                   vertabla=True)
# SALIDA
print('raíz en: ', raiz)

# GRAFICA
import matplotlib.pyplot as plt
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
xc = tabla[:,1]
yc = tabla[:,4]

plt.plot(xi,fi, label='f(x)')
plt.plot([a,b],[fx(a),fx(b)],'o',
         color='red',label='[[a,b],[f(a),f(b)]]')
plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Posición Falsa')
plt.grid()
plt.legend()
#plt.show()

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

xa = tabla[:,0]
ya = tabla[:,3]
xc = tabla[:,1]
yc = tabla[:,4]
xb = tabla[:,2]
yb = tabla[:,5]

# Inicializa parametros de trama/foto
narchivo = 'PosicionFalsa' # nombre archivo
retardo = 700 # milisegundos entre tramas
tramas = len(xa)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
graf_ani.set_xlim([a,b])
graf_ani.set_ylim([np.min(fi),np.max(fi)])
# Lineas y puntos base
lineafx, = graf_ani.plot(xi,fi,label ='f(x)')

puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                        color='red', label='a')
puntob, = graf_ani.plot(xb[0], yb[0],'o',
                        color='green', label='b')
puntoc, = graf_ani.plot(xc[0], yc[0],'o',
                        color='orange', label='c')

lineaab, = graf_ani.plot([xa[0],xb[0]],[ya[0],yb[0]],
                        color ='orange',
                        label='y=mx+b')
lineac0, = graf_ani.plot([xc[0],xc[0]],[0,yc[0]],
                        color='magenta',
                        linestyle='dashed')

# Configura gráfica
linea0 = graf_ani.axhline(0, color='k')
graf_ani.set_title('Posición Falsa')
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Cada nueva trama
def unatrama(i,xa,ya,xc,yc,xb,yb):
    # actualiza cada punto
    puntoa.set_xdata(xa[i])
    puntoa.set_ydata(ya[i])
    puntob.set_xdata(xb[i])
    puntob.set_ydata(yb[i])
    puntoc.set_xdata(xc[i])
    puntoc.set_ydata(yc[i])
    # actualiza cada linea
    lineaab.set_ydata([ya[i], yb[i]])
    lineaab.set_xdata([xa[i], xb[i]])
    lineac0.set_ydata([0, yc[i]])
    lineac0.set_xdata([xc[i], xc[i]])  
    return (puntoa, puntob, puntoc, lineaab,lineac0,)

# Cada nueva trama
def limpiatrama():
    puntoa.set_ydata(np.ma.array(xa, mask=True))
    puntob.set_ydata(np.ma.array(xb, mask=True))
    puntoc.set_ydata(np.ma.array(xc, mask=True))
    lineaab.set_ydata(np.ma.array([0,0], mask=True))
    lineac0.set_ydata(np.ma.array([0,0], mask=True))
    return (puntoa, puntob, puntoc, lineaab,lineac0,)

# contador de tramas
i = np.arange(0, tramas,1)
ani = animation.FuncAnimation(fig_ani,unatrama,
                              i ,
                              fargs=(xa, ya,
                                     xc, yc,
                                     xb,yb),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')
plt.show()

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]
..

Método de Newton-Raphson con gráfico animado en Python

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

def newton_raphson_tabla(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)
    tabla = []
    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]))

        tabla.append([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')
    tabla = np.array(tabla,dtype=float)
    return(xi,tabla)

# PROGRAMA ----------------------

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

a = 1
b = 4
x0 = 3
tolera = 0.001

# PROCEDIMIENTO
[raiz,tabla] = newton_raphson_tabla(fx, dfx, x0,
                             tolera, vertabla=True)
# SALIDA
print('raiz en :', raiz)

# GRAFICA
import matplotlib.pyplot as plt
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
xc = tabla[:,3]
yc = fx(xc)

plt.plot(xi,fi, label='f(x)')
plt.plot(x0,fx(x0),'o',
         color='red',label='[x0,f(x0)]')
plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
#plt.show()

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

xa = tabla[:,0]
ya = tabla[:,1]
xb = tabla[:,3]
dfi = tabla[:,2]
# Aproximacion con tangente
b0 = ya[0] - dfi[0]*x0
tangentei = dfi[0]*xi + b0
ci = -b0/dfi[0]

# Inicializa parametros de trama/foto
narchivo = 'NewtonRaphson' # nombre archivo
retardo = 700 # milisegundos entre tramas
tramas = len(xa)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
graf_ani.set_xlim([a,b])
graf_ani.set_ylim([np.min(fi),np.max(fi)])
# Lineas y puntos base
lineafx, = graf_ani.plot(xi,fi,label ='f(x)')

puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                        color='red', label='x[i]')
puntob, = graf_ani.plot(xb[0], 0,'o',
                        color='green', label='x[i+1]')

lineatanx, = graf_ani.plot(xi,tangentei,
                           color='orange',label='tangente')
lineaa, = graf_ani.plot([xa[0],xa[0]],
                        [ya[0],0], color='magenta',
                        linestyle='dashed')
lineab, = graf_ani.plot([xb[0],xb[0]],
                         [0,fx(xb[0])], color='magenta',
                         linestyle='dashed')
# Configura gráfica
linea0 = graf_ani.axhline(0, color='k')
graf_ani.set_title('Newton-Raphson')
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Cada nueva trama
def unatrama(i,xa,ya,xb,dfi):
    # actualiza cada punto
    puntoa.set_xdata(xa[i]) 
    puntoa.set_ydata(ya[i])
    puntob.set_xdata(xb[i])
    puntob.set_ydata(0)
    # actualiza cada linea
    lineaa.set_ydata([ya[i], 0])
    lineaa.set_xdata([xa[i], xa[i]])
    lineab.set_ydata([0, fx(xb[i])])
    lineab.set_xdata([xb[i], xb[i]])
    # Aproximacion con tangente
    b0 = ya[i] - dfi[i]*xa[i]
    tangentei = dfi[i]*xi+b0
    lineatanx.set_ydata(tangentei)
    
    return (puntoa, puntob, lineaa, lineab,lineatanx,)

# Limpia trama anterior
def limpiatrama():
    puntoa.set_ydata(np.ma.array(xa, mask=True))
    puntob.set_ydata(np.ma.array(xb, mask=True))
    lineaa.set_ydata(np.ma.array([0,0], mask=True))
    lineab.set_ydata(np.ma.array([0,0], mask=True))
    lineatanx.set_ydata(np.ma.array(xi, mask=True))
    return (puntoa, puntob, lineaa, lineab,lineatanx,)

# contador de tramas
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(fig_ani,unatrama,
                              i ,
                              fargs=(xa, ya,
                                     xb, dfi),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')
plt.show()

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

Método del Punto Fijo con gráfico animado en Python

Instrucciones en Python

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

import numpy as np

def puntofijo_tabla(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)
    tabla = [[a,b,tramo]]
    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]))
        tabla.append([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')
    tabla = np.array(tabla,dtype=float)
    return(respuesta,tabla)

# PROGRAMA ----------------------

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

a = 0
b = 1
tolera = 0.001

# PROCEDIMIENTO
[raiz,tabla] = puntofijo_tabla(gx,a,tolera, vertabla=True)

# SALIDA
print('raiz en :', raiz)

# GRAFICA
import matplotlib.pyplot as plt
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(raiz, color='magenta',
            linestyle='dotted')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Punto Fijo')
plt.grid()
plt.legend()
#plt.show()

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

xc = tabla[:,0]
yc = tabla[:,1]

# Inicializa parametros de trama/foto
narchivo = 'PuntoFijo' # nombre archivo
retardo = 700 # milisegundos entre tramas
tramas = len(xc)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
dx = np.abs(b-a)
dy = np.abs(gx(b)-gx(a))
graf_ani.set_xlim([a-0.05*dx,b+0.05*dx])
graf_ani.set_ylim([np.min(fi)-0.05*dy,np.max(fi)+0.05*dy])

# Lineas y puntos base
puntoa,  = graf_ani.plot(xc[0], yc[0], 'ro')
puntob,  = graf_ani.plot(xc[0], xc[0], 'go') # y=x
lineafx, = graf_ani.plot(xi,fi, color='blue',
                     linestyle='dashed', label='f(x)')
lineagx, = graf_ani.plot(xi,gi, color='orange', label='g(x)')
lineayx, = graf_ani.plot(xi,yi, color='green', label='y=x')
linearaiz = graf_ani.axvline(raiz, color='magenta',linestyle='dotted')

# Configura gráfica
linea0 = graf_ani.axhline(0, color='k')
graf_ani.set_title('Punto Fijo')
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Cada nueva trama
def unatrama(i,xc,yc):
    # actualiza cada punto
    puntoa.set_xdata(xc[i])
    puntoa.set_ydata(yc[i])
    puntob.set_xdata(xc[i]) # y=x
    puntob.set_ydata(xc[i])
    
    # actualiza cada flecha
    dx = xc[i+1]-xc[i]
    dy = yc[i]-xc[i]
    flecha01 = graf_ani.arrow(xc[i],xc[i], 0,dy,
              length_includes_head = True,
              head_width = 0.05*abs(dy),
              head_length = 0.1*abs(dy))
    flecha02 = graf_ani.arrow(xc[i],yc[i], dx,0,
              length_includes_head = True,
              head_width = 0.05*abs(dx),
              head_length = 0.1*abs(dx))
    return (puntoa, puntob, flecha01, flecha02)

# Limpia trama anterior
def limpiatrama():
    puntoa.set_ydata(np.ma.array(xc[0], mask=True))
    puntob.set_ydata(np.ma.array(xc[0], mask=True))
    return (puntoa, puntob)

# contador de tramas
i = np.arange(0, tramas-1,1)
ani = animation.FuncAnimation(fig_ani,unatrama,
                              i ,
                              fargs=(xc, yc),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')
plt.show()

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

Método de la Secante con gráfico animado en Python

Instrucciones en Python

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

def secante_raiz_tabla(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)
    tabla = []
    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)
        tabla.append([xi_1,xi,xi1,fi_1,fi])
        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
    tabla = np.array(tabla,dtype=float)
    return(respuesta,tabla)

# PROGRAMA ----------------------

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

a = 1
b = 4
c = (a+b)/2
tolera = 0.001
tramos = 100

# PROCEDIMIENTO
[raiz,tabla] = secante_raiz_tabla(fx,c,b,tolera,
                                  vertabla=True)
# SALIDA
print('raiz en :', raiz)

# GRAFICA
import matplotlib.pyplot as plt
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
xc = tabla[:,2]
yc = fx(xc)

plt.plot(xi,fi, label='f(x)')
plt.plot([a,b],[fx(a),fx(b)],'o',
         color='red',label='[x0,f(x0)]')
plt.scatter(xc,yc,color='orange', label='[c,f(c)]')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
#plt.show()


# GRÁFICO CON ANIMACION ######
import matplotlib.animation as animation

xa = tabla[:,0]
ya = tabla[:,3]
xb = tabla[:,1]
yb = tabla[:,4]
xc = tabla[:,2]

# Inicializa parametros de trama/foto
narchivo = 'SecanteMetodo' # nombre archivo
retardo = 700 # milisegundos entre tramas
tramas = len(xa)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
graf_ani.set_xlim([a,b])
graf_ani.set_ylim([np.min(fi),np.max(fi)])

# Lineas y puntos base
lineafx, = graf_ani.plot(xi,fi,label ='f(x)')

puntoa, = graf_ani.plot(xa[0], ya[0],'o',
                        color='red', label='x[i-1]')
puntob, = graf_ani.plot(xb[0], yb[0],'o',
                        color='green', label='x[i]')
puntoc, = graf_ani.plot(xc[0], yc[0],'o',
                        color='orange', label='x[i+1]')

lineatan1, = graf_ani.plot([xa[0],xb[0]],
                           [ya[0],yb[0]],
                           color='orange',label='secante ac')
lineatan2, = graf_ani.plot([xc[0],xb[0]],
                           [0,yb[0]],
                           color='orange',label='secante cb')
linea_a, = graf_ani.plot([xa[0],xa[0]],
                         [ya[0],0], color='magenta',
                         linestyle='dashed')
linea_b, = graf_ani.plot([xb[0],xb[0]],
                         [0,yb[0]], color='magenta',
                         linestyle='dashed')
# Configura gráfica
linea0 = graf_ani.axhline(0, color='k')
graf_ani.set_title('Método de la Secante')
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Cada nueva trama
def unatrama(i,xa,ya,xb,yb,xc):
    # actualiza cada punto
    puntoa.set_xdata(xa[i]) 
    puntoa.set_ydata(ya[i])
    puntob.set_xdata(xb[i])
    puntob.set_ydata(yb[i])
    puntoc.set_xdata(xc[i])
    puntoc.set_ydata(0)
    # actualiza cada linea
    linea_a.set_ydata([ya[i], 0])
    linea_a.set_xdata([xa[i], xa[i]])
    linea_b.set_ydata([0, yb[i]])
    linea_b.set_xdata([xb[i], xb[i]])
    lineatan1.set_ydata([ya[i], 0])
    lineatan1.set_xdata([xa[i], xc[i]])
    lineatan2.set_ydata([0, yb[i]])
    lineatan2.set_xdata([xc[i], xb[i]])

    return (puntoa, puntob, puntoc,
            linea_a, linea_b,
            lineatan1, lineatan2,)

# Limpia trama anterior
def limpiatrama():
    puntoa.set_ydata(np.ma.array(xa, mask=True))
    puntob.set_ydata(np.ma.array(xb, mask=True))
    puntoc.set_ydata(np.ma.array(xc, mask=True))
    linea_a.set_ydata(np.ma.array([0,0], mask=True))
    linea_b.set_ydata(np.ma.array([0,0], mask=True))
    lineatan1.set_ydata(np.ma.array([0,0], mask=True))
    lineatan2.set_ydata(np.ma.array([0,0], mask=True))
    return (puntoa, puntob, puntoc,
            linea_a, linea_b,
            lineatan1, lineatan2,)

# contador de tramas
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(fig_ani,unatrama,
                              i ,
                              fargs=(xa, ya,
                                     xb, yb,
                                     xc),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')
plt.show()

GIF animado: [ Bisección ] [ Posicion Falsa ] [ Newton-Raphson ] [ Punto Fijo ] [ Secante ]

2017-09-022024-09-28

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 valores 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 extensiones de los métodos abiertos para resolver ecuaciones 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: 
[ 6.5   1.5 ]
[           ]
[36.75  32.5]
determinante:  156.12499999999994
puntos xi,yi: 2.03602882305845 2.84387510008006
error: 0.656124899919936
iteración:  2
Jacobiano con puntos iniciales: 
[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: 
[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: 
[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: ')
    sym.pprint(J)
    print('determinante: ', determinante)
    print('puntos xi,yi:',xi,yi)
    print('error:',tramo)
    
# SALIDA
print('Resultado: ')
print(xi,yi)

2017-09-022024-09-30

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 ] [ 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, X₀ =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

X₀ = 4, X₁ = 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

X₁ =1.9742, X₂ =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 ]

2017-09-022024-10-08

2.5 Método de la Secante – Concepto

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

Método de la Secante

Referencia: Chapra 6.3 p154, Burden 2.3 p54,

El método de la secante busca evitar un posible inconveniente en el desarrollo el algoritmo para método de Newton-Raphson que es el implementar la evaluación de la expresión y valor de la primera derivada.

La derivada de la función f(x) también se puede aproximar mediante una diferencia finita dividida hacia atrás:

f'(x_i) = \frac{f(x_{i-1})-f(x_i)}{x_{i-1}-x_i}

la que se sustituye en la ecuación del método de Newton-Raphson para obtener:

x_{i+1}= x_i - f(x_i)\frac{(x_{i-1} - x_i)}{f(x_{i-1}) - f(x_i)}

Los métodos de la secante y de la falsa posición tienen ecuaciones idénticas, usan dos valores iniciales x[i-1] , x[i], para proyectar el nuevo valor de x[i+1].

Sin embargo, existe una diferencia importante entre ambos métodos en la forma en que uno de los valores iniciales se reemplaza por la nueva aproximación.

En el método de la falsa posición, la última aproximación de la raíz reemplaza cualquiera de los valores iniciales que dé un valor de la función con el mismo signo. En consecuencia, las dos aproximaciones siempre encierran a la raíz, es un método cerrado y siempre converge.

En el método de la secante se reemplaza los valores en secuencia estricta: con el nuevo valor x[i+1] se reemplaza a xi y xi reemplaza a x[i – 1]. Por lo que, algunas veces los dos valores están en el mismo lado de la raíz y en ciertos casos esto puede llevar a divergencias.

Observación: ¿Cuál es la diferencia con el método de Newton-Raphson?

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

2017-09-022024-10-04

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 ] [ 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}$

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 ]