3.3.3 Método con Matrices triangulares A=L.U

Referencia: Chapra 10.2 p287, pdf312

Se plantea resolver el sistema de ecuaciones usando matrices triangulares

A = \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix} B = [7.85,-19.3,71.4]

Para encontrar la solución al sistema de ecuaciones, se plantea resolver:
– sustitución hacia adelante de LY=B
– sustitución hacia atras para UX=Y

Instrucciones en Python

Al algoritmo de la sección anterior se añade los procedimientos correspondientes con los que se obtiene la solución:

[[ 3. ]
 [-2.5]
 [ 7. ]]
# Solución usando Matrices L y U
# continuación de ejercicio anterior

import numpy as np

# INGRESO
A = np.array([[ 3. , -0.1, -0.2],
              [ 0.1,  7. , -0.3],
              [ 0.3, -0.2, 10. ]], dtype=float)

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

# PROCEDIMIENTO
B  = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
AB = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
for i in range(0,n-1,1):

    # columna desde diagonal i en adelante
    columna = abs(AB[i:,i])
    dondemax = np.argmax(columna)

    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:] = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

AB1 = np.copy(AB)
A1 = np.copy(AB[:,:m-1])
B1 = np.copy(AB[:,m-1])

# eliminacion hacia adelante
# se inicializa L
L = np.identity(n,dtype=float)
for i in range(0,n-1,1):
    pivote = AB[i,i]
    adelante = i+1
    for k in range(adelante,n,1):
        factor = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor
        L[k,i] = factor

U = np.copy(AB[:,:m-1])

# Resolver LY = B
B1  = np.transpose([B1])
AB =np.concatenate((L,B1),axis=1)

# sustitución hacia adelante
Y = np.zeros(n,dtype=float)
Y[0] = AB[0,n]
for i in range(1,n,1):
    suma = 0
    for j in range(0,i,1):
        suma = suma + AB[i,j]*Y[j]
    b = AB[i,n]
    Y[i] = (b-suma)/AB[i,i]

Y = np.transpose([Y])

# Resolver UX = Y
AB =np.concatenate((U,Y),axis=1)

# sustitución hacia atrás
ultfila = n-1
ultcolumna = m-1
X = np.zeros(n,dtype=float)

for i in range(ultfila,0-1,-1):
    suma = 0
    for j in range(i+1,ultcolumna,1):
        suma = suma + AB[i,j]*X[j]
    b = AB[i,ultcolumna]
    X[i] = (b-suma)/AB[i,i]

X = np.transpose([X])

# SALIDA
print('Pivoteo parcial por filas: AB')
print(AB1)
print('eliminación hacia adelante')
print('Matriz U: ')
print(U)
print('matriz L: ')
print(L)
print('B1 :')
print(B1)
print('Y Sustitución hacia adelante')
print(Y)
print('X Sustitución hacia atras')
print(X)

2.6 Sistemas de Ecuaciones no lineales – Newton-Raphson

Referencia: Chapra 6.5 p162, Chapra Ejercicio 6.11 p.166/pdf190

Con el método de Newton-Raphson para múltiples ecuaciones, determine las raíces para:

x^2+xy =10 y + 3xy^2 = 57

Observe que un par correcto de raíces es x=2 y y=3.
Use como valore iniciales x=1.5, y=3.5

Planteamiento

Las ecuaciones se expresan de la forma f(x,y) = 0

x^2+xy -10 = 0 y + 3xy^2 -57 = 0

Se puede usar extensines de los métodos abiertos para resolver ecuacioens simples, por ejemplo Newton-Raphson.

u_{i+1} = u_i + (x_{i+1}-x_i)\frac{\partial u_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial u_i}{\partial y} v_{i+1} = v_i + (x_{i+1}-x_i)\frac{\partial v_i}{\partial x} + (y_{i+1}-y_i) \frac{\partial v_i}{\partial y}

ecuaciones que se pueden reordenar y encontrar la solución a partir de la matriz Jacobiano.

Instrucciones en Python

