3.6 Método de Jacobi con Python

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


1. El método iterativo de Jacobi

Referencia: Burden 7.3 p334, Rodríguez 5.1 p154, Chapra  11.2.1p312

La analogía presentadas entre la «norma como distancia 3D» y el «error de acoplamiento de aeronaves»,  

permite considerar desde un punto de partida o inicial las aproximaciones o iteraciones sucesivas hacia una solución del sistema de ecuaciones.

Las iteraciones pueden ser convergentes o no.

Los métodos iterativos para sistemas de ecuaciones, son semejantes al método de punto fijo para búsqueda de raíces, requieren un punto inicial para la búsqueda de la raíz o solución que satisface el sistema.

Para describir el método iterativo de Gauss-Seidel, se usa un sistema de 3 incógnitas y 3 ecuaciones, diagonalmente dominante.

\begin{cases} a_{0,0} x_0+a_{0,1}x_1+a_{0,2} x_2 = b_{0} \\ a_{1,0} x_0+a_{1,1}x_1+a_{1,2} x_2 = b_{1} \\ a_{2,0} x_0+a_{2,1}x_1+a_{2,2} x_2 = b_{2} \end{cases}

Para facilitar la escritura del algoritmo, note el uso de índices ajustado a la descripción de arreglos en Python para la primera fila i=0 y primera columna j=0.

Semejante a despejar una variable de la ecuación para representar un plano, se plantea despejar una variable de cada ecuación. Se obtiene así los valores de cada xi, por cada por cada fila i:

x_0 = \frac{b_{0} -a_{0,1}x_1 -a_{0,2} x_2 }{a_{0,0}} x_1 = \frac{b_{1} -a_{1,0} x_0 -a_{1,2} x_2}{a_{1,1}} x_2 = \frac{b_{2} -a_{2,0} x_0 - a_{2,1} x_1}{a_{2,2}}

Observe que el patrón de cada fórmula para determinar x[i], tiene la forma:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}

La parte de la sumatoria se realiza para cada término de A[i,j] en la fila i, excepto para el término de la diagonal A[i,i].

Si se tiene conocimiento del problema planteado y se puede «intuir» o «suponer» una solución para el vector X. Por ejemplo, iniciando con el vector cero, es posible calcular un nuevo vector X usando las ecuaciones para cada X[i] encontradas.

X = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}

Con cada nuevo valor actualiza el el vector X, en  Xnuevo, se calcula el vector diferencias entre X y el vector Xnuevo .

El error a prestar la atención es al mayor valor de las diferencias; se toma como condición para repetir la evaluación de cada vector.

     nuevo = [ 0, 0,  0]
         X = [x0, x1, x2]
diferencia = [|0 - x0|, |0 - x1|, |0 - x2|]
errado = max(diferencia)

Se observa los resultados de errado para cada iteración, relacionados con la convergencia. Si luego de «muchas» iteraciones se encuentra que (errado>tolera),  se detiene el proceso.

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


2. Ejercicio

Referencia:  Steven C. Chapra. Numerical Methods 7th Edition,  Ejercicio 12.39 p339; 1Eva_IT2018_T3 Temperatura en nodos de placa.

La temperatura en los nodos de la malla de una placa se puede calcular con el promedio de las temperaturas de los 4 nodos vecinos de la izquierda, derecha, arriba y abajo. Placa Temperatura

Una placa cuadrada de 3 m de lado tiene la temperatura en los nodos de los bordes como se indica en la figura,

a) Plantee el sistema de ecuaciones y resuelva con el método de Jacobi  para encontrar los valores de a, b, c, d.

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


3. Desarrollo Analítico

Solución Propuesta: s1Eva_IT2018_T3 Temperatura en nodos de placa

a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior:

a=\frac{50+c+100+b}{4} b=\frac{a+30+50+d}{4} c=\frac{a+60+100+d}{4} d=\frac{b+60+c+30}{4}

que reordenando se convierte en:

4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90

simplificando:

4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90

que a forma matricial se convierte en:

A = [[ 4, -1, -1, 0],
     [ 1, -4,  0, 1],
     [ 1,  0, -4, 1],
     [ 0,  1,  1,-4]]
B = [150, -80,-160,-90]

Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas.  Se aumentó el punto decimal a los valores de la matriz A y el vector B  para que sean considerados como números reales.

El número de condición es: np.linalg.cond(A) = 3.0

que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es «estable» y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.

b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial

X0 = [60,40,70,50] 

reemplazando los valores iniciales en cada ecuación sin cambios.

iteración 1
a=\frac{50+70+100+40}{4} = 65

