Resolviendo Sistemas Lineales De 3 Incógnitas Con Sympy Y Matplotlib

por | 1 febrero, 2014

Hola
Veamos cómo podemos dibujar con Matplotlib un sistema de ecuaciones lineales que son 3 planos, junto con su solución. El sistema de ejemplo va a ser el siguiente:
\left\{\begin{array}{lrr} x-2y+2z=0\\ 2x-y+z=0\\ x-y+z=0\end{array}\right\}
El objetivo es poder mostrar a nuestros alumnos de Bachiller de forma gráfica (y limpia) las soluciones de sistemas lineales, con Matplotlib y de forma rápida obtenemos dibujos de calidad que luego podemos mostrar a nuestros alumnos con una pizarra digital. De esa forma los alumnos tienen una mejor visión y comprensión de qué es la solución de estos sistemas; pueden comprender mejor lo de Incompatible, Compatible Determinado e Indeterminado.
El sistema al resolverlo nos damos cuenta que es Compatible Indeterminado, infinitas soluciones. Gráficamente son 3 planos que se cortan en una recta (0,\lambda,\lambda)\; \forall \; \lambda\in\mathbb{R}.
Si queremos podemos hacer que en nuestro mismo código de Python que utilizamos para dibujar nos solucione el sistema analíticamente utilizando las librerías matemáticas simbólicas Sympy de Python. Por ejemplo, para nuestro sistema sería:

  1. from sympy.solvers import *
  2. from sympy import *
  3.  
  4. x = Symbol('x')
  5. y = Symbol('y')
  6. z = Symbol('z')
  7.  
  8. fun1 = x -2*y + 2*z
  9. fun2 = 2*x - y + z
  10. fun3 = x - y + z
  11. solucion = solve([fun1, fun2, fun3], [x, y, z])
  12.  
  13. pprint(('Solución Del Sistema es: {}').format(solucion))

Si eres usuario de software matemático no debes tener problema alguno para entender el código. Por terminal nos da lo siguiente:

usuario@nombrepc:~$
Solución Del Sistema es: {y: z, x: 0}

Que es la misma solución que he puesto arriba pero escrita de otra forma. Tan sólo queda dibujar la figura en 3D con Matplotlib. Pongo todo el código con comentarios:

  1. #/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3.  
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7. from mpl_toolkits.mplot3d import Axes3D
  8. from sympy.solvers import *
  9. from sympy import *
  10. from matplotlib import rcParams
  11.  
  12.  
  13. # Activating LateX
  14. rcParams['text.latex.unicode'] = True
  15. rcParams['text.usetex'] = True
  16. rcParams['text.latex.preamble'] = '\\usepackage{amsthm}', '\\usepackage{amsmath}', '\\usepackage{amssymb}',
  17. '\\usepackage{amsfonts}', '\\usepackage[T1]{fontenc}', '\\usepackage[utf8]{inputenc}'
  18.  
  19. # Declaring the three planes as functions
  20. f1 = lambda x, y: (-x + 2 * y) / 2.
  21. f2 = lambda x, y: -2 * x + y
  22. f3 = lambda x, y: -x + y
  23.  
  24. # Declaring symbolic variables
  25. x = Symbol('x')
  26. y = Symbol('y')
  27. z = Symbol('z')
  28.  
  29. # Solving the linear system
  30. fun1 = x -2*y + 2*z
  31. fun2 = 2*x - y + z
  32. fun3 = x - y + z
  33. solucion = solve([fun1, fun2, fun3], [x, y, z])
  34.  
  35. # Printing the solution
  36. pprint(('Solución Del Sistema es: {}').format(solucion))
  37.  
  38.  
  39. # Stablishing our ranges for our variables
  40. x1 = y1 = np.arange(-6, 6, 0.25)
  41. ceros = np.zeros(len(x1))
  42.  
  43. # Stablishing our meshgrid
  44. x, y = np.meshgrid(x1, y1)
  45.  
  46. # Our 3D Canvas Figure Plot
  47. fig = plt.figure()
  48. ax = fig.add_subplot(111, projection='3d')
  49.  
  50. # Plotting the solution
  51. ax.plot(ceros, y1, zs=y1, zdir='z', color='#6F0382', linewidth=2.5)
  52.  
  53. # Plotting the 3 planes
  54. ax.plot_surface(x, y, f1(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='blue')
  55. ax.plot_surface(x, y, f2(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='red')
  56. ax.plot_surface(x, y, f3(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='green')
  57.  
  58. # Putting the limits in the axes
  59. ax.set_xlim(-10, 10)
  60. ax.set_ylim(-10, 10)
  61. ax.set_zlim(-15, 15)
  62.  
  63. # Writing the axis labels
  64. ax.set_xlabel('x', color='blue')
  65. ax.set_ylabel('y', color='blue')
  66. ax.set_zlabel('z', color='blue')
  67.  
  68. # Writing The Title of The Plot
  69. ax.set_title(r'$Graphical\; Resolution\; Linear\; System\; 3 \times 3$', fontsize=18)
  70.  
  71. # Stablishing the plots of our legend labels
  72. blue_proxy = plt.Rectangle((0, 0), 1, 1, fc='b')
  73. red_proxy = plt.Rectangle((0, 0), 1, 1, fc='r')
  74. green_proxy = plt.Rectangle((0, 0), 1, 1, fc='g')
  75. purple_proxy = plt.Rectangle((0, 0), 1, 1, fc='#6F0382')
  76.  
  77. # Drawing Our Legend
  78. ax.legend([blue_proxy,red_proxy, green_proxy, purple_proxy], [r'$x-2y+2z=0$',r'$2x-y+z=0$', r'$x-y+z=0$', r'$Sol.\; (0,\lambda, \lambda)$'], loc='upper left')
  79.  
  80. plt.show()

Una imagen del resultado final.

Creo que esta es una muy buena opción, mejor que los dibujos de muchos libros de texto. Y fácil de hacer, haces copy-paste de este mismo código, cambias las ecuaciones y arreando que es gerundio.


Saludos