Usando un algoritmo para resolver el Jacobiano y estimar los puntos luego de cada iteración se obtienen:

iteración:  1
Jacobiano con puntos iniciales: 
Matrix([[6.50000000000000, 1.50000000000000], [36.7500000000000, 32.5000000000000]])
determinante:  156.12499999999994
puntos xi,yi: 2.03602882305845 2.84387510008006
error: 0.656124899919936
iteración:  2
Jacobiano con puntos iniciales: 
Matrix([[6.91593274619696, 2.03602882305845], [24.2628767545662, 35.7412700376474]])
determinante:  197.78430344142245
puntos xi,yi: 1.99870060905582 3.00228856292451
error: 0.158413462844444
iteración:  3
Jacobiano con puntos iniciales: 
Matrix([[6.99968978103616, 1.99870060905582], [27.0412098452019, 37.0040558756713]])
determinante:  204.96962918261596
puntos xi,yi: 1.99999998387626 2.99999941338891
error: 0.00228914953559523
iteración:  4
Jacobiano con puntos iniciales: 
Matrix([[6.99999938114143, 1.99999998387626], [26.9999894410015, 36.9999926704397]])
determinante:  204.9999473486533
puntos xi,yi: 1.99999999999998 3.00000000000008
error: 5.86611161867978e-7
Resultado: 
1.99999999999998 3.00000000000008
>>> 

Algoritmo presentado para dos ecuaciones y dos incógnitas, en la unidad 3 se puede ampliar la propuesta. Revisar el método de Gauss-Seidel y Jacobi.

# Ejercicio Chapra Ej:6.11
# Sistemas de ecuaciones no lineales
# con método de Newton Raphson para xy

import numpy as np
import sympy as sym

def matrizJacobiano(variables, funciones):
    n = len(funciones)
    m = len(variables)
    # matriz Jacobiano inicia con ceros
    Jcb = sym.zeros(n,m)
    for i in range(0,n,1):
        unafi = sym.sympify(funciones[i])
        for j in range(0,m,1):
            unavariable = variables[j]
            Jcb[i,j] = sym.diff(unafi, unavariable)
    return Jcb

# PROGRAMA ----------
# INGRESO
x = sym.Symbol('x')
y = sym.Symbol('y')

f1 = x**2 + x*y - 10
f2 = y + 3*x*(y**2)-57

x0 = 1.5
y0 = 3.5

tolera = 0.0001

# PROCEDIMIENTO
funciones = [f1,f2]
variables = [x,y]
n = len(funciones)
m = len(variables)

Jxy = matrizJacobiano(variables, funciones)

# valores iniciales
xi = x0
yi = y0

# tramo inicial, mayor que tolerancia
itera = 0
tramo = tolera*2

while (tramo>tolera):
    J = Jxy.subs([(x,xi),(y,yi)])

    # determinante de J
    Jn = np.array(J,dtype=float)
    determinante =  np.linalg.det(Jn)

    # iteraciones
    f1i = f1.subs([(x,xi),(y,yi)])
    f2i = f2.subs([(x,xi),(y,yi)])

    numerador1 = f1i*Jn[n-1,m-1]-f2i*Jn[0,m-1]
    xi1 = xi - numerador1/determinante
    numerador2 = f2i*Jn[0,0]-f1i*Jn[n-1,0]
    yi1 = yi -numerador2/determinante
    
    tramo = np.max(np.abs([xi1-xi,yi1-yi]))
    xi = xi1
    yi = yi1

    itera = itera +1
    print('iteración: ',itera)
    print('Jacobiano con puntos iniciales: ')
    print(J)
    print('determinante: ', determinante)
    print('puntos xi,yi:',xi,yi)
    print('error:',tramo)
    
# SALIDA
print('Resultado: ')
print(xi,yi)

5.4.1 Cuadratura con dos puntos – Experimento

Para visualizar el concepto de Cuadratura de Gauss de dos puntos considere lo siguiente:

Se tiene un corte transversal del un recipiente rectangular lleno de líquido limitado en x entre [-1,1], al que se le ha depositado encima otro recipiente con perfil f(x) = x2 hasta que reposan sus extremos en x[-1,1].

La altura de ambos recipientes es la misma.

La superficie entre f(x) y el eje x es el integral de f(x) en el intervalo.

Si suponemos que la figura es el corte transversal de una vasija y la parte en amarillo es líquida, la vasija ha desplazado el líquido que ocupa ahora el «área» mostrada en la gráfica que corresponde al integral de f(x)=x2. entre [-1,1].

Ahora, suponga que se perfora el perfil de f(x) en dos puntos equidistantes cercanos  x=0. Los orificios permitirían el desplazamiento del liquido al interior de f(x) que dejando pasar suficiente tiempo, perimitiría tener todo el líquido en el recipiente rectangular entre [-1,1] como una línea horizontal.

Podría medir la altura que tiene el líquido y que tiene un equivalente en un punto f(x1). Debería encontrar el valor de x1 que permite disponer del mismo valor entre el área debajo de f(x) y el rectángulo del corte transversal amarillo ahora formado.

Se usa el resultado analítico del integral restando el área del rectángulo obtenido al evaluar la funcion f(x) entre [0,1], teniendo un problema de busqueda de raíces. Obtenemos el valor de x1.

Se muestra que el área bajo f(x) es equivalente al área del rectángulo conformado.

Si  utilizamos el desplazamiento horizontal desde el centro para un punto encontrado como un «factor», tendremos que el área del rectángulo se mantendría equivalente, y el desplazamiento proporcional a la mitad del intervalo si se varía el intervalo de observacion. Este factor coincide con el factor de Cuadratura de Gauss de dos puntos.

funcion  fx:   x**2
Integral Fx:   x**3/3
I analitico:   0.6666666666666666
I aproximado:  0.6666666666654123
desplaza centro:   0.5773502691890826
factor desplaza:   0.5773502691890826
Factor CuadGauss:  0.5773502691896258
erradoFactor:    1.2545520178264269e-12
error integral:  1.2543299732215019e-12

El error del integral es del orden de 10-12


Cambiamos la figura geométrica a un trapecio generado por la recta que pasa por los puntos xi desplazados desde el centro.

Usamos la función f(x) = x2 + x + 1, observaremos si los resultado son equivalentes.

La figura al inicio del experimento será:

Luego de realizar realizar el mismo cálculo anterior usando un equivalente a trapecio se tiene:

con valores numéricos:

funcion  fx:   x**2 + x + 1
Integral Fx:   x**3/3 + x**2/2 + x
I analitico:   2.6666666666666665
I aproximado:  2.6666666666654124
desplaza centro:   0.5773502691890826
factor desplaza:   0.5773502691890826
Factor CuadGauss:  0.5773502691896258
erradoFactor:    1.2545520178264269e-12
error integral:  1.2541079286165768e-12

El error del integral es también del orden de 10-12, además observe que el factor de cuadratura de Gauss se mantiene.


Tarea

Realice el experimento usando un polinomio de grado superior y observe los errores para el integral y las diferencia con el coeficiente de Cuatratura de 2 puntos.


Algoritmo con Python

Para resumir la cantidad de instrucciones, se usa el método de la bisección desde la librería scipy y el subgrupo de funciones de optimización.

Los cálculos para realizar las gráficas se tratan en un bloque luego de mostrar los resultado principales.

# Integración: Cuadratura de Gauss de dos puntos
# modelo con varios tramos entre [a,b]
# para un solo segmento.
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import scipy.optimize as op

# INGRESO
x = sym.Symbol('x')

# fx = (x)**2
fx = x**2 + x + 1
# fx = 0.2 + 25.0*x-200*(x**2)+675.0*(x**3)-900.0*(x**4)+400.0*(x**5)

a = -1
b = 1
muestras = 51
tolera = 1e-12
iteramax = 100

# PROCEDIMIENTO
# Desarrollo analítico con Sympy
Fx = sym.integrate(fx,x)
Fxn = sym.lambdify('x',Fx,'numpy')
Fiab = Fxn(b)-Fxn(a)

