Autor: Edison Del Rosario

1Eva_IIT2019_T1 Ecuación Recursiva

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

Tema 1. (30 puntos). Considere la sucesión
$\Big( x_n \Big)_{n=0}^{+ \infty}$
cuya ecuación recursiva es:
$x_n = g(x) = \sqrt{3 + x_{n-1}}$
para n ∈ Ν

a) ¿Se puede afirmar que ∀x ∈ [1,3], g(x) ∈ [1,3]?

b) Pruebe que g es una función contractiva en el intervalo [1,3] y estime el valor de la constante de Lipschitz (cota de la derivada de g)

c) Realice 5 iteraciones partiendo del dato inicial x₀ =2, y determine el orden de convergencia.

d) Encuentre el valor teórico de x^* al cual converge la sucesión y estime el error absolito en la iteración 5.

e) Realice 5 iteraciones con el método de bisección en el intervalo [1,3] para aproximar el punto fijo de la función g(x).

Rúbrica: literal a (3 puntos), literal b (3 puntos), literal c (10 puntos), literal d (4 puntos), literal e (10 puntos)

Referencia: Burden 9Ed. Definición 10.5 p633, Theorem 2.4 P62;
Contracción https://es.wikipedia.org/wiki/Contracci%C3%B3n_(espacio_m%C3%A9trico).
Función lipschitziana https://es.wikipedia.org/wiki/Funci%C3%B3n_lipschitziana

2019-11-27
1Eva_IIT2019_T2 Proceso Termodinámico

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

Tema 2. (20 puntos). Para simular la disminución de la temperatura en un proceso termodinámico,
un algoritmo evolutivo necesita usar un polinomio para aproximar en el intervalo [0,4] la función f con regla de correspondencia
$f(x)=e^{-0.5x}$
con constante k = 0.5

Para construir el mencionado polinomio, considere la tabla:

x 0 1 2 3 4

f(x) f(0) f(1) f(2) f(3) f(4)

a. Aplique interpolación polinomial y aproxime el valor de f(2.4) usando un polinomio de grado 2.

b. Encuentre una cota superior para el error de interpolación en la aproximación de f(1.7)}

Rúbrica: literal a (15 puntos), literal b (5 puntos)

2019-11-27
s1Eva_IIT2019_T1 Ecuación Recursiva
Ejercicio: 1Eva_IIT2019_T1 Ecuación Recursiva

la ecuación recursiva es:
x_n = g(x_{n-1}) = \sqrt{3 + x_{n-1}}

literal a y b

g(x) es creciente en todo el intervalo, con valor minimo en g(1) = 2, y máximo en g(3) =2.449. Por observación de la gráfica, la pendiente g(x) es menor que la recta identidad en todo el intervalo

Verifique la cota de g'(x)
g(x) = \sqrt{3 + x} =(3+x)^{1/2} g'(x) =\frac{1}{2}(3+x)^{-1/2} g'(x) =\frac{1}{2\sqrt{3+x}}
Tarea: verificar que g'(x) es menor que 1 en todo el intervalo.

Literal c

Usando el algoritmo del punto fijo, iniciando con el punto x₀=2
y tolerancia de 10^-4, se tiene que:

Iteración 1: x₀=2
g(x_0) = \sqrt{3 + 2} = 2.2361
error = |2.2361 - 2| = 0.2361

Iteración 2: x₁ = 2.2361
g(x_1) = \sqrt{3 + 2.2361} = 2.2882
error = |2.2882 - 2.2361| = 0.0522

Iteración 3: x₂ = 2.2882
g(x_2) = \sqrt{3 + 2.2882} = 2.2996
error = |2.2996 - 2.28821| = 0.0114

Iteración 4: x₃ = 2.2996
g(x_3) = \sqrt{3 + 2.2996} = 2.3021
error = |2.3021- 2.2996| = 0.0025

Iteración 5: x₄ = 2.3021
g(x_4) = \sqrt{3 + 2.3021} = 2.3026
error = |2.3021- 2.2996| = 5.3672e-04

con lo que determina que el error en la 5ta iteración es de 5.3672e-04 y el error se reduce en casi 1/4 entre iteraciones. El punto fijo converge a 2.3028

Se muestra como referencia la tabla resumen.
```
[[ x ,   g(x), 	 tramo  ] 
 [1.      2.      1.    ]
 [2.      2.2361  0.2361]
 [2.2361  2.2882  0.0522]
 [2.2882  2.2996  0.0114]
 [2.2996  2.3021  0.0025]
 [2.3021  2.3026  5.3672e-04]
 [2.3026  2.3027  1.1654e-04]
 [2.3027  2.3028  2.5305e-05]
raiz:  2.3027686193257098
```
con el siguiente comportamiento de la funcion:

literal e

Realizando el mismo ejercicio para el método de la bisección, se requiere cambiar a la forma f(x)=0
x = \sqrt{3 + x} 0 = \sqrt{3 + x} -x f(x) = \sqrt{3 + x} -x
tomando como intervalo el mismo que el inicio del problema [1,3], al realizar las operaciones se tiene que:
```
a = 1 ; f(a) = 1
b = 3 ; f(b) = -0.551
c = (a+b)/2 = (1+3)/2 = 2
f(c) = f(2) = (3 + 2)^(.5) +2 = 0.236
Siendo f(c) positivo, y tamaño de paso 2, se reduce a 1

a = 2 ; f(a) = 0.236
b = 3 ; f(b) = -0.551
c = (a+b)/2 = (2+3)/2 = 2.5
f(c) = f(2.5) = (3 + 2.5)^(.5) +2.5 = -0.155
Siendo fc(c) negativo y tamaño de paso 1, se reduce a .5

a = 2
b = 2.5
...
```
Siguiendo las operaciones se obtiene la siguiente tabla:
```
[ i, a,   c,   b,    f(a),  f(c),  f(b), paso]
 1 1.000 2.000 3.000 1.000  0.236 -0.551 2.000 
 2 2.000 2.500 3.000 0.236 -0.155 -0.551 1.000 
 3 2.000 2.250 2.500 0.236  0.041 -0.155 0.500 
 4 2.250 2.375 2.500 0.041 -0.057 -0.155 0.250 
 5 2.250 2.312 2.375 0.041 -0.008 -0.057 0.125 
 6 2.250 2.281 2.312 0.041  0.017 -0.008 0.062 
 7 2.281 2.297 2.312 0.017  0.005 -0.008 0.031 
 8 2.297 2.305 2.312 0.005 -0.001 -0.008 0.016 
 9 2.297 2.301 2.305 0.005  0.002 -0.001 0.008 
10 2.301 2.303 2.305 0.002  0.000 -0.001 0.004 
11 2.303 2.304 2.305 0.000 -0.001 -0.001 0.002 
12 2.303 2.303 2.304 0.000 -0.000 -0.001 0.001 
13 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
14 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
15 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
16 2.303 2.303 2.303 0.000  0.000 -0.000 0.000 
raiz:  2.302764892578125
```
Donde se observa que para la misma tolerancia de 10^-4, se incrementan las iteraciones a 16. Mientra que con punto fijo eran solo 8.

Nota: En la evaluación solo se requeria calcular hasta la 5ta iteración. Lo presentado es para fines didácticos
2019-11-27

s1Eva_IIT2019_T3 Circuito eléctrico

Ejercicio: 1Eva_IIT2019_T3 Circuito eléctrico

Las ecuaciones del problema son:

55 I_1 - 25 I_4 =-200

-37 I_3 - 4 I_4 =-250

-25 I_1 - 4 I_3 +29 I_4 =100

I_2 =-10

Planteo del problema en la forma A.X=B

A = [[ 55.0, 0,  0, -25],
     [  0  , 0,-37,  -4],
     [-25  , 0, -4,  29],
     [  0  ,  1, 0,   0]]

B = [-200,-250,100,-10]

El ejercicio se puede simplificar con una matriz de 3x3 dado que una de las corrientes I₂ es conocida con valor -10, queda resolver el problema para
[I₁ ,I₃ ,I₄ ]

A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

conformando la matriz aumentada

[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]

que se pivotea por filas para acercar a matriz diagonal dominante:

[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]

Literal a

Para métodos directos se aplica el método de eliminación hacia adelante.

Usando el primer elemento en la diagonal se convierten en ceros los números debajo de la posición primera de la diagonal

[[  55.    0.  -25.         -200.      ]
 [   0.  -37.   -4.         -250.      ]
 [   0.   -4.   17.636363      9.090909]]

luego se continúa con la segunda columna:

[[  55.    0.  -25.         -200.      ]
 [   0.  -37.   -4.         -250.      ]
 [   0.    0.   18.068796     36.117936]]

y para el método de Gauss se emplea sustitución hacia atrás
se determina el valor de I₄

18.068796 I_4 = 36.11793612

I_4 =\frac{36.11793612}{18.068796}= 1.99891216

-37 I_3 -4 I_4 = -250

-37 I_3= -250 + 4 I_4

I_3=\frac{-250 + 4 I_4}{-37}

I_3=\frac{-250 + 4 (1.99891216)}{-37} = 6.54065815

y planteando se obtiene el último valor

55 I_1 +25 I_4 = -200

55 I_1 = -200 -25 I_4

I_1 = \frac{-200 -25 I_4}{55}

I_1 = \frac{-200 -25(1.99891216)}{55} = -2.7277672

con lo que el vector solución es:

[-2.7277672   6.54065815  1.99891216]

sin embargo, para verificar la respuesta se aplica A.X=B

verificar que A.X = B, obteniendo nuevamente el vector B.
[-200.  -250.  100.]]

literal b

La norma de la matriz infinito se determina como:

||x|| = max\Big[ |x_i| \Big]

considere que en el problema el término en A de magnitud mayor es 55.
El vector suma de filas es:

[[| 55|+|  0|+|-25|],    [[80],
 [|  0|+|-37|+| -4|],  =  [41],
 [[-25|+| -4|+| 29|]]     [58]]

por lo que la norma ∞ ejemplo ||A||∞ 
es el maximo de suma de filas: 80

para revisar la estabilidad de la solución, se observa el número de condición

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

En éste caso no está muy alejado de 1. De resultar alejado del valor ideal de uno, la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.

Tarea: Matriz de transición de Jacobi

Literal c

En el método de Gauss-Seidel acorde a lo indicado, se inicia con el vector cero. Como no se indica el valor de tolerancia para el error, se considera tolera = 0.0001

las ecuaciones para el método son:

I_1 =\frac{-200 + 25 I_4}{55}

I_3 = \frac{-250+ 4 I_4}{-37}

I_4 =\frac{100 +25 I_1 + 4 I_3}{29}

Como I₂ es constante, no se usa en las iteraciones

I_2 =-10

teniendo como resultados de las iteraciones:

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
      diferencia,errado
0 [0. 0. 0.] 2e-05
1 [-3.6363636  6.7567568  1.2454461]
   [3.6363636 6.7567568 1.2454461] 6.756756756756757
2 [-3.0702518  6.6221139  1.7149021]
   [0.5661119 0.1346428 0.469456 ] 0.5661118513783094
3 [-2.8568627  6.5713619  1.891858 ]
   [0.2133891 0.050752  0.1769559] 0.2133891067340583
4 [-2.7764282  6.5522316  1.9585594]
   [0.0804345 0.0191304 0.0667014] 0.08043447732439457
5 [-2.7461094  6.5450206  1.9837016]
   [0.0303188 0.007211  0.0251423] 0.030318816370094925
6 [-2.7346811  6.5423025  1.9931787]
   [0.0114283 0.0027181 0.0094771] 0.011428316023939011
7 [-2.7303733  6.541278   1.996751 ]
   [0.0043078 0.0010246 0.0035723] 0.004307767346479974
8 [-2.7287495  6.5408918  1.9980975]
   [0.0016238 0.0003862 0.0013465] 0.001623761494915943
9 [-2.7281375  6.5407462  1.9986051]
   [0.0006121 0.0001456 0.0005076] 0.0006120575185017962
10 [-2.7279068  6.5406913  1.9987964]
   [2.3070778e-04 5.4871039e-05 1.9131760e-04] 0.00023070777766820427
11 [-2.7278198  6.5406707  1.9988685]
   [8.6962544e-05 2.0682983e-05 7.2114885e-05] 8.696254366213907e-05
12 [-2.727787   6.5406629  1.9988957]
   [3.2779493e-05 7.7962038e-06 2.7182845e-05] 3.277949307367578e-05
13 [-2.7277747  6.5406599  1.998906 ]
   [1.2355839e-05 2.9386860e-06 1.0246249e-05] 1.235583874370505e-05
14 [-2.72777    6.5406588  1.9989098]
   [4.6573860e-06 1.1077026e-06 3.8622013e-06] 4.6573859666665385e-06
numero de condición: 4.997509004325604
respuesta con Gauss-Seidel
[-2.72777    6.5406588  1.9989098]
>>>

con lo que el vector resultante es:

respuesta con Gauss-Seidel
[-2.72777 6.5406588 1.9989098]

que para verificar, se realiza la operación A.X
observando que el resultado es igual a B

[[-200.00002751]
 [-249.9999956 ]
 [ 100.0000125 ]]

Solución alterna

Usando la matriz de 4x4, los resultados son iguales para las corrientes
[I₁ ,I₂ , I₃ ,I₄ ]. Realizando la matriz aumentada,

[[  55.    0.    0.  -25. -200.]
 [   0.    0.  -37.   -4. -250.]
 [ -25.    0.   -4.   29.  100.]
 [   0.    1.    0.    0.  -10.]]

que se pivotea por filas para acercar a matriz diagonal dominante:

[[  55.    0.    0.  -25. -200.]
 [   0.    1.    0.    0.  -10.]
 [   0.    0.  -37.   -4. -250.]
 [ -25.    0.   -4.   29.  100.]]

Literal a

Para métodos directos se aplica el método de eliminación hacia adelante.

Usando el primer elemento en la diagonal.

[[  55.     0.     0.   -25.         -200.        ]
 [   0.     1.     0.     0.          -10.        ]
 [   0.     0.   -37.    -4.         -250.        ]
 [   0.     0.    -4.    17.63636364    9.09090909]]

para el segundo no es necesario, por debajo se encuentran valores cero.
Por lo que se pasa al tercer elemento de la diagonal

[[  55.     0.     0.     -25.         -200.        ]
 [   0.     1.     0.      0.          -10.        ]
 [   0.     0.   -37.     -4.         -250.        ]
 [   0.     0.     0.     18.06879607    36.11793612]]

y para el método de Gauss se emplea sustitución hacia atras.
para x₄:

18.06879607 x_4 = 36.11793612

x_4 = 1.99891216

para x₃:

-37 x_3 -4 x_3 = -250

37 x_3 = 250-4 x_4 = 250-4(1.99891216)

x_3 = 6.54065815

como ejercicio, continuar con x₁, dado que x₂=-10

55 x_1 + 25 x_4 = -200

El vector solución obtenido es:

el vector solución X es:
[[ -2.7277672 ]
 [-10.        ]
 [  6.54065815]
 [  1.99891216]]

sin embargo, para verificar la respuesta se aplica A.X=B.

[[-200.]
 [-250.]
 [ 100.]
 [ -10.]]

Se revisa el número de condición de la matriz:

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

Y para éste caso, el número de condición se encuentra alejado del valor 1, contrario a la respuesta del la primera forma de solución con la matriz 3x3. De resultar alejado del valor ideal de uno, la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.

Algoritmo en Python

Presentado por partes para revisión:

Para el método de Gauss, los resultados del algoritmo se muestran como:

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Elimina hacia adelante:
 fila 0 pivote:  55.0
   factor:  0.0  para fila:  1
   factor:  -0.45454545454545453  para fila:  2
 fila 1 pivote:  -37.0
   factor:  0.10810810810810811  para fila:  2
 fila 2 pivote:  18.06879606879607
[[  55.            0.          -25.         -200.        ]
 [   0.          -37.           -4.         -250.        ]
 [   0.            0.           18.06879607   36.11793612]]
solución: 
[-2.7277672   6.54065815  1.99891216]

Instrucciones en Python usando las funciones creadas en la unidad:

# 1Eva_IIT2019_T3 Circuito eléctrico
# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
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
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    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(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

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

def gauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss elimina hacia adelante
    tarea: verificar términos cero
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    if vertabla==True:
        print('Elimina hacia adelante:')
    for i in range(0,n,1):
        pivote = AB[i,i]
        adelante = i+1
        if vertabla==True:
            print(' fila',i,'pivote: ', pivote)
        for k in range(adelante,n,1):
            if (np.abs(pivote)>=casicero):
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    if vertabla==True:
        print(AB)
    return(AB)

def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss sustituye hacia atras
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    # Sustitución hacia atras
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    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]
        X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
    return(X)

# INGRESO
A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_sustituyeAtras(AB,vertabla=True)

# SALIDA
print('solución: ')
print(X)

literal c

Resultados con el algoritmo de Gauss Seidel

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
      diferencia,errado
0 [0. 0. 0.] 2e-05
1 [-3.6363636  6.7567568  1.2454461]
   [3.6363636 6.7567568 1.2454461] 6.756756756756757
2 [-3.0702518  6.6221139  1.7149021]
   [0.5661119 0.1346428 0.469456 ] 0.5661118513783094
3 [-2.8568627  6.5713619  1.891858 ]
   [0.2133891 0.050752  0.1769559] 0.2133891067340583
4 [-2.7764282  6.5522316  1.9585594]
   [0.0804345 0.0191304 0.0667014] 0.08043447732439457
5 [-2.7461094  6.5450206  1.9837016]
   [0.0303188 0.007211  0.0251423] 0.030318816370094925
6 [-2.7346811  6.5423025  1.9931787]
   [0.0114283 0.0027181 0.0094771] 0.011428316023939011
7 [-2.7303733  6.541278   1.996751 ]
   [0.0043078 0.0010246 0.0035723] 0.004307767346479974
8 [-2.7287495  6.5408918  1.9980975]
   [0.0016238 0.0003862 0.0013465] 0.001623761494915943
9 [-2.7281375  6.5407462  1.9986051]
   [0.0006121 0.0001456 0.0005076] 0.0006120575185017962
10 [-2.7279068  6.5406913  1.9987964]
   [2.3070778e-04 5.4871039e-05 1.9131760e-04] 0.00023070777766820427
11 [-2.7278198  6.5406707  1.9988685]
   [8.6962544e-05 2.0682983e-05 7.2114885e-05] 8.696254366213907e-05
12 [-2.727787   6.5406629  1.9988957]
   [3.2779493e-05 7.7962038e-06 2.7182845e-05] 3.277949307367578e-05
13 [-2.7277747  6.5406599  1.998906 ]
   [1.2355839e-05 2.9386860e-06 1.0246249e-05] 1.235583874370505e-05
14 [-2.72777    6.5406588  1.9989098]
   [4.6573860e-06 1.1077026e-06 3.8622013e-06] 4.6573859666665385e-06
numero de condición: 4.997509004325604
respuesta con Gauss-Seidel
[-2.72777    6.5406588  1.9989098]
>>>

Instrucciones en Python

# 1Eva_IIT2019_T3 Circuito eléctrico
# Algoritmo Gauss-Seidel
# solución de matrices
# métodos iterativos
# Referencia: Chapra 11.2, p.310,
#      Rodriguez 5.2 p.162
import numpy as np

def gauss_seidel(A,B,X0,tolera, iteramax=100,
                 vertabla=False, precision=4):
    ''' Método de Gauss Seidel, tolerancia, vector inicial X0
        para mostrar iteraciones: vertabla=True
    '''
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float)
    X0 = np.array(X0, dtype=float)
    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]')
        print('      diferencia,errado')
        print(itera, X, errado)
        np.set_printoptions(precision)
    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)
            print('  ', diferencia, errado)
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas, entrega matriz aumentada AB
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    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(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

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

# INGRESO
A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

X0  = [0.,0.,0.]

tolera = 0.00001
iteramax = 100
verdecimal = 7

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

AB = pivoteafila(A,B,vertabla=True)
# separa matriz aumentada en A y B
[n,m] = np.shape(AB)
A = AB[:,:m-1]
B = AB[:,m-1]

respuesta = gauss_seidel(A,B,X0, tolera,
                         vertabla=True, precision=verdecimal)

# SALIDA
print('numero de condición:', ncond)
print('respuesta con Gauss-Seidel')
print(respuesta)

2019-11-27

s1Eva_IIT2019_T2 Proceso Termodinámico

Ejercicio: 1Eva_IIT2019_T2 Proceso Termodinámico

la ecuación para el problema se describe como:

f(x)=e^{-0.5x}

ecuación que se usa para describir los siguientes puntos:

x	0	1	2	3	4
f(x)	1	0.60653065	0.36787944	0.22313016	0.13533528

Como el polinomio es de grado 2, se utilizan tres puntos. Para cubrir el intervalo los puntos seleccionados incluyen los extremos y el punto medio.

literal a

Con los puntos seleccionados se escriben las ecuaciones del polinomio:

p_2(x)= a_0 x^2 + a_1 x + a_2

usando los valores de la tabla:

p_2(0)=a_0 (0)^2 + a_1 (0) + a_2 = 1

p_2(2)=a_0 (2)^2 + a_1 (2) + a_2 = 0.36787944

p_2(4)=a_0 (4)^2 + a_1 (4) + a_2 = 0.13533528

con la que se escribe la matriz Vandermonde con la forma A.x=B

A= [[ 0.,  0.,  1.,]
    [ 4.,  2.,  1.,]
    [16.,  4.,  1.,]]

B= [[1.        ],
    [0.36787944],
    [0.13533528]])

matriz aumentada

[[ 0.,  0.,  1.,  1.        ]
 [ 4.,  2.,  1.,  0.36787944]
 [16.,  4.,  1.,  0.13533528]]

matriz pivoteada

[[16.,  4.,  1.,  0.13533528]
 [ 4.,  2.,  1.,  0.36787944]
 [ 0.,  0.,  1.,  1.        ]]

Resolviendo por algún método directo, la solución proporciona los coeficientes del polinomio

Tarea: escribir la solución del método directo, semejante a la presentada en el tema 3

[ 0.04994705 -0.41595438  1.        ]

con lo que el polinomio de interpolación es:

p_2(x) = 0.04994705 x^2 - 0.41595438 x + 1.0

en el enunciado se requiere la evaluación en x=2.4

p_2(2.4) = 0.04994705 (2.4)^2 - 0.41595438 (2.4) + 1.0

f(2.4)=e^{-0.5(2.4)}

error = |f(2.4)-p_2(2.4)|

Evaluando en X1:  2.4
Evaluando p(x1):  0.2894044975129779
Error en x1:      0.011789714399224216
 Error relativo:  0.039143230291095066

La diferencia entre la función y el polinomio de interpolación se puede observar en la gráfica:

literal b

Tarea: Encontrar la cota de error con f(1.7)

Algoritmo en Python

Resultado con el algoritmo

Matriz Vandermonde: 
[[ 0.  0.  1.]
 [ 4.  2.  1.]
 [16.  4.  1.]]
los coeficientes del polinomio: 
[ 0.04994705 -0.41595438  1.        ]
Polinomio de interpolación: 
0.049947050111716*x**2 - 0.415954379637711*x + 1.0

 formato pprint
                   2                            
0.049947050111716*x  - 0.415954379637711*x + 1.0

Evaluando en X1:  2.4
Evaluando p(x1):  0.2894044975129779
Error en x1:      0.011789714399224216
 Error relativo:  0.039143230291095066

Evaluando en X2:  1.7
Evaluando p(x2):  0.2894044975129779
Error en x2:      0.011789714399224216
 Error relativo:  0.039143230291095066

Presentado por secciones, semejante a lo desarrollado en clases

# 1Eva_IIT2019_T2 Proceso Termodinámico
# El polinomio de interpolación
import numpy as np
import sympy as sym

# INGRESO
fx = lambda x: np.exp(-0.5*x)
xi =np.array([0,2,4],dtype=float)

# determina vector
fi= fx(xi)

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)

B = fi
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype=float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda x:
# para evaluación con vectores de datos xin
muestras = 21
px = sym.lambdify(x,polinomio)

# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)


# literal b
x1 = 2.4
px1 = px(x1)
fx1 = fx(x1)
errorx1 = np.abs(px1-fx1)
errorx1rel = errorx1/fx1
x2 = 1.7
px2 = px(x1)
fx2 = fx(x1)
errorx2 = np.abs(px1-fx1)
errorx2rel = errorx1/fx1
print()
print('Evaluando en X1: ',x1)
print('Evaluando p(x1): ',px1)
print('Error en x1:     ',errorx1)
print(' Error relativo: ', errorx1rel)
print()
print('Evaluando en X2: ',x2)
print('Evaluando p(x2): ',px2)
print('Error en x2:     ',errorx2)
print(' Error relativo: ', errorx2rel)


# GRAFICA
import matplotlib.pyplot as plt
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)

# Usando evaluación simbólica
##yin = np.zeros(muestras,dtype=float)
##for j in range(0,muestras,1):
##    yin[j] = polinomio.subs(x,xin[j])

plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.plot(xin,fx(xin), label='f(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()

2019-11-27

1Eva_IIT2019_T3 Circuito eléctrico
1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

Tema 3. (30 puntos) El sistema de ecuaciones que sigue se generó por medio de aplicar la ley de malla de corriente al circuito de la figura.
\begin{cases} 55 I_1 - 25 I_4 = -200 \\ -37 I_3 - 4 I_4 = -250 \\ -25 I_1 - 4 I_3 + 29 I_4 = 100 \end{cases}
a) Use el método de eliminación de Gauss para calcular I₁, I₃, I₄, I₁ observando que
I₂ = -10

b) Encuentre la norma infinita de la matriz de transición T en el método de Jacobi y comente.

c) Con el método de Gauss-Seidel realice tres iteraciones comenzando con el vector cero. Además en la tercera iteración, encuentre una cota para el error relativo.

Rúbrica: literal a (12 puntos), literal b (6 puntos), literal c (12 puntos)
```
A = [[ 55.0, 0,  0, -25],
     [  0  , 0,-37,  -4],
     [-25  , 0, -4,  29],
     [  0  ,  1, 0,   0]]

B = [-200,-250,100,-10]
```
2019-11-27
1Eva_IIT2019_T4 Concentración de químico

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

Tema 4. (20 puntos) La siguiente ecuación permite calcular la concentración de un químico en un reactor, donde se tiene una mezcla completa.

$C = C_{ent} ( 1 - e^{-0.04t})+C_{0} e^{-0.03t}$

Si la concentración inicial es C₀ = 4 y la concentración de entrada es C_ent = 10, use el método de Newton-Raphson con t₀ = 0, para aproximar el tiempo requerido para que el valor de C sea 93% de C_ent.

Encuentre un intervalo en donde la convergencia está garantizada.

Rúbrica: intervalo (5 puntos), iteraciones (10 puntos), convergencia (5 puntos)

2019-11-27
s1Eva_IIT2019_T4 Concentración de químico
Ejercicio: 1Eva_IIT2019_T4 Concentración de químico

formula a usar:
C = C_{ent} ( 1 - e^{-0.04t})+C_{0} e^{-0.03t}
Se sustituyen los valores dados con:
C₀ = 4, C_ent = 10, C = 0.93 C_ent.
0.93(10) = 10 ( 1 - e^{-0.04t}) + 4 e^{-0.03t}
igualando a cero para forma estandarizada del algoritmo,
10( 1 - e^{-0.04t}) + 4 e^{-0.03t} - 9.3 = 0 0.7 - 10 e^{-0.04t} + 4 e^{-0.03t} = 0
se usas las funciones f(t) y f'(t) para el método de Newton-Raphson,
f(t) = 0.7 - 10 e^{-0.04t} + 4 e^{-0.03t} f'(t) = - 10(-0.04) e^{-0.04t} + 4(-0.03) e^{-0.03t} f'(t) = 0.4 e^{-0.04t} - 0.12 e^{-0.03t}
con lo que se pueden realizar los cálculos de forma iterativa.
t_{i+1} = t_i -\frac{f(t_i)}{f'(t_i)}
De no disponer de la gráfica de f(t), y se desconoce el valor inicial para t₀ se usa 0. Como no se indica la tolerancia, se estima en 10^-4

Iteración 1

t₀ = 0
f(0) = 0.7 - 10 e^{-0.04(0)} + 4 e^{-0.03(0)} = 5.3 f'(0) = 0.4 e^{-0.04(0)} - 0.12 e^{-0.03(0)} = -0.28 t_{1} = 0 -\frac{5.3}{-0.28} = 18.92 error = |18.92-0| = 18.92
Iteración 2
f(18.92) = 0.7 - 10 e^{-0.04(18.92)} + 4 e^{-0.03(18.92)} = -1.72308 f'(18.92) = 0.4 e^{-0.04(18.92)} - 0.12 e^{-0.03(18.92)} = 0.119593 t_{2} = 18.92 -\frac{-1.723087}{0.119593} = 33.3365 error = |33.3365 - 18.92| = 14.4079
Iteración 3
f(33.3365) = 0.7 - 10 e^{-0.04(33.3365)} + 4 e^{-0.03(33.3365)} = -0.4642 f'(33.3365) = 0.4 e^{-0.04(33.3365)} - 0.12 e^{-0.03(33.3365)} = 0.06128 t_{3} = 33.3365 -\frac{-0.46427}{-5.8013} = 40.912 error = |40.912 - 33.3365| = 7.5755
Observando que los errores disminuyen entre cada iteración, se encuentra que el método converge.

y se forma la siguiente tabla:
```
['xi' ,  'xnuevo', 'error']
[ 0.      18.9286  18.9286]
[18.9286  33.3365  14.4079]
[33.3365  40.912    7.5755]
[40.912   42.654     1.742]
[42.654   42.7316   0.0776]
[4.2732e+01 4.2732e+01 1.4632e-04]
raiz en:  42.731721341402796
```
Observando la gráfica de la función puede observar el resultado:

Algoritmo en Python
```
# 1Eva_IIT2019_T4
# Método de Newton-Raphson

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda t: 0.7-10*np.exp(-0.04*t)+4*np.exp(-0.03*t)
dfx = lambda t:0.40*np.exp(-0.04*t)-0.12*np.exp(-0.03*t)

x0 = 0
tolera = 0.001

a = 0
b = 60
muestras = 21

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

tabla = np.array(tabla)
n=len(tabla)

# para la gráfica
xk = np.linspace(a,b,muestras)
fk = fx(xk)

# SALIDA
print(['xi', 'xnuevo', 'error'])
np.set_printoptions(precision = 4)
for i in range(0,n,1):
    print(tabla[i])
print('raiz en: ', xi)

# grafica
plt.plot(xk,fk)
plt.axhline(0, color='black')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.show()
```
2019-11-27
3Eva_IT2019_T3 EDP Difusión en sólidos

3ra Evaluación I Término 2019-2020. 10/Septiembre/2019. MATG1013

Tema 3. (30 Puntos). En el año 1855, los experimentos de Adolf Fick tratan sobre la medición de concentraciones y sus flujos, también ahora aplicados a la difusión en sólidos que en ese tiempo no se consideraba posible.

La gráfica muestra los cambios en el tiempo de concentración Φ de un gas en un sólido (estado no-estacionario) para un sólido semi infinito (eje y).

La segunda ley de Fick predice la forma en que la difusión causa que la concentración cambie con el tiempo. Se trata de una ecuación diferencial parcial que en una dimensión se escribe:
$\frac{\partial \phi}{\partial t} = D\frac{\partial ^2 \phi}{\partial x^2}$

Φ(0, t) = 5
Φ(L, t) = 0
Φ(x,0) = 0 D = 0.16
L =0.1

a. Plantee las ecuaciones, la malla, desarrolle y obtenga el modelo Φ(x_i,t_j)

b. Aproxime la solución con Δx = 0.02, Δt = Δx/100. Realice al menos tres iteraciones en el eje tiempo.

c. Estime el error de Φ(x_i,t_j)

Rúbrica: Construir la malla (5 puntos), plantear la ecuación en el nodo i,j (5 puntos), modelo de ecuación (5 puntos), literal b (10 puntos), literal c (5 puntos).

Referencia: https://es.wikipedia.org/wiki/Leyes_de_Fick;
Difusión 2ª Ley de Fick|7/22|UPV (2011) https://www.youtube.com/watch?v=HHBvZDNvTic

2019-09-10
s3Eva_IT2019_T3 EDP Difusión en sólidos
Ejercicio: 3Eva_IT2019_T3 EDP Difusión en sólidos

Siguiendo el procedimiento planteado en la sección EDP parabólicas, se plantea la malla del ejercicio:

Para plantear la ecuación en forma discreta:
\frac{\phi_{i,j+1}-\phi_{i,j}}{\Delta t}=D\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{(\Delta x)^2}
y resolver usando el método explícito para ecuaciones parabólicas, obteniendo el siguiente resultado:
\phi_{i,j+1}-\phi_{i,j}=D\frac{\Delta t }{\Delta x^2}(\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \lambda = D\frac{\Delta t }{\Delta x^2} \phi_{i,j+1}-\phi_{i,j}=\lambda (\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \phi_{i,j+1} =\lambda \phi_{i+1,j}-2\lambda\phi_{i,j}+\lambda\phi_{i-1,j}+\phi_{i,j} \phi_{i,j+1} =\lambda \phi_{i+1,j}(1-2\lambda)\phi_{i,j}+\lambda\phi_{i-1,j} \phi_{i,j+1} =P \phi_{i+1,j}+Q\phi_{i,j}+R\phi_{i-1,j}
siendo:
P = λ = 0.16 (Δx/100)/Δx² = 0.0016/Δx = 0.0016/0.02=0.08
Q = 1-2λ = 1-2*(0.08) = 0.84
R = λ =0.08
\phi_{i,j+1} =0.08 \phi_{i+1,j}+ 0.84\phi_{i,j}+0.08\phi_{i-1,j}
Iteración 1 en tiempo:
i=1, j=0
\phi_{1,1} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0} \phi_{1,1} =0.08 (0)+ 0.84(0)+0.08(5)=0.4
i=2,j=0
\phi_{2,1} =0.08 \phi_{3,0}+ 0.84\phi_{2,0}+0.08\phi_{1,0} = 0
Para los proximos valores i>2, todos los resultados son 0

Iteración 2 en tiempo
i=1, j=1
\phi_{1,2} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0}
$\phi_{1,2} =0.08 (0)+ 0.84(0.4)+0.08(5)=0.736$
i=2, j=1
\phi_{2,2} =0.08 \phi_{3,1}+ 0.84\phi_{2,1}+0.08\phi_{1,1} \phi_{2,2} =0.08(0)+ 0.84(0)+0.08(0.4) = 0.032
i=3, j=1
\phi_{3,2} =0.08\phi_{4,1}+ 0.84\phi_{3,1}+0.08\phi_{2,1}=0
Para los proximos valores i>3, todos los resultados son 0

Tarea: Desarrollar la iteración 3 en el tiempo.

siguiendo las iteraciones se tiene la siguiente tabla:
```
[[5.0, 0.000, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.400, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.736, 0.032, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.021, 0.085, 0.00256, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.264, 0.153, 0.00901, 0.00020, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
...
]
```
Con lo que se obtiene la siguiente gráfica.

El resultado se interpreta mejor con una animación: (tarea)

Tarea: Presentar el orden de error de la ecuación basado en las fórmulas de diferenciación

Algorirmo en Python
```
# 3Eva_IT2019_T3 Difusión en sólidos
# EDP parabólicas. método explícito,usando diferencias finitas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 5
Tb = 0
T0 = 0
# longitud en x
a = 0
b = 0.1
# Constante K
K = 1/(1.6e-1)
# Tamaño de paso
dx = 0.02
dt = dx/100
# iteraciones en tiempo
n = 50

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0
u[ultimox,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
    u[0,j+1] = Ta
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    u[m-1,j+1] = Tb
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('phi[i,j]')
plt.title('Solución EDP parabólica')
plt.show()
```
La animación se complementa con lo mostrado en la sección de Unidades.
2019-09-10

x	0	1	2	3	4
f(x)	f(0)	f(1)	f(2)	f(3)	f(4)

Autor: Edison Del Rosario

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

literal a y b

Literal c

literal e

Literal a

literal b

Literal c

Solución alterna

Literal a

Algoritmo en Python

literal c

literal a

literal b

Algoritmo en Python

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

1ra Evaluación II Término 2019-2020. 26/Noviembre/2019. MATG1013

Algoritmo en Python

3ra Evaluación I Término 2019-2020. 10/Septiembre/2019. MATG1013

Algorirmo en Python