9.2 EDP Elípticas

Referencia: Chapra 29.1 p866 pdf890, Rodriguez 10.3 p425, Burden 12.1 p694 pdf704

Las Ecuaciones Diferenciales Parciales tipo elípticas semejantes a la mostrada:

\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{ \delta y^2} = 0

(ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)

se interpreta como una placa metálica de dimensiones Lx y Ly, delgada con aislante que recubren las caras de la placa, y sometidas a condiciones en las fronteras:

Lx = dimensión x de placa metálica
Ly = dimensión y de placa metálica
u[0,y]  = Ta
u[Lx,y] = Tb
u[x,0]  = Tc
u[x,Ly] = Td

Para el planteamiento se usa una malla en la que cada nodo corresponden a los valores u[xi,yj]. Para simplificar la nomenclatura se usan los subíndices i para el eje de las x y j para el eje t, quedando u[i,j].

Se discretiza la ecuación usando diferencias divididas que se sustituyen en la ecuación, por ejemplo:

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^2}=0

Se agrupan los términos de los diferenciales:

\frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0

con lo que se simplifican los valores como uno solo \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1 . Por facilidad de lo que se realiza se supone que lambda tiene valor de 1 o los delta son iguales.

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

Obteniendo así la solución numérica conceptual. La forma de resolver el problema determina el nombre del método a seguir.

9.1 EDP Parabólicas

Referencia:  Chapra 30.2 p.888 pdf.912, Burden 9Ed p714, Rodriguez 10.2 p.406

Las Ecuaciones Diferenciales Parciales tipo parabólicas semejantes a la mostrada, representa la ecuación de calor para una barra aislada sometida a fuentes de calor en cada extremo.

La temperatura se representa en el ejercicio como u[x,t]

\frac{d ^2 u}{dx ^2} = K\frac{d u}{d t}

Para la solución numérica, se discretiza la ecuación usando diferencias finitas divididas que se sustituyen en la ecuación,

\frac{d^2 u}{dx^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta X)^2} \frac{du}{dt} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t}

con lo que la ecuación continua se convierte a discreta:

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = K\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

Para interpretar mejor el resultado, se usa una malla que en cada nodo representa la temperatura como los valores u[xi,tj].
Para simplificar nomenclatura se usan los subíndices i para el eje de las x y j para el eje t, quedando u[ i , j ].

En el enunciado del problema habían establecido los valores en las fronteras:
– temperaturas en los extremos Ta, Tb
– la temperatura inicial de la barra T0,
– El parámetro para la barra K.

El resultado obtenido se interpreta como los valores de temperatura a lo largo de la barra luego de transcurrido un largo tiempo. Las temperaturas en los extremos de la barra varían entre Ta y Tb a lo largo del tiempo.

Tomando como referencia la malla, existirían algunas formas de plantear la solución, dependiendo de la diferencia finita usada: centrada, hacia adelante, hacia atrás.


Resultado en una animación con la variable tiempo:


Tarea: Revisar ecuación para difusión de gases, segunda ley de Fick.

La difusión molecular desde un punto de vista microscópico y macroscópico.

Sistemas EDO. modelo predador-presa

Referencia: Chapra 28.2 p831 pdf855, Rodriguez 9.2.1 p263
https://es.wikipedia.org/wiki/Ecuaciones_Lotka%E2%80%93Volterra
https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations

Modelos depredador- presa y caos. Ecuaciones Lotka-Volterra. En el sistema de ecuaciones:

\frac{dx}{dt} = ax - bxy \frac{dy}{dt} = -cy + dxy

x =  número de presas
y = depredadores
a = razón de crecimiento de la presa, (0.5)
c = razón de muerte del depredador (0.35)
b = efecto de la interacción depredador-presa sobre la muerte de la presa (0.7)
d = efecto de la interacción depredador-presa sobre el crecimiento del depredador, (0.35)

Los términos que multiplican xy hacen que las ecuaciones sean no lineales.