# Busca igualar trapecio con Integral analitico
fxn = sym.lambdify('x',fx,'numpy')
base = b-a
mitad = base/2
xc = (a+b)/2  # centro
diferencia = lambda x: Fiab-base*(fxn(xc-x)+fxn(xc+x))/2
desplazado =  op.bisect(diferencia,0,mitad,
                        xtol=tolera,maxiter=iteramax)
factor = desplazado/mitad

# Integral aproximando con trapecio
x0 = xc - factor*mitad
x1 = xc + factor*mitad
Faprox = base*(fxn(x0)+fxn(x1))/2

# Integral cuadratura Gauss
xa = xc + mitad/np.sqrt(3)
xb = xc - mitad/np.sqrt(3)
FcuadG = base*(fxn(xa)+fxn(xb))/2
erradofactor = np.abs(FcuadG - Faprox)
erradoIntegral = np.abs(Fiab-Faprox)
# SALIDA
print('funcion  fx:  ', fx)
print('Integral Fx:  ', Fx)
print('I analitico:  ', Fiab)
print('I aproximado: ', Faprox)
print('desplaza centro:  ', desplazado)
print('factor desplaza:  ', factor)
print('Factor CuadGauss: ', 1/np.sqrt(3))
print('erradoFactor:   ', erradofactor)
print('error integral: ', erradoIntegral) 

# Grafica
# Para GRAFICAR 
# Para gráfica f(x)
xi = np.linspace(a,b,muestras)
fi = fxn(xi)

# Para gráfica Trapecio
m = (fxn(x1)-fxn(x0))/(x1-x0)
trapeciof = lambda x: fxn(x0)+m*(x-x0)
trapecioi = trapeciof(xi)

# Areas Trapecio para cada punto que busca
k = int(muestras/2)
xicg = xi[k:muestras-1]
Fcg = [base*(fxn(xi[k+0])+fxn(xi[k-0]))/2]
for i in range(1,k,1):
    untrapecio = base*(fxn(xi[k+i])+fxn(xi[k-i]))/2
    Fcg.append(untrapecio)

# Punto buscado
Fiaprox = base*(fxn(x1)+fxn(x0))/2

Fi = Fxn(xi)-Fxn(a)

# Areas de curvas y trapecio

plt.subplot(211) # Grafica superior
plt.xlim(a,b)
plt.plot(xi, fi, label='f(x)')
# Solo fi
# plt.fill_between(xi,0, fi,
#                label='integral fi',
#                 color='yellow')
# usando cuadratura
plt.fill_between(xi,0, trapecioi,
                 label='Cuadratura 2 puntos',
                 color='yellow')
plt.axvline(x0,color='white')
plt.axvline(x1,color='white')
plt.plot([x0,x1],[fxn(x0),fxn(x1)],
         'ro', label='x0,x1')
plt.axvline(0,color='black')
plt.xlabel('x')
plt.ylabel('f(x) y Cuadratura de 2 puntos')
plt.legend()

# Valores de integrales
plt.subplot(212) # Grafica inferior
plt.xlim(a,b)
plt.axhline(Fiab, label='F[a,b]')
# plt.plot(xi,Fi,label='F(x)')
plt.plot(xicg,Fcg,color='orange',label='Aprox Trapecio')

plt.axvline(x1,color='yellow')
plt.axvline((1/np.sqrt(3))*(b-a)/2 + xc ,color='magenta')
plt.plot(x1,Fiaprox,'ro', label='x0,x1')
plt.axvline(0,color='black')

plt.xlabel('x')
plt.legend()
plt.ylabel('Integrando')
plt.show()

3.6.1 Método de Gauss-Seidel. Ejemplo01

Referencia: Ejemplo 01 Chapra Ejemplo 11.3  p311 pdf335

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fué implementado el procedimiento.


Planteamiento

Para el desarrollo con lápiz y papel, se despeja una variable de cada ecuación, se empieza con la primera expresión para obtener x0

3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 3 x_0 = 7.85 + 0.1 x_1 + 0.2 x_2 x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3}

con la segunda ecuación se obtiene x1

0.1 x_0 +7x_1 -0.3 x_2 = -19.3 7x_1 = -19.3 - 0.1 x_0 +0.3 x_2 x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7}

usando la tercera ecuación se obtiene x2

0.3 x_0 - 0.2 x_1 + 10 x_2 = 71.4 10 x_2 = 71.4 - 0.3 x_0 + 0.2 x_1 x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}

Vector inicial

Como el vector inicial no se indica en el enunciado, se considera usar el vector de ceros para iniciar  las iteraciones.

X = [0,0,0]

con tolerancia de 0.00001


Iteraciones

Con las ecuaciones obtenidas en el planteamiento, se desarrolla usando el vector inicial presentado, hasta que el |error|<tolera

itera = 0

x_0 = \frac{7.85 + 0.1 (0) + 0.2 (0)}{3} = 2.61 x_1 = \frac{ -19.3 - 0.1 (2.61) +0.3 (0)}{7} = -2.79 x_2 = \frac{71.4 - 0.3 (2.61) + 0.2 (-2.79)}{10} = 7.00

X = [2.61, -2.79, 7.00]

diferencias = [2.61, -2.79, 7.00] – [0,0,0]
diferencias = [2.61, -2.79, 7.00]
error = ||diferencias|| , se usa la norma de max(diferencias)
error = 7.00

itera = 1

X = [2.61, -2.79, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.79) + 0.2 (7.00)}{3} = 2.99 x_1 = \frac{ -19.3 - 0.1 (2.99) +0.3 (7.00)}{7} = -2.49 x_2 = \frac{71.4 - 0.3 (2.99) + 0.2 (-2.49)}{10} = 7.00

X = [2.99, -2.49, 7.00]

diferencias = [2.99, -2.49, 7.00] – [2.61, -2.79, 7.00]

diferencias = [0.38, 0.3 , 0. ]

error = ||diferencias||, se usa la norma de max(diferencias)

error = 0.38

itera = 2

X = [2.99, -2.49, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.49) + 0.2 (7.00)}{3} = 3.00 x_1 = \frac{ -19.3 - 0.1 (3.00) +0.3 (7.00)}{7} = -2.5 x_2 = \frac{71.4 - 0.3 (3.00) + 0.2 (-2.5)}{10} = 7.00

X = [3.00, -2.50, 7.00]

diferencias = [3.00, -2.50, 7.00] – [2.99, -2.49, 7.00]

diferencias = [ 0.01, -0.01, 0. ]

error = ||diferencias||, se usa la norma de max(diferencias)

error = 0.01

El error aún es mayor que tolera, por lo que es necesario continuar con las iteraciones.

Observaciones

El error disminuye en cada iteración, por lo que el método converge hacia la raiz.


Algoritmo con Python

Con la descripción dada para el método de Gauss-Seidel, se muestra una forma básica de implementar el algoritmo.

Ejemplo 01 Chapra Ejemplo 11.3 p311 pdf335

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fué implementado el procedimiento. Para generalizar el algoritmo, incluya como tarea aumentar el procedimiento de pivoteo por filas.

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

# 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])

X0  = np.array([0.,0.,0.])

tolera = 0.00001
iteramax = 100

# PROCEDIMIENTO

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        X[i] = nuevo
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])

# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

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

que da como resultado:

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

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


Tarea

Completar el algoritmo para usar una matriz diagonal dominante, usando al menos el pivoteo parcial por filas.

3.7 Método de Jacobi

Referencia: Burden 9Ed 7.3 p450, Rodriguez 5.1 p154

El método de Jacobi realiza operaciones semejantes al método de Gauss-Seidel.

El método de Jacobi también usa el vector inicial X0, la diferencia consiste en que la actualización del vector X en cada iteración se realiza cuando se ha calculado el vector nuevo completo.


Tarea

Realice el algoritmo de Jacobi en Python, las modificaciones al algoritmo de Gauss-Seidel, incluyendo el pivoteo por filas de la matriz.