Categoría: Soluciones

Ejercicio: 1Eva_IIT2007_T2 Aplicar Gauss-Seidel 6×6

1. Desarrollo Analítico

– Verificar que la matriz es diagonal dominante

No es necesario realizar el pivoteo por filas, ya la matriz tiene la diagonal dominante.

A = [[7.63, 0.30, 0.15,  0.50, 0.34, 0.84],
     [0.38, 6.40, 0.70,  0.90, 0.29, 0.57],
     [0.83, 0.19, 8.33,  0.82, 0.34, 0.37],
     [0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
     [0.71, 0.30, 0.85,  0.82, 5.95, 0.55],
     [0.43, 0.54, 0.59,  0.66, 0.31, 9.25]]

B = [ -9.44, 25.27, -48.01, 19.76, -23.63, 62.59]

– revisar el número de condición

cond(A) = || A||.||A^-1||

El número de condición no es «muy alto»,  los valores de la diagonal son los mayores en toda la fila, por lo que el sistema converge.

>>> np.linalg.cond(A)
2.04518539662910119

– realizar iteraciones con las expresiones:

x_0 = \Big(-9.44 -0.3 x_1 -0.15x_2 -0.50 x_3 -0.34 x_4 -0.84 x_5\Big)\frac{1}{7.63} x_1 = \Big( 25.27 -0.38 x_0 -0.70 x_2 -0.90 x_3 -0.29 x_4 -0.57 x_5 \Big) \frac{1}{6.40} x_2 = \Big(-48.01 -0.83x_0 -0.19 x_1 -0.82x_3 -0.34x_4 -0.37x_5 \Big)\frac{1}{8.33} x_3 = \Big(19.76 -0.50 x_0 -0.68 x_1 -0.86 x_2 -0.53 x_4 -0.70x_5 \Big)\frac{1}{10.21} x_4 = \Big( -23.63 - 0.71 x_0 -0.30 x_1 -0.85 x_2 -0.82 x_3 -0.55x_5 \Big)\frac{1}{5.95} x_5 = \Big( 62.59 - 0.43x_0 -0.54x_1 - 0.59 x_2 -0.66 x_3 -0.31 x_4 \Big)\frac{1}{9.25}

Dado que no se establece en el enunciado el vector inicial, se usará el vector cero. La tolerancia requerida es 10^-5
X0 = [ 0. 0. 0. 0. 0. 0.] = [x₀,x₁,x₂,x₃,x₄,x₅]

iteración 0

x_0 = \Big(-9.44 -0.3 (0) -0.15(0) -0.50 (0) -0.34(0) -0.84 (0)\Big)\frac{1}{7.63} = -1.23722 x_1 = \Big( 25.27 -0.38 (-1.23722) -0.70 (0) -0.90 (0) -0.29 (0) -0.57 (0) \Big) \frac{1}{6.40} = 4.0219 x_2 = \Big(-48.01 -0.83 (-1.23722) -0.19 (4.0219) -0.82(0) -0.34(0) -0.37(0) \Big)\frac{1}{8.33} = -5.73196 x_3 = \Big( 19.76 -0.50 (-1.23722) -0.68 (4.0219) -0.86 (-5.73196) -0.53 (0) -0.70(0) \Big)\frac{1}{10.21} = 2.21089 x_4 = \Big( -23.63 - 0.71 (-1.23722) -0.30 (4.0219) -0.85 (-5.73196) -0.82 (2.21089) -0.55(0) \Big)\frac{1}{5.95} = -3.51242 x_5 = \Big( 62.59 - 0.43(-1.23722) -0.54(4.0219) - 0.59 (-5.73196) -0.66 (2.21089) -0.31 (-3.51242) \Big)\frac{1}{9.25} = 6.91478 X_1 = [-1.23722, 4.0219, -5.73196, 2.21089, -3.51242, 6.91478 ] diferencia = X1 - [0,0,0,0,0,0] errado = max(|diferencia|) = 6.91478

iteración 1

x_0 = \Big(-9.44 -0.3 (4.0219) -0.15(4.0219) -0.50 (2.21089) -0.34 (-3.51242) -0.84(6.91478) \Big)\frac{1}{7.63} = -2.0323 x_1 = \Big( 25.27 -0.38 (-2.0323) -0.70 ( -5.73196) -0.90 (2.21089) -0.29 (-3.51242) -0.57 (6.91478) \Big) \frac{1}{6.40} = 3.92844 x_2 = \Big(-48.01 -0.83(-2.0323) -0.19 (3.92844) -0.82(2.21089) -0.34(-3.51242) -0.37(6.91478) \Big)\frac{1}{8.33} = -6.03203 x_3 = \Big(19.76 -0.50 (-2.0323) -0.68 (3.92844) -0.86 (-6.03203) -0.53 (-3.51242) -0.70(6.91478) \Big)\frac{1}{10.21} = 1.98958 x_4 = \Big( -23.63 - 0.71 (-2.0323) -0.30 (3.92844) -0.85 (-6.03203) -0.82 (1.98958) -0.55(6.91478) \Big)\frac{1}{5.95} = -3.97865 x_5 = \Big( 62.59 - 0.43(-2.0323) -0.54(3.92844) - 0.59 (-6.03203) -0.66 (1.98958) -0.31 (-3.97865) \Big)\frac{1}{9.25} = 7.00775 X_2 = [-2.0323, 3.92844, -6.03203, 1.98958, -3.97865, 7.00775] diferencia = X2 - [-1.23722, 4.0219, -5.73196, 2.21089, -3.51242, 6.91478 ] errado = max(|diferencia|) = 0.79507

iteración 2

Desarrollar como tarea.

---

2. Algoritmo en Python

La tabla de aproximaciones sucesivas para el vector X es:

Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X],[errado]
0 [0. 0. 0. 0. 0. 0.] [1.]
1 [-1.23722  4.0219  -5.73196  2.21089 -3.51242  6.91478] [6.91478]
2 [-2.0323   3.92844 -6.03203  1.98958 -3.97865  7.00775] [0.79507]
3 [-1.99768  4.00317 -6.00049  1.99808 -4.00082  6.9999 ] [0.07473]
4 [-1.99994  4.00037 -5.99979  2.      -4.00005  6.99996] [0.00281]
5 [-2.00001  3.99998 -6.       2.00001 -4.       7.     ] [0.00038]
6 [-2.  4. -6.  2. -4.  7.] [1.60155e-05]
7 [-2.  4. -6.  2. -4.  7.] [1.74554e-06]
Matriz aumentada y pivoteada:
[[  7.63   0.3    0.15   0.5    0.34   0.84  -9.44]
 [  0.38   6.4    0.7    0.9    0.29   0.57  25.27]
 [  0.83   0.19   8.33   0.82   0.34   0.37 -48.01]
 [  0.5    0.68   0.86  10.21   0.53   0.7   19.76]
 [  0.71   0.3    0.85   0.82   5.95   0.55 -23.63]
 [  0.43   0.54   0.59   0.66   0.31   9.25  62.59]]