Para resolver el sistema, se plantean las ecuaciones de forma simplificada para el algoritmo:

f =  a x - b x y
g = -c y + d x y
con valores iniciales:
t = 0,x = 2, y = 1

Planteamiento que se ingresan al algoritmo con el algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en:

8.2.2 Runge-Kutta d2y/dx2

http://blog.espol.edu.ec/matg1013/8-2-2-runge-kutta-d2y-dx2/

Al ejecutar el algoritmo se obtienen los siguientes resultados:

 [ ti, xi, yi] 
[[  0.         2.         1.      ]
 [  0.5        1.754875   1.16975 ]
 [  1.         1.457533   1.302069]
 [  1.5        1.167405   1.373599]
 [  2.         0.922773   1.381103]
 [  2.5        0.734853   1.33689 ]
 [  3.         0.598406   1.258434]
 ...
 [ 49.         1.11309    1.389894]
 [ 49.5        0.876914   1.38503 ]
 [ 50.         0.698717   1.331107]
 [ 50.5        0.570884   1.246132]]
>>> 

Los resultados se pueden observar de diferentes formas:
a) Cada variable xi, yi versus ti, es decir cantidad de animales de cada especie durante el tiempo de observación

b) Independiente de la unidad de tiempo, xi vs yi, muestra la relación entre la cantidad de presas y predadores. Relación que es cíclica y da la forma a la gráfica.

Las instrucciones del algoritmo en Python usadas en el problema son:

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,3),dtype=float)
    tabla[0] = [t0,x0,y0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi]
    tabla = np.array(tabla)
    return(tabla)

# Programa
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35
# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y
# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1
# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi]')
print(tabla)

Los resultados numéricos se usan para generar las gráficas presentadas, añadiendo las instrucciones:

# Grafica tiempos vs población
import matplotlib.pyplot as plt
plt.plot(ti,xi, label='xi presa')
plt.plot(ti,yi, label='yi predador')
plt.title('Modelo predador-presa')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()
# gráfica xi vs yi
plt.plot(xi,yi)
plt.title('Modelo presa-predador [xi,yi]')
plt.xlabel('x presa')
plt.ylabel('y predador')
plt.grid()
plt.show()

Tarea: Añadir la animación de la gráfica usando la variable tiempo para mostrar un punto que marque el par ordenado[x,y] al variar t.

8.2.2 Runge-Kutta d2y/dx2

Para sistemas de ecuaciones diferenciales ordinarias con condiciones de inicio en x0, y0, y’0

\frac{\delta ^2 y}{\delta x^2} = \frac{\delta y}{\delta x} + etc

Forma estandarizada de la ecuación:

y'' = y' + etc

se convierte a:

y' = z = f_x(x,y,z) z' = (y')' = z + etc = g_x(x,y,z)

con las condiciones de inicio en x0, y0, z0

y se pueden reutilizar los métodos para primeras derivadas, por ejemplo Runge-Kutta de 2do y 4to orden.

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)
Runge-Kutta 4do Orden tiene error de truncamiento O(h5)

Runge-Kutta 2do Orden d2y/dx2 en Python:

import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0,z0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

Runge-Kutta 4do Orden d2y/dx2 en Python:

def rungekutta4_fg(fx,gx,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
        K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
        
        K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
        K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)

        K4y = h * fx(xi+h, yi + K3y, zi + K3z)
        K4z = h * gx(xi+h, yi + K3y, zi + K3z)

        yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
        zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

Una aplicación del algoritmo en Señales y Sistemas:

http://blog.espol.edu.ec/telg1001/01-respuesta-entrada-cero-ltic/

8.2.1 Runge-Kutta 4to Orden dy/dx

Referencia: Chapra 25.3.3 p746 pdf 770, Rodriguez 9.1.8 p358

Para una ecuación diferencial de primer orden con una condición de inicio, la fórmula de Runge-Kutta de 4to orden se obtiene de la expresión con cinco términos:

y_{i+1} = y_i + aK_1 + bK_2 + cK_3 + dK_4