b=\frac{60+30+50+50}{4} = 47.5 c=\frac{60+60+100+50}{4} = 67.5 d=\frac{40+60+70+30}{4} = 50
X1 = [65.  47.5 67.5 50. ]

vector de error = 
     [|65-60|,|47.5-40|, |67.5-70|, |50-50|]
  = [|5|,|7.5|,|-2.5|, |0|]
errormax = 7.5

iteración 2
a=\frac{50+67.5+100+47.5}{4} = 66.25

b=\frac{65+30+50+50}{4} = 48.75 c=\frac{65+60+100+50}{4} = 68.75 d=\frac{47.5+60+67.5+30}{4} = 51.25
X(2) = [66.25 48.75 68.75 51.25]

vector de error = 
       [|66.25-65|,|48.75-47.5|,
        |68.75-67.5|,|51.25-50|] 
  = [|1.25|,|1.25|,|1.25|,|1.25|]
errormax = 1.25

iteración 3
a=\frac{50+68.75+100+48.75}{4} = 66.875

b=\frac{66.25+30+50+51.25}{4} = 49.375 c=\frac{66.25+60+100+51.25}{4} = 69.375 d=\frac{48.75+60+68.75+30}{4} = 51.875
X(3) = [66.875 49.375 69.375 51.875]
vector de error = 
      [|66.875-66.25|,|49.38-48.75|,
       |69.3875-68.75|,|51.875-51.25|]
    = [|0.625|,|0,63|, |0.6375|,|0.625|]
errormax = 0.625

siguiendo las iteraciones se debería llegar a:

>>> np.linalg.solve(A,B)
array([67.5, 50. , 70. , 52.5])

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..


4. Algoritmo en Python

Tarea: Realice el algoritmo de Jacobi en Python, incluyendo el pivoteo por filas de la matriz.

Resultados con el algoritmo

itera,[X],errado
0 [60. 40. 70. 50.] 1.0
1 [65.  47.5 67.5 50. ] 7.5
2 [66.25 48.75 68.75 51.25] 1.25
3 [66.875 49.375 69.375 51.875] 0.625
4 [67.1875 49.6875 69.6875 52.1875] 0.3125
5 [67.34375 49.84375 69.84375 52.34375] 0.15625
6 [67.421875 49.921875 69.921875 52.421875] 0.078125
7 [67.4609375 49.9609375 69.9609375 52.4609375] 0.0390625
8 [67.48046875 49.98046875 69.98046875 52.48046875] 0.01953125
9 [67.49023438 49.99023438 69.99023438 52.49023438] 0.009765625
10 [67.49511719 49.99511719 69.99511719 52.49511719] 0.0048828125
11 [67.49755859 49.99755859 69.99755859 52.49755859] 0.00244140625
12 [67.4987793 49.9987793 69.9987793 52.4987793] 0.001220703125
13 [67.49938965 49.99938965 69.99938965 52.49938965] 0.0006103515625
14 [67.49969482 49.99969482 69.99969482 52.49969482] 0.00030517578125
15 [67.49984741 49.99984741 69.99984741 52.49984741] 0.000152587890625
16 [67.49992371 49.99992371 69.99992371 52.49992371] 7.62939453125e-05
numero de condición: 3.000000000000001
respuesta con Jacobi
[[67.49992371]
 [49.99992371]
 [69.99992371]
 [52.49992371]]

 verificar A.X:
[[ 149.99984741]
 [ -79.99984741]
 [-159.99984741]
 [ -89.99984741]]
max|A.X - B|
0.000152587890625
>>> 

Instrucciones en Python

# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada
# Método de Jacobi
import numpy as np

def jacobi(A,B,tolera,X0,iteramax=100, vertabla=False):
    ''' Método de Jacobi, tolerancia, vector inicial X0
        para mostrar iteraciones: vertabla=True
    '''
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)
    xnuevo = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('itera,[X],errado')
        print(itera, xnuevo, errado)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        X = np.copy(xnuevo)
        itera = itera + 1
        if vertabla==True:
            print(itera, xnuevo, errado)
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)


# INGRESO
A = np.array([[ 4, -1, -1,  0],
              [ 1, -4,  0,  1],
              [ 1,  0, -4,  1],
              [ 0,  1,  1, -4]],dtype=float)
B = np.array([150, -80,-160,-90],dtype=float)
X0  = np.array([60,40,70,50],dtype=float)

tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

respuesta = jacobi(A,B,tolera,X0,vertabla=True)

verifica = np.dot(A,respuesta)
verificaErrado = np.max(np.abs(B-np.transpose(verifica)))

# SALIDA
print('numero de condición:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('\n verificar A.X:')
print(verifica)
print('max|A.X - B|')
print(verificaErrado)

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]