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

Como punto de partida de la unidad se presenta una solución a un sistema de 3 ecuaciones y 3 incógnitas. El ejercicio se presenta desde la perspectiva que la solución es un punto de intersección de los planos en el espacio dados por cada ecuación.

Referencia: Chapra 5Ed 9.1 p247, pdf 271.
Ejercicio: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

Considere el siguiente sistema de ecuaciones:

\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}

Se pueden representar como planos en el espacio despejando la variable z para cada ecuación, de la forma:

\begin{cases} z=(1+ 2x - 5y)/9\\z=(6 -7x-y)/1\\z=(-26+3x-7y)/(-1)\end{cases}

Ecuaciones como Planos en el espacio

Si observamos la primera ecuación en una gráfica, encontramos que para graficar el plano se requiere definir los intervalos para los ejes X, Y, por ejemplo:

z=\frac{1+ 2x - 5y}{9} -5 \leq x \leq 5 -7 \leq y \leq 7

Al incorporar el plano de la segunda ecuación encontramos que la solución al sistema es la recta interseción de ambos planos.

Usando la tercera ecuación, la intersección de los planos genera un punto cuyas coordenadas corresponden a la solución del sistema.


Algoritmo en Python

Para observar mejor del resultado en tres dimensiones, se plantean algunas intrucciones básicas para graficar y con el resultado obtener una mejor vista rotando el gráfico con el cursor.

Las ecuaciones del problema se pueden escribir en forma matricial Ax=B, usando solamente los coeficientes.

\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}
# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

B = np.array([1,6,-26])

Con la forma matricial se generan las ecuaciones de cada plano usando los coeficientes al despejar la variable z.

Para visualizar los planos generados por cada ecuación se escribe cada ecuación como funciones(x,y) de dos variables. Para simplificar la escritura se usa el formato lambda.

# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]

Intervalos de observación

Para la gráfica se requiere usar intervalos y muestras para las variables independientes x,y de forma semejante a las gráficas en 2D.

Las muestras en cada eje generan puntos en el plano X,Y y luego se usan todas las combinaciones que se generan.

– Primero se generan las muestras para cada eje en los vectores xi, yi.

ax = -5     # Intervalos
bx = 5
ay = ax-2
by = bx+2
muestras = 21

xi = np.linspace(ax,bx, muestras)
yi = np.linspace(ay,by, muestras)

– Las combinaciones entre las muestras de cada eje se obtienen generando las matrices Xi, Yi que representan la malla de muestras en el plano para evaluar cada punto.

Xi, Yi = np.meshgrid(xi,yi)

– Se evaluan los puntos Xi,Yi en cada ecuación, generando las matrices Zi

Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)

Si es necesario, revise la sección de Gráficas 3D wireframe en la sección de Recursos/Resumen Python.

Se grafica cada ecuación como una «malla de alambre» o wireframe, usando librerias 3D, usando Xi,Yi, Z0, luego con Z1 y Z2.

Por simplicidad del algoritmo, pues la revisión del concepto es sobre intersección de planos, la solución se obtiene con las funciones de numpy.

punto = np.linalg.solve(A,B)

El punto se incorpora a la gráfica como un punto usando «dispersión» o «scatter».

El algoritmo completo se muestra a continuación, e incluye la parte gráfica:

# 1ra Evaluación II Término 2011
# Tema 2. Sistema de ecuaciones 3x3
# Concepto como interseccion de Planos

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

B = np.array([1,6,-26])

ax = -5     # Intervalo X
bx = 5
ay = ax-2   # Intervalo Y
by = bx+2

muestras = 11

# PROCEDIMIENTO --------
# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]

xi = np.linspace(ax,bx, muestras)
yi = np.linspace(ay,by, muestras)
Xi, Yi = np.meshgrid(xi,yi)

Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)

# solución al sistema
punto = np.linalg.solve(A,B)

# SALIDA
print('respuesta de A.X=B : ')
print(punto)

# Interseccion entre ecuacion 1 y 2
# PlanoXZ, extremo inferior de y
Aa  = np.copy(A[0:2,[0,2]])
Ba  = np.copy(B[0:2])
Ba  = Ba-ay*A[0:2,1]
pta = np.linalg.solve(Aa,Ba)
pa  = np.array([ay])
pxza = np.array([pta[0],ay,pta[1]])

# PlanoXZ, extremo superior de y
Ba  = Ba-by*A[0:2,1]
ptb = np.linalg.solve(Aa,Ba)
pb  = np.array([by])
pxzb = np.array([ptb[0],by,ptb[1]])

# GRAFICA de planos
figura = plt.figure()
grafica = figura.add_subplot(111, projection='3d')

grafica.plot_wireframe(Xi,Yi,Z0,
                       color ='blue',
                       label='Ecuación 1')
grafica.plot_wireframe(Xi,Yi,Z1,
                       color ='green',
                       label='Ecuación 2')
grafica.plot_wireframe(Xi,Yi,Z2,
                       color ='orange',
                       label='Ecuación 3')
# recta intersección planos 1 y 2
grafica.plot([pxza[0],pxzb[0]],
             [pxza[1],pxzb[1]],
             [pxza[2],pxzb[2]],
             label='Sol 1y2',
             color = 'violet',
             linewidth = 4)
# Punto solución del sistema 3x3
grafica.scatter(punto[0],punto[1],punto[2],
                color = 'red',
                marker='o',
                label ='punto',
                linewidth = 6)

grafica.set_title('Sistema de ecuaciones 3x3')
grafica.set_xlabel('x')
grafica.set_ylabel('y')
grafica.set_zlabel('z')
grafica.legend()
grafica.view_init(45, 45)
# rotacion de ejes
for angulo in range(45, 360+45, 5 ):
    grafica.view_init(45, angulo)
    plt.draw()
    plt.pause(.001)
plt.show()