siendo:

y'(x) = f(x_i,y_i) y(x_0) = y_0

debe ser equivalente a la serie de Taylor de 5 términos:

y_{i+1} = y_i + h f(x_i,y_i) + + \frac{h^2}{2!} f'(x_i,y_i) + \frac{h^3}{3!} f''(x_i,y_i) + +\frac{h^4}{4!} f'''(x_i,y_i) + O(h^5) x_{i+1} = x_i + h

que usando aproximaciones de derivadas, se obtienen:

def rungekutta4(d1y,x0,y0,h,muestras):
    # Runge Kutta de 4do orden
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h/2, yi + K1/2)
        K3 = h * d1y(xi+h/2, yi + K2/2)
        K4 = h * d1y(xi+h, yi + K3)

        yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
        xi = xi + h
        
        estimado[i] = [xi,yi]
    return(estimado)

Note que el método de Runge-Kutta de 4to orden es similar a la regla de Simpson 1/3. La ecuación representa un promedio ponderado para establecer la mejor pendiente.

Diferenciación numérica

Referencia: Chapra 23.1 p668 pdf692, Rodriguez 8.2 p324

La fórmula de Taylor muestra una aproximación de una función f(x):

P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 + ...

Primera Derivada con Diferencias divididas hacia adelante

Usando un polinomio de Taylor alrededor de xi en para un punto a la derecha xi+1 a una distancia h = xi+1xi

f(x_{i+1}) = f(x_i)+\frac{f'(x_i)}{1!} (h) + \frac{f''(x_i)}{2!}(h)^2 + ...

se puede simplificar en un polinomio de grado uno y un término de error:

f_{i+1} = f_i + f'_i (h) + O(h^2) ...
f'_i = \frac{f_{i+1}-f_i}{h} = \frac{\Delta f_i}{h}

La expresión es la primera diferencia finita dividida con un error del orden O(h).
Revise que el término de la expresión queda O(h2)/h con lo que se disminuye el exponente en uno.

Primera derivada con diferencias divididas centradas

Se realiza el mismo procedimiento que el anterior, usando un punto xi+1 y xi-1 alrededor de xi. En el término xi-1 el valor de h es negativo al invertir el orden de la resta.

f_{i+1} = f_i+\frac{f'_i}{1!}(h) + \frac{f''_i}{2!}(h)^2 f_{i-1} = f_i-\frac{f'_i}{1!}(h) + \frac{f''_i}{2!}(h)^2

restando la ecuaciones se tiene que

f_{i+1} - f_{i-1} = f'_i (h)+f'_i (h) f_{i+1} - f_{i-1} = 2h f'_i
f'_i = \frac{f_{i+1} - f_{i-1}}{2h}

con un error del orden O(h2)

Segunda derivadas

Al continuar con el procedimiento mostrado se pueden obtener las fórmulas para segundas derivadas, las que se resumen en las tablas de Fórmulas de diferenciación por diferencias divididas.

7.4 Cuadratura de Gauss

Referencia: Chapra 22.3 p655 pdf679, Rodriguez 7.3 p294, Burden 4.7 p220 pdf230

La cuadratura de Gauss aproxima el integral de una función en un intervalo [a,b] centrado en cero mediante un cálculo numérico con menos operaciones y evaluaciones de la función. Se representa como una suma ponderada:

I \cong c_0f(x_0) + c_1f(x_1)

para la fórmula de dos puntos se tiene obtiene:

c_0 = c_1 = 1 x_0 = -\frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}

para un intervalo de evaluación desplazado en el eje x se requiere convertir los puntos al nuevo rango. Se desplaza el punto cero al centro del intervalo [a,b] y se obtiene:

x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{b-a}{2}(f(x_a) + f(x_b))

cuya fórmula es semejante a una mejor aproximación de un trapecio, cuyos promedios de alturas son puntos internos de [a,b], concepto mostrado en la gráfica.


Algoritmo

El algoritmo se desarrolla en un tramo en el intervalo [a,b] junto a la gráfica para mostar el concepto. Para el ejemplo el integral buscado es:

\int_1^3 \sqrt{x} \sin(x) dx

que usando cuadratura de Gauss con un solo intervalo[a,b] tiene como resultado:

Integral:  2.03821975687
>>>

el resultado se puede mejorar aumentando el número de tramos en el intervalo [a,b]. Por ejemplo, el resultado usando 4 tramos el resultado es semejante al usar el método del trapecio con 128 tramos, lo que muestra el ahorro en calculos entre los métodos

Integral:  2.05357719003
>>> 

Las intrucciones en Python son:

# Integración: Cuadratura de Gauss de dos puntos
# modelo con varios tramos entre [a,b]
import numpy as np
import matplotlib.pyplot as plt

# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    return(area)

#Función a evaluar
def funcionx(x):
    fxi = np.sqrt(x)*np.sin(x)
    return(fxi)

# PROGRAMA
a = 1 # float(input('a: '))
b = 3 # float(input('b: '))
tramos = 4 # int(input('tramos: '))

# PROCEDIMIENTO
muestras = tramos+1
xi = np.linspace(a,b,muestras)
area = 0
for i in range(0,muestras-1,1):
    deltaA = integraCuadGauss2p(funcionx,xi[i],xi[i+1])
    area = area + deltaA
# SALIDA
print('Integral: ', area)

Gráfica por tramos

La gráfica que complementa el resultado anterior se realiza añadiendo las intrucciones presentadas a continuación.

Considere que la gráfica es útil con pocos tramos en el intervalo[a,b]

# Gráfica por cada subtramo
x0 = -1/np.sqrt(3) 
x1 = 1/np.sqrt(3)
muestrastramo = 20
subtramo = np.linspace(a,b,tramos+1)
xi = np.array([])
fi = np.array([])
xat = np.array([])
xbt = np.array([])
recta = np.array([])
for i in range(0,tramos,1):
    at = subtramo[i]
    bt = subtramo[i+1]
    xit = np.linspace(at,bt,muestrastramo)
    fit = funcionx(xit)
    xi = np.concatenate((xi,xit))
    fi = np.concatenate((fi,fit))
    # puntos xa y xb por tramo
    xa = (bt+at)/2 + (bt-at)/2*(x0)
    xb = (bt+at)/2 + (bt-at)/2*(x1)
    xat = np.concatenate((xat,[xa]))
    xbt = np.concatenate((xbt,[xb]))
    # Recta entre puntos x0 y x1 por tramo
    pendiente = (funcionx(xb)-funcionx(xa))/(xb-xa)
    b0 = funcionx(xa)-pendiente*xa
    linea = b0 + pendiente*xit
    recta = np.concatenate((recta,linea))
# Marcadores 'o' de xa y xb por tramos
puntox = np.concatenate((xat,xbt))
puntoy = funcionx(puntox)
    
# Trazado de lineas
plt.plot(xi,recta, label = 'grado 1', color = 'tab:orange')
plt.fill_between(xi,0,recta, color='tab:olive')
plt.plot(xi,fi, label='f(x)', color = 'b')
# Verticales para dividir los tramos
for j in range(0,len(subtramo),1):
    plt.axvline(subtramo[j], color='tab:gray')
# Marcadores de puntos xa y xb por tramos
for j in range(0,len(xat),1):
    plt.axvline(xat[j], color='w')
    plt.axvline(xbt[j], color='w')
plt.plot(puntox,puntoy, 'o', color='g')
plt.title('Integral: Cuadratura Gauss')
plt.xlabel('x')
plt.legend()
plt.show()

3.4 Gauss-Jordan Método

Referencia: Chapra 9.7 p277 pdf301

Para la solución con Gauss-Jordan se procede con la matriz aumentada para aplicar:

  • eliminación hacia adelante de incógnitas
  • eliminación hacia atras de incógnitas

sigue el ejemplo01 para matrices, usando la forma de función para el algoritmo.

A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

B = np.array([[60.70],
              [92.90],
              [56.30]])

Se crea la matriz aumentada AB.
Para el procedo de eliminación hacia adelante, en cada fila el elemento diagonal será el pivote, con el que se crea el coeficiente para eliminación de términos hacia adelante.

Para la eliminación hacia atras, se normaliza la fila haciendo el elemento diagonal uno, para calcular el coeficiente.

matriz aumentada:
[[  4.    2.    5.   60.7]
 [  2.    5.    8.   92.9]
 [  5.    4.    3.   56.3]]
 *** Gauss elimina hacia adelante ***
coeficiente:  2.0
[[   4.     2.     5.    60.7]
 [   0.     8.    11.   125.1]
 [   5.     4.     3.    56.3]]
coeficiente:  0.8
[[   4.      2.      5.     60.7 ]
 [   0.      8.     11.    125.1 ]
 [   0.      1.2    -2.6   -15.66]]
coeficiente:  6.66666666667
[[   4.            2.            5.           60.7       ]
 [   0.            8.           11.          125.1       ]
 [   0.            0.          -28.33333333 -229.5       ]]
 *** Gauss-Jordan elimina hacia atras *** 
coeficiente:  0.0909090909091
[[  4.           2.           5.          60.7       ]
 [  0.           0.72727273   0.           3.27272727]
 [ -0.          -0.           1.           8.1       ]]
coeficiente:  0.2
[[ 0.8         0.4         0.          4.04      ]
 [ 0.          0.72727273  0.          3.27272727]
 [-0.         -0.          1.          8.1       ]]
coeficiente:  2.5
[[ 2.   0.   0.   5.6]
 [ 0.   1.   0.   4.5]
 [-0.  -0.   1.   8.1]]
aumentada: 
[[ 1.   0.   0.   2.8]
 [ 0.   1.   0.   4.5]
 [-0.  -0.   1.   8.1]]
X: 
[[ 2.8]
 [ 4.5]
 [ 8.1]]
>>> 

El algoritmo desarrollado por partes:

# Gauss Jordan
# Supone que los elementos de diagonal no son cero
# Tarea: verificar diagonal no tiene ceros
# Muestra proceso
import numpy as np

# INGRESO
A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

B = np.array([[60.70],
              [92.90],
              [56.30]])

# PROCEDIMIENTO
# Matriz aumentada
AB = np.concatenate((A,B), axis=1)
print('matriz aumentada:')
print(AB)

print(' *** Gauss elimina hacia adelante ***')
casicero = 0 # 1e-15
# Gauss elimina hacia adelante
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(AB[k,i])>=casicero):
            coeficiente = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*coeficiente - AB[i,:]
        else:
            coeficiente= 'division para cero'
        print('coeficiente: ',coeficiente)
        print(AB)

print(' *** Gauss-Jordan elimina hacia atras *** ')
# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            coeficiente = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*coeficiente - AB[i,:]
        else:
            coeficiente= 'division para cero'
        print('coeficiente: ', coeficiente)
        print(AB)
X = AB[:,ultcolumna]
X = np.transpose([X])

# SALIDA
print('aumentada: ')
print(AB)
print('X: ')
print(X)

Tarea: implementar caso cuando aparecen ceros en la diagonal para dar respuesta, convertir a funciones cada parte

4.1.2 Normas funciones

Referencia: Chapra 10.3, p299, pdf323

Algunas normas vectoriales y matriciales. Cálculo del número de condición.

# Normas vectoriales y matriciales
# Referencia: Chapra 10.3, p299, pdf323
import numpy as np

def norma_p(X,p):
    Xp = (np.abs(X))**p
    suma = np.sum(Xp)
    norma = suma**(1/p)
    return(norma)

def norma_euclidiana(X):
    norma = norma_p(X,2)
    return(norma)

def norma_filasum(X):
    sfila = np.sum(np.abs(X),axis=1)
    norma = np.max(sfila)
    return(norma)

