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.
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.
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, en la sección de pivoteo se añade un contador cambiofilas en el condicional de cambio de filas.
Para el resultado del operador multiplicación, se usan todas las casillas de la diagonal al acumular las multiplicaciones.
Se aplica la operación de la fórmula planteada para el método, y se presenta el resultado.
# 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. ]])
# PROCEDIMIENTO# Matriz aumentada
AB = np.copy(A)
# Pivoteo parcial por filas
cambiofila = 0 # contador
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (dondemax !=0):
# intercambia filas
cambiofila = cambiofila +1
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
AB1 = np.copy(AB)
# eliminación hacia adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i + 1
for k inrange(adelante,n,1):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - AB[i,:]*factor
# calcula determinante
multiplica = 1
for i inrange(0,n,1):
multiplica = multiplica*AB[i,i]
determinante = ((-1)**cambiofila)*multiplica
# SALIDAprint('Pivoteo parcial por filas')
print(AB1)
print('eliminación hacia adelante')
print(AB)
print('determinante: ')
print(determinante)
El método de Gauss opera sobre la matriz aumentada y pivoteada por filas, añadiendo el proceso de «eliminación hacia adelante» mediante la operación entre filas. Se continúa entonces desde el resultado del tema de 3.2 pivoteo parcial por filas para matrices:
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.
Los índices de fila y columna, 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.
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:
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:
realice las operaciones con los valores encontrados para X2 y X3
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 atras
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]]
Método de Gauss con Algoritmo en Python
El algoritmo en su primera parte reutiliza lo desarrollado en Python para la 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=Bimport 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 ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (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 adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i + 1
for k inrange(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 inrange(ultfila,0-1,-1):
suma = 0
for j inrange(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])
# SALIDAprint('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
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:
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
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 Ecuacionesimport 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 ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (dondemax !=0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
# SALIDAprint('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.
defpivoteafila(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 ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(M[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (dondemax != 0):
# intercambia filas
temporal = np.copy(M[i,:])
M[i,:] = M[dondemax+i,:]
M[dondemax+i,:] = temporal
return(M)
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
# Método de la secante# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)import numpy as np
defsecante_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 inrange(0,n,1):
print(tabla[i])
print('raiz en: ', raiz)
En el caso de añadir la gráfica para la primera iteración:
# GRAFICAimport 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
Un posible inconveniente en el método de Newton Raphson es implementar la evaluación de la derivada. La derivada 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:
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.
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 = toleraimport numpy as np
defpuntofijo(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 respuestaif (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)
# SALIDAprint(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 método del punto fijo es un método abierto, también llamado de iteración de un punto o sustitución sucesiva, que reordena la ecuación
f(x)=0
de la forma en que x esté del lado izquierdo de la ecuación, para buscarla intersección entre la recta identidad y la curva g(x), como se muestran en los siguientes ejemplos.
y = x
x=g(x)
Observe que la raiz de f(x) se encuentra en el mismo valor de x donde ocurre la intersección entre la recta identidad en color verde y la función g(x) en color naranja. Se usa la linea vertical en color morado en x=raiz como referencia de lo indicado.
El método consiste en establecer un punto inicia x0 para la búsqueda, que se usa para calcular el valor g(x0).
En la siguiente iteración el nuevo valor para x es g(x0), que se refleja en la recta identidad y nuevamente se usa para calcular g(x).
El resultado iterativo se muestra en la figura animada, donde se observa que el resultado es convergente.
Ejemplo 1
f(x):e^{-x} - x = 0
se reordena para tener:
x = e^{-x} g(x) = e^{-x}
Ejemplo 2
f(x): x^2 - 2x -3 = 0
se reordena para tener:
x = \frac {x^2 - 3}{2} g(x) = \frac {x^2 - 3}{2}
Ejemplo 3
f(x): \sin (x) = 0
puede ser complicado despejar x, por lo que se simplifica el proceso sumando x en ambos lados.
x = \sin (x) + x g(x) = \sin (x) + x
El método proporciona una fórmula para predecir un valor nuevo de x en función del valor anterior: