Resolver Sistema Tridiagonal Disperso Con Scipy, Python

Hola , como paso previo para poder resolver una E.D.P del calor se necesita saber resolver sistemas tridiagonales lineales dispersos. Con Python es relativamente sencillo si utilizamos las librerías Scipy, con las cuales podemos implementar matrices dispersas con funciones específicas como dia_matrix.
Una vez implementada la matriz dispersa de nuestro sistema tan sólo hay que elegir el método para resolverlo: descomposición LU, Choleski, etc. En mi caso he utilizado la descomposición LU (también conocida como Algoritmo de Crout) y que viene implementada con Scipy. Éste método es eficiente para resolver la ecuación del calor, además; nuestra matriz de coeficientes dispersa de la ecuación del calor es tridiagonal y cumpliendo las condiciones necesarias para que funcione el algoritmo de Crout: semidefinida positiva, etc.
He realizado un pequeño algoritmo en Python 3 y Scipy en el que vemos cómo resolver primero un sistema lineal cualquiera 5 x 5 con Crout y luego cómo aplicarlo con un sistema 7 x 7 en el que nuestra matriz de coeficientes es una matriz dispersa tridiagonal.
Pongo el algoritmo aquí con todo explicado en el mismo código, lo he escrito en inglés, pero cualquier conocedor del tema no debe tener complicaciones para entenderlo.

  1. from scipy import linalg, sparse, array
  2.  
  3.  
  4. def crout(A, b):
  5.     """
  6.    Find out the solution of a linear system using the LU decomposition, also knnown
  7.    as The Crout Algorithm.
  8.    :param A: Coefficient Matrix
  9.    :param b: Independent column matrix terms
  10.    :return: A solution matrix of a linear system
  11.    """
  12.     lu = linalg.lu_factor(A)
  13.     solution = linalg.lu_solve(lu, b)
  14.     return solution
  15.  
  16.  
  17. def tridiagonal_system(data, diagonals, rows, columns, b):
  18.     """
  19.    Find out the solution of a sparse tri-diagonal linear system using The Crout Algorithm
  20.    :param data: the diagonal elements of a coefficient tri-diagonal matrix
  21.    :param diagonals: a vector to specify the positions of the diagonals in a tri-diagonal matrix
  22.    :param rows: the number of equations of the linear system
  23.    :param columns: the number of unknowns of the linear system
  24.    :param b: Independent column matrix terms
  25.    :return: A solution matrix of a sparse tri-diagonal linear system
  26.    """
  27.     A = sparse.dia_matrix((data, diagonals), shape=(rows, columns)).todense()
  28.     return crout(A, b)
  29.  
  30.  
  31. def main():
  32.     A = array([[4, 5, 9, 2, 8], [9, 7, 8, 12, 3], [2, 5, 7, 2, 4], [2, 6, 8, 4, 8], [1, 7, 4, 3, 7]])
  33.     b = array([2, 9, 7, 4, 5])
  34.  
  35.     print('The solution of the linear system using The Crout Algorithm is:')
  36.     print(crout(A, b))
  37.  
  38.     diag_down = array([1]).repeat(7)  # The -1 diagonal
  39.     diag_main = array([2]).repeat(7)  # The main diagonal
  40.     diag_up = array([3]).repeat(7)  # The 1 diagonal
  41.     diagonals = array([-1, 0, 1])  # The positions of our diagonals in the coefficient mattrix
  42.     data = array([diag_down, diag_main, diag_up])
  43.  
  44.     b2 = array([2, 9, 7, 4, 5, 6, 3])
  45.     print('The solution of the tridiagonal linear system using The Crout Algorithm is:')
  46.     print(tridiagonal_system(data, diagonals, 7, 7, b2))
  47.  
  48. main()

Se puede apreciar que el código es corto y lo suficientemente claro para poder aplicarlo con eficiencia. La salida del mismo es la siguiente:

usuario@nombrepc:~$
The solution of the linear system using The Crout Algorithm is:
[-0.46742021 1.76130319 0.64827128 -0.02393617 -1.34042553]
The solution of the tridiagonal linear system using The Crout Algorithm is:
[-5.16071429 4.10714286 1.98214286 -0.35714286 0.91071429 1.17857143
0.91071429]

En cuanto tenga tiempo escribiré sobre cómo resolver y dibujar la ecuación del calor, reutilizando éste código, aplicando el Algoritmo de Crank-Nicolson y las librerías Matplotlib de Python 3.
Saludos

Pingbacks/Trackbacks

  1. Bitacoras.com - 4 septiembre, 2014

    Información Bitacoras.com

    Valora en Bitacoras.com: Hola , como paso previo para poder resolver una E.D.P del calor se necesita saber resolver sistemas tridiagonales lineales dispersos. Con Python es relativamente sencillo si utilizamos las librerías Scipy, con las cuales pod…

Comments are closed.

Proudly powered by WordPress   Premium Style Theme by www.gopiplus.com