Vector Xi: 
[-2.  4. -6.  2. -4.  7.]
>>> 

el gráfico de los iteraciones vs errores es:

La tabla obtiene aplicando la función de Gauss-Seidel, tomando como vector inicial el vector de ceros.

Tarea: X=TX+C

Instrucciones en Python

# 1Eva_IIT2007_T2 Aplicar Gauss-Seidel
# Algoritmo Gauss-Seidel
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    M = np.copy(AB)
    
    # Pivoteo por filas AB
    tamano = np.shape(M)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(M[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(M[i,:])
            M[i,:] = M[dondemax+i,:]
            M[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(pivoteado, 'intercambiar: ',i,dondemax)
    if vertabla==True and pivoteado==0:
        print('Pivoteo por filas NO requerido')
    return(M)

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

    itera = 0
    if vertabla==True:
        print('Iteraciones Gauss-Seidel')
        print('itera,[X],[errado]')
        np.set_printoptions(precision)
        print(itera, X, np.array([errado]))
    while (errado>tolera and itera<iteramax):
        for i in range(0,n,1):
            xi = B[i]
            for j in range(0,m,1):
                if (i!=j):
                    xi = xi-A[i,j]*X[j]
            xi = xi/A[i,i]
            diferencia[i] = np.abs(xi-X[i])
            X[i] = xi
        errado = np.max(diferencia)
        itera = itera + 1
        if vertabla==True:
            print(itera, X, np.array([errado]))        
    
    if (itera>iteramax): # No converge
        X = itera
        print('iteramax superado, No converge')
    return(X)

# Programa de prueba #######
# INGRESO
A = np.array([[7.63, 0.30, 0.15,  0.50, 0.34, 0.84],
              [0.38, 6.40, 0.70,  0.90, 0.29, 0.57],
              [0.83, 0.19, 8.33,  0.82, 0.34, 0.37],
              [0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
              [0.71, 0.30, 0.85,  0.82, 5.95, 0.55],
              [0.43, 0.54, 0.59,  0.66, 0.31, 9.25]])

B = np.array([-9.44,25.27,-48.01,19.76,-23.63,62.59])

tolera = 0.00001
X = np.zeros(len(A), dtype=float)

# PROCEDIMIENTO
n = len(A)
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)


AB = pivoteafila(A,B, vertabla=True)
A = AB[:,:n]
B = AB[:,n]
respuesta = gauss_seidel(A,B, tolera, X, vertabla=True)

# SALIDA
print('Matriz aumentada y pivoteada:')
print(AB)
print('Vector Xi: ')
print(respuesta)

En el caso de la norma infinito, para la matriz A, se puede usar el algoritmo desarrollado en clase.
Como valor para verificar su algoritmo, se obtuvo:

>>> np.linalg.norm(A, np.inf)
13.479999999999999

Tarea: incluir la norma infinito para T

Ejercicio: 1Eva_IIT2007_T1 Distribución binomial acumulada

Dado F=0.4, dado que n=5 y k=1

F = \sum_{t=0}^{k} \binom{n}{t} p^t (1-p)^{n-t}

La fórmula para el ejercicio se convierte en:

F = \Bigg( \begin{array}{c} 5 \\ 0 \end{array} \Bigg) p ^0 (1-p)^{(5-0)} + \Bigg( \begin{array}{c} 5 \\ 1 \end{array} \Bigg) p ^1 (1-p)^{(5-1)} = 0.4

Los valores de las combinatorias se calculan como:

>>> import scipy.special as sts
>>> sts.comb(5,0,repetition=False)
1.0
>>> sts.comb(5,1,repetition=False)
5.0
>>> 

(1-p)^{5} + 5p (1-p)^{4} = 0.4

La expresión para el ejercicio se convierte en:

f(p) = (1-p)^{5} + 5p (1-p)^{4} - 0.4

como referencia se revisa la gráfica para f(p)

f(p) = (1-p)^5 + 5p(1-p)^4 - 0.4 = (1-p)^4 (1 - p + 5p) - 0.4 = (1-p)^4 (1 + 4p) - 0.4 = (1-p)^2 (1-p)^2 (1 + 4p) - 0.4 = (1-2p+p^2) (1-2p+p^2) (1 + 4p) - 0.4 = (1 - 4p + 6p^2 - 4p^3 +p^4 ) (1 + 4p) - 0.4 = 1 - 10p^2 + 20p^3 + 15p^4 + 4p^5 - 0.4 f(p) = 0.6 - 10p^2 + 20p^3 + 15p^4 + 4p^5

y su derivada:

f'(p) = - 20p + 60p^2 + 60p^3 +20p^4

con lo que se puede desarrollar el método Newton-Raphson.

Verificando el polinomio  obtenido a partir de la expresión inicial usando Sympy:

>>> import sympy as sp
>>> p = sp.Symbol('p')
>>> poli = (1-p)**5 + 5*p*((1-p)**4) - 0.4
>>> pol = poli.expand()
>>> pol
4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
>>> pol.diff(p,1)
20*p**4 - 60*p**3 + 60*p**2 - 20*p

A partir de la gráfica, un punto inicial cercano a la raíz es X0 = 0.2

itera = 0

f(0.2) = 0.6 - 10(0.2)^2 + 20(0.2)^3 + 15(0.2)^4 + 4(0.2)^5 = 0.3373 f'(0.2) = - 20(0.2) + 60(0.2)^2 + 60(0.2)^3 +20(0.2)^4 = -2.048 x_1 = 0.2 -\frac{0.3373}{-2.048} = 0.3647 errado = |0.2 - 0.36469| = 0.1647

itera = 1

f(0.36469) = 0.6 - 10(0.36469)^2 + 20(0.36469)^3 + 15(0.36469)^4 + 4(0.36469)^5 = 0.0005566 f'(0.36469) = - 20(0.36469) + 60(0.36469)^2 + 60(0.36469)^3 +20(0.36469)^4 = -1.8703 x_1 = 0.36469 -\frac{0.0005566}{-1.8703} = 0.36499 errado = |0.36469 - 0.36499| = 0.000297

itera = 3

f(0.36499) = 0.6 - 10(0.36499)^2 + 20(0.36499)^3 + 15(0.36499)^4 + 4(0.36499)^5 = 1.6412291237166698e-07 f'(0.36499) = - 20(0.36499) + 60(0.36499)^2 + 60(0.36499)^3 +20(0.36499)^4 = -1.869204616112814 x_1 = 0.36469 -\frac{1.6412291237166698e-07}{ -1.8692} = 0.36499 errado = |0.36499 - 0.36499| = 8.7804e-08

verificar con raíz: 0.3649852264049102

Instrucciones para la gráfica

# 1ra Eval II Término 2007
# Tema 1. Distribución binomial acumulada
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sts

fp = lambda p: (1-p)**5 + 5*p*((1-p)**4) - 0.4

a = 0
b = 1
pasos = 100

# PROCEDIMIENTO
xi = np.linspace(a,b,pasos+1)
p_i = fp(xi)

# SALIDA
plt.plot(xi,p_i)
plt.axhline(0)
plt.show()

---

Algoritmo en Python

el resultado usando el algoritmo es:

i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 0.2     0.3373 -2.048   0.3647  0.1647]
1 [ 3.6469e-01  5.5668e-04 -1.8703e+00  3.6499e-01  2.9764e-04]
2 [ 3.6499e-01  1.6412e-07 -1.8692e+00  3.6499e-01  8.7804e-08]
raiz en:  0.3649852264049102

con error de: 8.780360960525257e-08

Instrucciones en Python para Newton-Raphson

# 1Eva_IIT2007_T1 Distribución binomial acumulada
# Método de Newton-Raphson

import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx,dfx,xi, tolera, iteramax=100, vertabla=False, precision=4):
    '''
    fx y dfx en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('método de Newton-Raphson')
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1

    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera, 'No converge,se alcanzó el máximo de iteraciones')

    return(xi)

# INGRESO
fx  = lambda p: 4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
dfx = lambda p: 20*p**4 - 60*p**3 + 60*p**2 - 20*p

x0 = 0.2
tolera = 0.0000001

# PROCEDIMIENTO
respuesta = newton_raphson(fx,dfx,x0, tolera, vertabla=True)
# SALIDA
print('raiz en: ', respuesta)

# GRAFICA
a = 0
b = 1
muestras = 21

xi = np.linspace(a,b,muestras)
fi = fx(xi)
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()