def norma_frobenius(X):
    tamano = np.shape(X)
    n = tamano[0]
    m = tamano[1]
    norma = 0
    for i in range(0,n,1):
        for j in range(0,m,1):
            norma =  norma + np.abs(X[i,j])**2
    norma = np.sqrt(norma)
    return(norma)

def num_condicion(X):
    M = np.copy(X)
    Mi = np.linalg.inv(M)
    nM = norma_filasum(M)
    nMi= norma_filasum(Mi)
    ncondicion = nM*nMi
    return(ncondicion)

# Programa de prueba #######
# INGRESO
A = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])

B = np.array([7.85,-19.3,71.4])

p = 2

# PROCEDIMIENTO
normap = norma_p(B, p)
normaeucl = norma_euclidiana(B)
normafilasuma = norma_filasum(A)
numerocondicion = num_condicion(A)

# SALIDA
print('vector:',B)
print('norma p: ',2)
print(normap)

print('norma euclididana: ')
print(normaeucl)

print('******')
print('matriz: ')
print(A)
print('norma suma fila: ',normafilasuma)

print('número de condición:')
print(numerocondicion)

cuyos resultados del ejercicio serán:

vector: [  7.85 -19.3   71.4 ]
norma p:  2
74.3779033047
norma euclididana: 
74.3779033047
******
matriz: 
[[  3.   -0.1  -0.2]
 [  0.1   7.   -0.3]
 [  0.3  -0.2  10. ]]
norma suma fila:  10.5
número de condición:
3.61442432483
>>> 

compare sus resultados con las funciones numpy:

np.linalg.norm(A)
np.linalg.cond(A)

http://www.numpy.org/devdocs/reference/generated/numpy.linalg.norm.html

http://www.numpy.org/devdocs/reference/generated/numpy.linalg.cond.html

4.2 Gauss-Seidel – Algoritmo

El método se describe con un sistema de 3 incógnitas y 3 ecuaciones.
Note el uso de índices ajustado a la descripción de python para la primera fila=0 y primera columna=0:

a00x0 + a01x1+ a02x2 = b0
a10x0 + a11x1+ a12x2 = b1
a20x0 + a21x1+ a22x2 = b2

La forma matricial del sistema original de ecuaciones es: [A]{X} = {B}

Se usa los indices i para fila y j para columna, asi un elemento de la matriz es aij .
Para obtener los valores de cada xi, se despeja uno por cada por cada fila i:

x0 = (b0        - a01x1 - a02x2) /a00
x1 = (b1 - a10x0        - a12x2) /a11
x2 = (b2 - a20x0 - a21x1         ) /a22

Observe el patrón de cada fórmula usada para determinar x[i]:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}
la sumatoria 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 = [0, 0, 0]

Con cada nuevo valor se calcula el vector diferencias entre el vector X y cada nuevo valor calculado X[i] .

El error que llama la atención es al mayor valor de las diferencias, y 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.

Con lo descrito, se muestra una forma de implementar el algoritmo como una función y un ejemplo de prueba con datos del problema conocido.

# Algoritmo Gauss-Seidel,
# matrices, métodos iterativos
# Referencia: Chapra 11.2, p.310, pdf.334
#      Rodriguez 5.2 p.162
# ingresar iteramax si requiere más iteraciones

import numpy as np

def gauss_seidel(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    
    itera = 0
    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]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=0
    return(X)

# Programa de prueba #######
# INGRESO
A = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])

B = np.array([7.85,-19.3,71.4])

tolera = 0.00001

# PROCEDIMIENTO
n = len(B)
X = np.zeros(n, dtype=float)
respuesta = gauss_seidel(A,B,tolera,X)
verifica = np.dot(A,respuesta)

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

que da como resultado:

respuesta de A.X=B : 
[[ 3. ]
 [-2.5]
 [ 7. ]]
verificar A.X: 
[[  7.84999999]
 [-19.3       ]
 [ 71.4       ]]
>>> 

que es la respuesta del problema obtenida durante el desarrollo del ejemplo con otros métodos.