Resolución Gráfica Sistema De Vandermonde Python 3

por | 6 febrero, 2014

Hola
En ésta ocasión he realizado un código en Python 3 utilizando Sympy y Matplotlib para solucionar analítica y gráficamente un sistema generalizado 3×3 de Vandermonde, que tiene la siguiente expresión:
\left\{\begin{array}{lr} x+y+z=1\\ ax+by+cz=d\\ a^2x+b^2y+c^2z=d^2\end{array}\right\}\;\forall\; a,b,c,d\in \mathbb{R}
Para solucionarlo he realizado un estudio de rangos con Sympy y luego para los casos en que la matriz de coeficientes y la ampliada son 2 ó 3 he programado la Regla de Cramer. Cuando ambos grados son 1 he insertado directamente en el código el conjunto de soluciones. Cuando los rangos son diferentes devuelvo como solución la cadena “Sistema Incompatible”. Las respectivas soluciones las imprimo por terminal.
Para que el programa tenga su utilidad en la docencia dibujo los planos y soluciones en 3D utilizando la librería Matplotlib, diferenciando casos respecto los rangos.
El código si estás acostumbrado a Python no tiene mayor complicación, un poco largo sólo.

  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.  
  20. # Declaring the three planes as functions
  21. f1 = lambda x, y: 1 - x - y
  22. f2 = lambda x, y, a, b, c, d: (d - a * x - b * y) / c
  23. f3 = lambda x, y, a, b, c, d: (d ** 2 - a ** 2 * x - b ** 2 * y) / c ** 2
  24.  
  25. # Declaring symbolic variables
  26. x, y, z = symbols('x y z')
  27.  
  28. # Declaring symbolic parameters
  29. a, b, c, d, alpha, beta = symbols('a b c d alpha beta')
  30.  
  31.  
  32. def sistema(a, b, c, d):
  33.     ''' Solve the Vandermonde System with Cramer Rule'''
  34.  
  35.     'Declaring coef matrix and sol matrices'
  36.     Acoef = Matrix([[1, 1, 1], [a, b, c], [a**2, b**2, c**2]])
  37.     Aampl = Matrix([[1, 1, 1, 1], [a, b, c, d], [a**2, b**2, c**2, d**2]])
  38.     rangAcoef = Acoef.rank()
  39.     rangAampl = Aampl.rank()
  40.  
  41.     print(('La solución del sistema es:'))
  42.  
  43.     if rangAcoef == 2 and rangAampl == 2 and a == b:
  44.         Ay = Matrix([[1-alpha, 1], [d-a*alpha, c]])
  45.         Az = Matrix([[1, 1-alpha], [b, d-a*alpha]])
  46.         Acoef = Matrix([[1, 1], [b, c]])
  47.         'Figuring out the determinants'
  48.         detcoef = factor(det(Acoef))
  49.         dety = factor(det(Ay))
  50.         detz = factor(det(Az))
  51.         'Figuring out the three coordinate solutions'
  52.         xsol = alpha
  53.         ysol = dety / detcoef
  54.         zsol = detz / detcoef
  55.         'The vector generalized solution'
  56.         solsist = [xsol, ysol, zsol]
  57.         print((solsist))
  58.         return solsist
  59.  
  60.     elif rangAcoef == 2 and rangAampl == 2 and a == c:
  61.         Ay = Matrix([[1-alpha, 1], [d-a*alpha, c]])
  62.         Az = Matrix([[1, 1-alpha], [b, d-a*alpha]])
  63.         Acoef = Matrix([[1, 1], [b, c]])
  64.         'Figuring out the determinants'
  65.         detcoef = factor(det(Acoef))
  66.         dety = factor(det(Ay))
  67.         detz = factor(det(Az))
  68.         'Figuring out the three coordinate solutions'
  69.         xsol = alpha
  70.         ysol = dety / detcoef
  71.         zsol = detz / detcoef
  72.         'The vector generalized solution'
  73.         solsist = [xsol, ysol, zsol]
  74.         print((solsist))
  75.         return solsist
  76.  
  77.     elif rangAcoef == 2 and rangAampl == 2 and b == c:
  78.         Ax = Matrix([[1-alpha, 1], [d-c*alpha, b]])
  79.         Ay = Matrix([[1, 1-alpha], [a, d-c*alpha]])
  80.         Acoef = Matrix([[1, 1], [a, b]])
  81.         'Figuring out the determinants'
  82.         detcoef = factor(det(Acoef))
  83.         detx = factor(det(Ax))
  84.         dety = factor(det(Ay))
  85.         'Figuring out the three coordinate solutions'
  86.         xsol = detx / detcoef
  87.         ysol = dety / detcoef
  88.         zsol = alpha
  89.         'The vector generalized solution'
  90.         solsist = [xsol, ysol, zsol]
  91.         print((solsist))
  92.         return solsist
  93.  
  94.     elif rangAcoef == 1 and rangAampl == 1:
  95.         xsol = 1 - alpha - beta
  96.         ysol = alpha
  97.         zsol = beta
  98.         solsist = [xsol, ysol, zsol]
  99.         print((solsist))
  100.         return solsist
  101.  
  102.     elif rangAcoef != rangAampl:
  103.         solsist = 'Sistema Incompatible'
  104.         print((solsist))
  105.         return solsist
  106.  
  107.     elif rangAcoef == 3 and rangAampl == 3:
  108.         Ax = Matrix([[1, 1, 1], [d, b, c], [d**2, b**2, c**2]])
  109.         Ay = Matrix([[1, 1, 1], [a, d, c], [a**2, d**2, c**2]])
  110.         Az = Matrix([[1, 1, 1], [a, b, d], [a**2, b**2, d**2]])
  111.         'Figuring out the determinants'
  112.         detcoef = factor(det(Acoef))
  113.         detx = factor(det(Ax))
  114.         dety = factor(det(Ay))
  115.         detz = factor(det(Az))
  116.  
  117.         'Figuring out the three coordinate solutions'
  118.         xsol = detx / detcoef
  119.         ysol = dety / detcoef
  120.         zsol = detz / detcoef
  121.  
  122.         'The vector generalized solution'
  123.         solsist = [xsol, ysol, zsol]
  124.  
  125.         print((solsist))
  126.         return solsist
  127.  
  128.  
  129. def dibuja(a, b, c, d, sol):
  130.     Acoef = Matrix([[1, 1, 1], [a, b, c], [a**2, b**2, c**2]])
  131.     Aampl = Matrix([[1, 1, 1, 1], [a, b, c, d], [a**2, b**2, c**2, d**2]])
  132.     rangAcoef = Acoef.rank()
  133.     rangAampl = Aampl.rank()
  134.  
  135.     # Stablishing our ranges for our variables
  136.     x1 = y1 = np.arange(-7, 7, 0.25)
  137.     ceros = np.zeros(len(x1))
  138.  
  139.     # Stablishing our meshgrid
  140.     x, y = np.meshgrid(x1, y1)
  141.  
  142.     # Our 3D Canvas Figure Plot
  143.     fig = plt.figure()
  144.     ax = fig.add_subplot(111, projection='3d')
  145.  
  146.     # Putting the limits in the axes
  147.     ax.set_xlim(-15, 15)
  148.     ax.set_ylim(-15, 15)
  149.     ax.set_zlim(-15, 15)
  150.  
  151.     # Writing the axis labels
  152.     ax.set_xlabel('x', color='blue')
  153.     ax.set_ylabel('y', color='blue')
  154.     ax.set_zlabel('z', color='blue')
  155.  
  156.     # Writing The Title of The Plot
  157.     ax.set_title(r'$Graphical\; Resolution\; Vandermonde\; System$', fontsize=18, loc='right')
  158.  
  159.     # Plotting the 3 planes
  160.     if rangAcoef == 1 and rangAampl == 1:
  161.         ax.plot_surface(x, y, f1(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='blue')
  162.         # Stablishing the plots of our legend labels
  163.         blue_proxy = plt.Rectangle((0, 0), 1, 1, fc='b')
  164.         # Drawing Our Legend
  165.         text1 = r'${}x+{}y+{}z={}$'.format(a, b, c, d)
  166.         textsol = r'$( {},\; {},\; {} ), \; \forall \; \alpha , \; \beta \; reales$'.format(latex(sol[0]), latex(sol[1]), latex(sol[2]))
  167.         ax.legend([blue_proxy, blue_proxy], [r'$x+y+z=1$', textsol], loc='upper left')
  168.  
  169.     elif rangAcoef == 2 and rangAampl == 2:
  170.         ax.plot_surface(x, y, f1(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='blue')
  171.         ax.plot_surface(x, y, f2(x, y, a, b, c, d), rstride=1, cstride=1, linewidth=0, antialiased=True, color='red')
  172.  
  173.         # Stablishing the plots of our legend labels
  174.         blue_proxy = plt.Rectangle((0, 0), 1, 1, fc='b')
  175.         red_proxy = plt.Rectangle((0, 0), 1, 1, fc='r')
  176.         black_proxy = plt.Line2D([0], [0], color='#6F0382', lw=2)
  177.         # Drawing Our Legend
  178.         text1 = r'${}x+{}y+{}z={}$'.format(a, b, c, d)
  179.         textsol = r'$( {},\; {},\; {} ), \; \forall \; \alpha \; real$'.format(latex(sol[0]), latex(sol[1]), latex(sol[2]))
  180.         ax.legend([blue_proxy, red_proxy, black_proxy], [r'$x+y+z=1$', text1, textsol], loc='upper left')
  181.  
  182.     elif rangAcoef != rangAampl:
  183.         ax.plot_surface(x, y, f1(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='blue')
  184.         ax.plot_surface(x, y, f2(x, y, a, b, c, d), rstride=1, cstride=1, linewidth=0, antialiased=True, color='red')
  185.         ax.plot_surface(x, y, f3(x, y, a, b, c, d), rstride=1, cstride=1, linewidth=0, antialiased=True, color='green')
  186.         # Stablishing the plots of our legend labels
  187.         blue_proxy = plt.Rectangle((0, 0), 1, 1, fc='b')
  188.         red_proxy = plt.Rectangle((0, 0), 1, 1, fc='r')
  189.         green_proxy = plt.Rectangle((0, 0), 1, 1, fc='g')
  190.         black_proxy = plt.Line2D([0], [0], linestyle="none", marker='o', alpha=0.6, markersize=10, markerfacecolor='black')
  191.  
  192.         # Drawing Our Legend
  193.         text1 = r'${}x+{}y+{}z={}$'.format(a, b, c, d)
  194.         text2 = r'${}x+{}y+{}z={}$'.format(a ** 2, b ** 2, c ** 2, d ** 2)
  195.         textsol = r'$Sistema\; Incompatible$'
  196.         ax.legend([blue_proxy, red_proxy, green_proxy, black_proxy], [r'$x+y+z=1$', text1, text2, textsol], numpoints=1, loc='upper left')
  197.  
  198.     elif rangAcoef == 3 and rangAampl == 3:
  199.         ax.plot_surface(x, y, f1(x, y), rstride=1, cstride=1, linewidth=0, antialiased=True, color='blue')
  200.         ax.plot_surface(x, y, f2(x, y, a, b, c, d), rstride=1, cstride=1, linewidth=0, antialiased=True, color='red')
  201.         ax.plot_surface(x, y, f3(x, y, a, b, c, d), rstride=1, cstride=1, linewidth=0, antialiased=True, color='green')
  202.         ax.plot([sol[0]], [sol[1]], [sol[2]], markerfacecolor='k', markeredgecolor='k', marker='o', markersize=5, alpha=0.6)
  203.         # Stablishing the plots of our legend labels
  204.         blue_proxy = plt.Rectangle((0, 0), 1, 1, fc='b')
  205.         red_proxy = plt.Rectangle((0, 0), 1, 1, fc='r')
  206.         green_proxy = plt.Rectangle((0, 0), 1, 1, fc='g')
  207.         black_proxy = plt.Line2D([0], [0], linestyle="none", marker='o', alpha=0.6, markersize=10, markerfacecolor='black')
  208.  
  209.         # Drawing Our Legend
  210.         text1 = r'${}x+{}y+{}z={}$'.format(a, b, c, d)
  211.         text2 = r'${}x+{}y+{}z={}$'.format(a ** 2, b ** 2, c ** 2, d ** 2)
  212.         textsol = r'$Sol.\; ({},{},{})$'.format(sol[0], sol[1], sol[2])
  213.         ax.legend([blue_proxy,red_proxy, green_proxy, black_proxy], [r'$x+y+z=1$', text1, text2, textsol], numpoints=1, loc='upper left')
  214.  
  215.     plt.show()
  216.  
  217.  
  218. #COMP. INDETDO. RANGO 1
  219. #sol = sistema(1, 1, 1, 1)
  220. #dibuja(1, 1, 1, 1, sol)
  221.  
  222. #COMP. INDETDO. RANGO 2 a = b
  223. #sol = sistema(1, 1, 2, 1)
  224. sol = sistema(1/2, 1/2, -4, 1/2)
  225. dibuja(1/2, 1/2, -4, 1/2, sol)
  226.  
  227. #COMP. INDETDO. RANGO 2 a = c
  228. #sol = sistema(1, 2, 1, 1)
  229. #dibuja(1, 2, 1, 1, sol)
  230.  
  231. #COMP. INDETDO. RANGO 2 b = c
  232. #sol = sistema(2, 1, 1, 1)
  233. #dibuja(2, 1, 1, 1, sol)
  234.  
  235. #INCOMPATIBLE
  236. #sol = sistema(1, 1, 1, 2)
  237. #dibuja(1, 1, 1, 2, sol=None)
  238.  
  239. #COMPATIBLE DETERMINADO
  240. #sol = sistema(1, 2, 3, 1)
  241. #sol = sistema(1, 2, 3, 1)
  242. #dibuja(1, 2, 3, 1, sol)

Una imagen para un caso indeterminado:

Puedes descargar el código en:


Ahora con éste código podemos mostrar a nuestros alumnos de forma rápida y eficiente qué ocurre con éste conjunto de sistemas lineales de Vandermonde según los valores de a,b,c,d reales.
Desde luego resolver problemas matemáticos en Python 3 utilizando librerías como Sympy, Numpy, SciPy o Matplotlib es relativamente sencillo y útil, es como trabajar con Octave, Matplotlib, Mathematica y Scilab todo junto pero teniendo yo el mango por la sartén.
Saludos

Un pensamiento en “Resolución Gráfica Sistema De Vandermonde Python 3

  1. Pingback: Bitacoras.com

Los comentarios están cerrados.