Ejercicio: 1Eva_IIT2008_T1 Distribuidores de productos

Siguiendo las instrucciones del enunciado, el promedio de precios del nodo A, se conforma de los precios en los nodos aledaños menos el costo de transporte.

precio en X₁ para A = precio en nodo_A – costo de transporte_A

siguiendo el mismo procedimiento,

precio en X₁ para A: (3.1-0.2)
precio de X₁ para B: (2.8-0.3)
precio de X₁ para C: (2.7-0.4)
precio de X₁ para X₂: (X₂-0.1)
precio de X₁ para X₃: (X₃-0.5)

x_1 = \frac{1}{5} \Big[ (3.1-0.2)+(2.8-0.3)+(2.7-0.4)+ +(x_2-0.1)+(x_3-0.5)\Big] x_1 = \frac{1}{5} \Big[ 2.9+2.5 +2.3+x_2+x_3-0.6\Big] x_1 = \frac{1}{5} (7.1+x_2+x_3) 5x_1 = 7.1+x_2+x_3 5x_1-x_2-x_3 = 7.1

Se continua con el mismo proceso para los siguientes nodos:

x_2 = \frac{1}{4} \Big[ (3.2-0.5)+(3.4-0.3) +(x_1-0.1)+(x_3-0.2)\Big] 4x_2 = (3.2-0.5)+(3.4-0.3) +(x_1-0.1)+(x_3-0.2) 4x_2 = 2.7+3.1 +x_1+x_3-0.3 -x_1+4x_2-x_3 = 5.5

Para X₃

x_3 = \frac{1}{4} \Big[ (3.3-0.3)+(2.9-0.2) +(x_1-0.5)+(x_2-0.2)\Big] 4x_3 = 3.0+2.7+x_1+x_2-0.7 4x_3 = 5+x_1+x_2 -x_1-x_2+4x_3= 5

El sistema de ecuaciones se convierte en:

\begin{cases} 5x_1-x_2-x_3 = 7.1 \\ -x_1+4x_2-x_3 = 5.5 \\-x_1-x_2+4x_3= 5 \end{cases}

Para resolve se convierte a la forma Ax=B

\begin{bmatrix} 5 && -1 && -1 \\ -1 && 4 && -1 \\ -1 && -1 && 4 \end{bmatrix} . \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 7.1 \\ 5.5 \\5 \end{bmatrix}

En los métodos directos se usa la forma de matriz aumentada

\begin{bmatrix} 5 && -1 && -1 && 7.1 \\ -1 && 4 && -1 && 5.5 \\ -1 && -1 && 4 && 5 \end{bmatrix}

pivoteo: no es necesario, pues la matriz ya está ordenada de forma diagonalmente dominante.

Eliminación hacia adelante
1ra Iteración

\begin{bmatrix} 5 && -1 && -1 && 7.1 \\ -1-5\big(\frac{-1}{5}\big) && 4 -(-1)\big(\frac{-1}{5}\big) && -1 -(-1)\big(\frac{-1}{5}\big) && 5.5-7.1\big(\frac{-1}{5}\big) \\ -1-5\big(\frac{-1}{5}\big) && -1-(-1)1\big(\frac{-1}{5}\big) && 4-(-1)\big(\frac{-1}{5}\big) && 5-7.1\big(\frac{-1}{5}\big) \end{bmatrix} \begin{bmatrix} 5 && -1 && -1 && 7.1 \\ 0 && 3.8 && -1.2 && 6.92 \\ 0 && -1.2 && 3.8 && 6.42 \end{bmatrix}

Eliminación hacia adelante
2da Iteración

Elimina hacia adelante
[[ 5.         -1.         -1.          7.1       ]
 [ 0.          3.8        -1.2         6.92      ]
 [ 0.          0.          3.42105263  8.60526316]]

Eliminación hacia atras, continuando el desarrollo de forma semejante a los pasos anteriores se obtiene:

Elimina hacia atras
[[ 1.         -0.         -0.          2.44615385]
 [ 0.          1.         -0.          2.61538462]
 [ 0.          0.          1.          2.51538462]]
el vector solución X es:
[[2.44615385]
 [2.61538462]
 [2.51538462]]
verificar que A.X = B
[[7.1]
 [5.5]
 [5. ]]

---

Para el literal b se usa como referencia el número de condición:

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

El número de condición es cercano a 1, dado que la matriz A es diagonalmente dominante pues los valores mayores de la fila se encuentran en la diagonal. Como el número de condición es cercano a 1 el sistema converge usando métodos iterativos.

La selección del vector inicial para las iteraciones siguiendo el enunciado del problema, se evita el vector cero dado que el precio de un producto para una fabrica no puede ser cero. Se observa los valores de los precios, y se encuentra que el rango de existencia en los nodos es [ 2.7, 3.4] que restando el valor de transporte podrían ser un valor menor a 2.7. Por lo que un buen vector inicial será [2,2,2]

Queda como tarea las iteraciones por un método Jacobi.

Ejercicio: 1Eva_IIT2008_T1_MN Bacterias contaminantes

a.  Planteamiento del problema

estandarizar la fórmula a la forma f(x) = 0

c = 70 e^{-1.5t} + 25 e^{-0.075t}

el valor que debe tomar c = 9, por lo que la función a desarrollar se convierte en:

9 = 70 e^{-1.5t} + 25 e^{-0.075t}

y la que se usará en el algoritmo de búsqueda de raices es:

70 e^{-1.5t} + 25 e^{-0.075t} -9 = 0 f(t) = 70 e^{-1.5t} + 25 e^{-0.075t} -9 = 0

---

b. Intervalo de búsqueda [a,b]

Como la variable t representa el tiempo, se inicia el análisis con cero por la izquierda o un poco mayor. a=0

Para verificar que existe raíz se evalua f(a) y f(b) buscando valores donde se presenten cambios de signo.

La función depende de tiempo, por lo que a = 0, y b seleccionar un valor mayor.

Para el caso de b,  si fuesen minutos, se pueden probar con valores consideranto el tiempo que duraría un «experimento» en el laboratorio durante una clase. Por ejemplo 60 minutos, 30 minutos, etc, para obtener un punto «b» donde se garantice un cambio de signo f(a) y f(b)

Para el ejemplo, se prueba con b = 15

a = 0
b = 15
f(a) = f(0) = 70*e0 + 25 e0 -9 = 86
f(b) = f(15) = 70*e-1.5*15. + 25 e-0.075*15 -9 = - 0.0036

---

c. Número de iteraciones

Calcular el número de iteraciones para llegar a la raíz con el error tolerado

error = 0.001

\frac{|15-0|}{2^n} = 0.001 15/0.001 = 2^n log(15/0.001) = nlog(2) n = \frac{log(15/0.001)}{log(2)} = 13.87 n = 14

Verificando la selección usando una gráfica, usando 50 tramos entre [a,b] o 51 muestras en el rango.

---

d. Calcular la raíz:

Usando el algoritmo se encuentra que la raiz está en:

raiz en:  13.62213134765625
f(raiz) = -7.733799068e-05

si realiza la tabla:

[i,a,b,c,sign(fa),sign(fb),sign(fc),paso]
[[1, 0, 15, 7.5, 1.0, -1.0, 1.0, 15], 
 [2, 7.5, 15, 11.25, 1.0, -1.0, 1.0, 7.5], 
 [3, 11.25, 15, 13.125, 1.0, -1.0, 1.0, 3.75], 
 [4, 13.125, 15, 14.0625, 1.0, -1.0, -1.0, 1.875], 
 [5, 13.125, 14.0625, 13.59375, 1.0, -1.0, 1.0, 0.9375], 
 [6, 13.59375, 14.0625, 13.828125, 1.0, -1.0, -1.0, 0.46875], 
 [7, 13.59375, 13.828125, 13.7109375, 1.0, -1.0, -1.0, 0.234375], 
 [8, 13.59375, 13.7109375, 13.65234375, 1.0, -1.0, -1.0, 0.1171875], 
 [9, 13.59375, 13.65234375, 13.623046875, 1.0, -1.0, -1.0, 0.05859375], 
 [10, 13.59375, 13.623046875, 13.6083984375, 1.0, -1.0, 1.0, 0.029296875],
 [11, 13.6083984375, 13.623046875, 13.61572265625, 1.0, -1.0, 1.0, 0.0146484375], 
 [12, 13.61572265625, 13.623046875, 13.619384765625, 1.0, -1.0, 1.0, 0.00732421875], 
 [13, 13.619384765625, 13.623046875, 13.6212158203125, 1.0, -1.0, 1.0, 0.003662109375], 
 [14, 13.6212158203125, 13.623046875, 13.62213134765625, 1.0, -1.0, -1.0, 0.0018310546875]]

---

Cálculos con  Newton –  Raphson

Se encuentra la derivada f'(t) y se aplica el algoritmo Newton-Raphson con valor inicial cero.

f'(t) = 70(-1.5) e^{-1.5t} + 25(-0.075) e^{-0.075t}

raiz en:  13.622016772385583

Se usa el algoritmo en python para encontrar el valor. El algoritmo newton Raphson mostrado es más pequeño que por ejemplo la bisección, pero requiere realizar un trabajo previo para encontrar la derivada de la función.

# 1ra Eval II Termino 2008 tema 1
# Método de Newton-Raphson
import numpy as np

def newtonraphson(funcionx, fxderiva, c, tolera):
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = c - funcionx(c)/fxderiva(c)
        tramo = abs(xnuevo-c)
        c = xnuevo
    return(c)


# PROGRAMA ---------------------------------
funcionx = lambda t: 70*(np.e**(-1.5*t)) +25*(np.e**(-0.075*t))-9
fxderiva = lambda t: -70*1.5*(np.e**(-1.5*t)) -25*0.075*(np.e**(-0.075*t))

# INGRESO
c = 0.1
tolera = 0.001

# PROCEDIMIENTO
raiz = newtonraphson(funcionx, fxderiva, c, tolera)

# SALIDA
print('raiz en: ', raiz)

# Gráfica
import matplotlib.pyplot as plt
xi = np.linspace(0,15,100)
yi = funcionx(xi)
plt.plot(xi,yi)
plt.axhline(0, color = 'k')
plt.show()

Tarea: Para el problema, realice varios métodos y compare el número de iteraciones y el trabajo realizado al plantear el problema al implementar cada uno.

Ejercicio: 1Eva_IT2008_T1 Raíz de funcion(f)

Pasos a Seguir usando: BISECCION

1. Plantear la fórmula estandarizada f(x) = 0
2. Seleccionar el rango de análisis [a,b] donde exista cambio de signo.
3. Calcular el número de iteraciones para llegar a la raíz con el error tolerado
4. Calcular la raíz:
   4.1 Solución manual en papel: Realizar la tabla que muestre las iteraciones del método:
   4.2 Solución usando el algoritmo: Usar el algoritmo para encontrar la raiz.

---

1. Plantear la fórmula estandarizada f(x) = 0

\sqrt{f} . ln \Big( R\frac{\sqrt{f}}{2.51} \Big) - 1.1513 = 0

La función depende de la variable f, por lo que por facilidad de trabajo se cambiará a x para no ser confundida con una función f(x).
El valor de R también se reemplaza con R=5000 como indica el problema.
El error tolerado para el problema es de 0.00001

\sqrt{x} . ln \Big( 5000\frac{\sqrt{x}}{2.51} \Big) - 1.1513 = 0

2. Seleccionar el rango de análisis [a,b] donde exista cambio de signo

La función tiene un logaritmo, por lo que no será posible iniciar con cero, sin o con  valor un poco mayor a=0.01. Para el límite superior se escoge para prueba b=2. y se valida el cambio de signo en la función.

>>> import numpy as np
>>> fx = lambda x: np.sqrt(x)*np.log((5000/2.51)*np.sqrt(x))-1.1513
>>> fx(0.01)
-0.62186746547214999
>>> fx(2)
10.082482845673042

validando el rango por cambio de signo en la función [0.01 ,2]

3. Calcular el número de iteraciones para llegar a la raíz con el error tolerado

error = 0.00001

\frac{|2-0.01|}{2^n} = 0.001 1.99/0.00001 = 2^n log(1.99/0.00001) = nlog(2) n = \frac{log(1.99/0.00001)}{log(2)} = 17.6 n = 18

4. Calcular la raíz:

Usando el algoritmo se encuentra que la raiz está en:

raiz:  0.0373930168152
f(raiz) =  -2.25294254252e-06

Primeras iteraciones de la tabla resultante:

 
a,b,c, fa, fb, fc, error
[[  1.00000000e-02   2.00000000e+00   1.00500000e+00  -6.21867465e-01
    6.46707903e+00   6.46707903e+00   1.99000000e+00]
 [  1.00000000e-02   1.00500000e+00   5.07500000e-01  -6.21867465e-01
    4.01907320e+00   4.01907320e+00   9.95000000e-01]
 [  1.00000000e-02   5.07500000e-01   2.58750000e-01  -6.21867465e-01
    2.36921961e+00   2.36921961e+00   4.97500000e-01]
 [  1.00000000e-02   2.58750000e-01   1.34375000e-01  -6.21867465e-01
    1.26563722e+00   1.26563722e+00   2.48750000e-01]
 [  1.00000000e-02   1.34375000e-01   7.21875000e-02  -6.21867465e-01
    5.36709904e-01   5.36709904e-01   1.24375000e-01]
 [  1.00000000e-02   7.21875000e-02   4.10937500e-02  -6.21867465e-01
    6.51903790e-02   6.51903790e-02   6.21875000e-02]
...

el número de iteraciones es filas-1 de la tabla

>>> len(tabla)
19

con lo que se comprueba los resultados.

Ejercicio: 1Eva_IIIT2007_T3 Factorar polinomio

P_3(x) = 2x^3-5x^2 + 3x-0.1

la gráfica se obtuvo con Python, con lo que se puede observar la primera raíz…

# 1ra Eval III Término 2007
# Tema 3. Factorar polinomio
import numpy as np
import matplotlib.pyplot as plt
def newtonraphson(funcionx, fxderiva, c, tolera):
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = c - funcionx(c)/fxderiva(c)
        tramo = abs(xnuevo-c)
        c = xnuevo
    return(c)

# Literal a)
p3 = lambda x: 2*x**3 - 5*x**2 + 3*x -0.1
dp3 = lambda x: 6*x**2 - 10*x +3
a = 0
b = 2
pasos = 100
c = 2
tolera = 0.0001

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

raiz1 = newtonraphson(p3, dp3, c, tolera)

# SALIDA
print('primera raiz: ',raiz1)
plt.plot(xi,p_i)
plt.axhline(0)
plt.show()

para el literal b)
se añade:

# Literal b)
# PROCEDIMIENTO
p2 =  lambda x: (2*x**3 - 5*x**2 + 3*x -0.1)/(x-raiz1)
# SALIDA
pol2 = p2(xi)
plt.plot(xi,pol2)
plt.show()

Ejercicio: 1Eva_IIIT2007_T2 Función Cobb-Douglas

usando los algoritmos para encontrar los polinomios de lagrange y búsqueda de raíz por la bisección, se obtiene:

 **** literal a:
Polinomio de Lagrange
[0, 6.16333333333368e-5*x**3 - 0.00527800000000012*x**2 + 0.254746666666668*x + 11.0, 0.00010365*x**3 - 0.008876*x**2 + 0.428425*x + 18.4997, 0.000140483333333337*x**3 - 0.0120304999999998*x**2 + 0.580686666666667*x + 25.0746]
Puntos f[25]:
[0, 15.0329375000000, 25.2823562500000, 34.2677562500002]
Polinoio con f[25]:
0.000586583333333373*x**3 - 0.0415150937500013*x**2 + 1.85978635416668*x
Estimado de [la,ka]:
29.7130898437501
 **** Literal b
La raiz se encuentra en:  25.3295898438

Revisar el algoritmo por partes: literal a y literal b.

# Interpolacion de Lagrange
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

def lagrange(xi,fi):
    n = len(xi)
    x = sp.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if j!=i:
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*fi[i]
    # Expande el polinomio
    px = polinomio.expand()
    return(px)

# INGRESO , Datos de prueba
M = np.array([[0, 0, 0, 0],
              [11.0000, 13.0813, 14.4768, 15.5563],
              [18.4997, 22.0000, 24.3470, 26.1626],
              [25.0746, 29.8189, 33.0000, 35.4608]])
li = np.array([0, 10, 20, 30])
kj = np.array([10, 20, 30, 40])

la = 25
ka = 25

# PROCEDIMIENTO
tamano = np.shape(M)
n = tamano[0]
m = tamano[1]
x = sp.Symbol('x')
poli = []
for i in range(0,n,1):
    xi = li
    fi = M[i,:]
    p = lagrange(xi,fi)
    poli.append(p)
# literal a, evalua en la
f_la =[]
for i in range(0,n,1):
    puntola = poli[i].subs(x,la)
    f_la.append(puntola)
poli_la = lagrange(li,f_la)
# literal a, evalua en ka
puntolaka = poli_la.subs(x,ka)

# SALIDA
print(' **** literal a:')
print('Polinomio de Lagrange')
print(poli)
print('Puntos f['+str(la)+']:')
print(f_la)
print('Polinoio con f['+str(ka)+']:')
print(poli_la)
print('Estimado de [la,ka]:')
print(puntolaka)


# ------------------------------------
# literal b, usa f_ka con resultado 30

def biseccion(funcionx,a,b,tolera,iteramax = 20):
    '''
    Algoritmo de Bisección
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    error = tolera
    '''
    fa = funcionx(a)
    fb = funcionx(b)

    itera = 0
    tramo = np.abs(b-a)
    while (tramo>=tolera and itera<=iteramax):
        c = (a+b)/2
        fc = funcionx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if (cambia<0):
            b = c
            fb = fc
        else:
            a = c
            fa = fc
        tramo = np.abs(b-a)
        itera = itera + 1
    respuesta = c
    # Valida respuesta
    if (i>=iteramax):
        respuesta = np.nan
    return(respuesta)


# INGRESO
a = kj[0]
b = kj[m-1]
tolera = 0.01

# PROCEDIMIENTO
# Busca raiz
fx = sp.lambdify(x,poli_la - 30)
fxi = fx(kj)

raiz = biseccion(fx,a,b,tolera)

# SALIDA
print(' **** Literal b')
print('La raiz se encuentra en: ',raiz)

# GRAFICA
plt.plot(kj,fxi)
plt.axhline(0)
plt.show()

Ejercicio: 1Eva_IIIT2007_T1 Container: Cocinas y Refrigeradoras

x: cantidad de refrigeradoras
y: cantidad de cocinas

ecuaciones:

200 x + 100  y = 1000
  2 x + 1.05 y =   10.4

literal a)
Nota: se han omitido los pasos para la solución del sistema, usando la función gauss desarrollada en python.

Si las matrices son:
A = np.array([[200, 100],
              [2,1.05]])

B = np.array([[1000],
              [10.4]])

el vector solución X es:
[ 1.  8.]
verificar que A.X sea igual a B
[ 1000.     10.4]

literal b)

>>> 
el vector solución X es:
[ 3.  4.]
verificar que A.X sea igual a B
[ 1000.     10.4]
>>> 

literal c)

Observación: el pequeño cambio de volumen de la cocina no es consistente con los resultados.

El asunto es que la forma de la refrigeradora o cocina no se adapta al volumen disponible, pues son objetos rígidos. Por lo que el sistema de ecuaciones estaría mal planteado.
Las ecuaciones tendrían sentido si esta diseñando el mejor «tamaño» para que entren la mayor cantidad dentro de un container, sin embargo los tamaños de las refrigeradoras y cocinas se encuentran estandarizados.

Revisamos el número de condición, que resulta ser del orden de 10⁴, lo que confirma que el sistema está mal condicionado.

Usando la el valor de 1.05

>>> C = np.concatenate((A,B),axis=1)
>>> C
array([[  200. ,   100. ,  1000.  ],
       [    2. ,     1.05,    10.4]])
>>> np.linalg.cond(C)
12926.000640466344

Ejercicio: 1Eva_IIT2007_T3 Interpolación inversa

f(0.50) = 1.648
f(0.65) = 1.915
f( x  ) = 2.117
f(0.80) = 2.225
f(0.95) = 2.5857

Para el algoritmo se intercambian las variables previo a usarlo.
Luego se evalua en el punto buscado, en éste caso: fi=2.117, obteniendo que x es: 0.750321134121361

fi = np.array([0.50 , 0.65 , 0.80,  0.95   ])
x = np.array([1.648, 1.915, 2.225, 2.5857 ])

Polinomio de Lagrange
0.924124055152463*(-3.74531835205992*x + 7.17228464419475)*(x - 2.5857)*(x - 2.225) + 3.12624749298999*(x - 2.5857)*(x - 2.225)*(3.74531835205992*x - 6.17228464419475) - 7.15454716188057*(x - 2.5857)*(x - 1.915)*(1.73310225303293*x - 2.85615251299827) + 3.92689380344011*(x - 2.225)*(x - 1.915)*(1.06643915964594*x - 1.75749173509651)
Expandiendo: 
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226
>>> px
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226
>>> px.subs(x,2.117)
0.750321134121361

El algortmo modificado usado es:

# Interpolacion de Lagrange
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

# INGRESO , Datos de prueba
fi = np.array([0.50 , 0.65 , 0.80,  0.95   ])
xi = np.array([1.648, 1.915, 2.225, 2.5857 ])

# PROCEDIMIENTO
n = len(xi)
x = sp.Symbol('x')
# Polinomio
polinomio = 0
for i in range(0,n,1):
    # Termino de Lagrange
    termino = 1
    for j  in range(0,n,1):
        if j!=i:
            termino = termino*(x-xi[j])/(xi[i]-xi[j])
    polinomio = polinomio + termino*fi[i]
    
# Expande el polinomio
px = polinomio.expand()

# Salida
print('Polinomio de Lagrange')
print(polinomio)
print('Expandiendo: ')
print(px)

para visualizar el resultado, intercambiando nuevamente los ejes:

>>> px
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226
>>> px.subs(x,2.117)
0.750321134121361
>>> px
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226

para revisar con la gráfica, se añaden las líneas,

# GRAFICA
polinomio = sp.lambdify(x,px)
y = np.linspace(1,3,100)
pyi= polinomio(y)
plt.plot(pyi,y)
plt.show()

Ejercicio: 1Eva_IIT2007_T2 Aplicar Gauss-Seidel

Probando solución con Jacobi, enviado como tarea

Desarrrollo 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.0451853966291011

– realizar iteraciones

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.]

iteración 1

x0 = (-9.44 -0.30(0) -0.15(0) 
      -0.50(0) -0.34(0) -0.84(0))/7.63 
   = -9.44/7.63 = -1.23722149
x1 = (25.27 -0.38(0) -0.70(0) 
      -0.90(0) -0.29(0) -0.57(0))/6.40 
   = 25.27/6.40 = 3.9484375
x2 = (-48.01 -0.83(0) -0.19(0) -0.82(0)
      -0.34(0) -0.37(0))/8.33 
   = -48.01/8.33 = -5.7635054
x3 = (19.76 -0.50(0) -0.68(0) -0.86(0)
      -0.53(0) -0.70(0))/10.21 
   = 19.76/10.21 = 1.93535749
x4 = (-23.63 - 0.71(0) -0.30(0) -0.85(0)
      -0.82(0) -0.55(0))/5.95 
   = -23.63/5.95 = -3.97142857
x5 = (62.59 - 0.43(0) -0.54(0) -0.59(0)
      -0.66(0) -0.31(0))/9.25 
   = 62.59/9.25 = 6.76648649
X1 = [-1.23722149  3.9484375  -5.7635054   1.93535749 -3.97142857  6.76648649]
diferencia = X1 - X0 = X1
errado = max(|diferencia|) = 6.76648649

iteración 2

x0 = (-9.44 -0.30(3.9484375) -0.15(-5.7635054) 
      -0.50(1.93535749) -0.34(-3.97142857) -0.84(6.76648649))/7.63 
   = -9.44/7.63 = -1.23722149
x1 = (25.27 -0.38(-1.23722149) -0.70(-5.7635054) 
      -0.90(1.93535749) -0.29(-3.97142857) -0.57(6.76648649))/6.40 
   = 25.27/6.40 = 3.9484375
x2 = (-48.01 -0.83(-1.23722149) -0.19(3.9484375)
      -0.82(1.93535749) -0.34(-3.97142857) -0.37(6.76648649))/8.33 
   = -48.01/8.33 = -5.7635054
x3 = (19.76 -0.50(-1.23722149) -0.68(3.9484375)
      -0.86(-5.7635054) -0.53(-3.97142857) -0.70(6.76648649))/10.21 
   = 19.76/10.21 = 1.93535749
x4 = (-23.63 - 0.71(-1.23722149) -0.30(3.9484375) 
      -0.85(-5.7635054) -0.82(1.93535749) -0.55(6.76648649))/5.95 
   = -23.63/5.95 = -3.97142857
x5 = (62.59 - 0.43(-1.23722149) -0.54(3.9484375) 
      -0.59(-5.7635054) -0.66(1.93535749) -0.31(-3.97142857))/9.25 
   = 62.59/9.25 = 6.76648649
X1 = [-1.97395113  3.95743644 -6.05925771  1.96068604 -4.09171178  6.95612152]
diferencia = X1 - X0 = [-0.73672964,  0.00899894, -0.29575231,  0.02532855, -0.12028321, 0.18963504]
errado = max(|diferencia|) = 0.736729635697

iteración 3

x0 = (-9.44 -0.30(3.95743644) -0.15(-6.05925771) 
      -0.50(1.96068604) -0.34(-4.09171178) -0.84(6.95612152))/7.63 
   = -9.44/7.63 = -1.23722149
x1 = (25.27 -0.38(-1.97395113) -0.70(-6.05925771) 
      -0.90(1.96068604) -0.29(-4.09171178) -0.57(6.95612152))/6.40 
   = 25.27/6.40 = 3.9484375
x2 = (-48.01 -0.83(-1.97395113) -0.19(3.95743644) 
      -0.82(1.96068604) -0.34(-4.09171178) -0.37(6.95612152))/8.33 
   = -48.01/8.33 = -5.7635054
x3 = (19.76 -0.50(-1.97395113) -0.68(3.95743644) 
     -0.86(-6.05925771) -0.53(-4.09171178) -0.70(6.95612152))/10.21 
   = 19.76/10.21 = 1.93535749
x4 = (-23.63 - 0.71(-1.97395113) -0.30(3.95743644) 
      -0.85(-6.05925771) -0.82(1.96068604) -0.55(6.95612152))/5.95 
   = -23.63/5.95 = -3.97142857
x5 = (62.59 - 0.43(-1.97395113) -0.54(3.95743644) 
      -0.59(-6.059257710) -0.66(1.96068604) -0.31(-4.09171178))/9.25 
   = 62.59/9.25 = 6.76648649
X1 = [-1.98566781  4.0185268  -5.9920623   2.01431955 -3.98302285  7.01093224]
diferencia = X1 - X0 = [-0.01171668,  0.06109037,  0.0671954 ,  0.05363351,  0.10868893, 0.05481072]
errado = max(|diferencia|) = 0.108688931048

Desarrollo numérico con Python

Se verifica el resultado obtenido realizando A.Xi y comparando con el vecto B
en la tabla se usa el signo de errado para la gráfica.

X0:  [ 0.  0.  0.  0.  0.  0.]
Xi:  [-1.23722149  3.9484375  -5.7635054   1.93535749 -3.97142857  6.76648649]
errado:  6.76648648649
Xi:  [-1.97395113  3.95743644 -6.05925771  1.96068604 -4.09171178  6.95612152]
errado:  -0.736729635697
Xi:  [-1.98566781  4.0185268  -5.9920623   2.01431955 -3.98302285  7.01093224]
errado:  0.108688931048
Xi:  [-2.00378293  3.99452422 -6.00443878  1.99576482 -4.0067623   6.9961552 ]
errado:  -0.0240025781157
Xi:  [-1.99869529  4.00195452 -5.99863447  2.00153846 -3.99769932  7.00130746]
errado:  0.00906298532211
Xi:  [-2.00045097  3.99933614 -6.00047801  1.99948184 -4.00078219  6.99955127]
errado:  -0.00308287453133
Xi:  [-1.99984629  4.00022733 -5.99983706  2.00017793 -3.99973154  7.00015339]
errado:  0.0010506528253
Xi:  [-2.00005265  3.9999222  -6.00005579  1.99993915 -4.00009178  6.9999475 ]
errado:  -0.0003602430897
Xi:  [-1.99998199  4.00002662 -5.99998091  2.00002082 -3.99996859  7.00001796]
errado:  0.000123195668206
Xi:  [-2.00000616  3.99999089 -6.00000653  1.99999287 -4.00001075  6.99999385]
errado:  -4.21623029112e-05

respuesta de A.X=B : 
[-2.00000616  3.99999089 -6.00000653  1.99999287 -4.00001075  6.99999385]
iteraciones:  10

A.Xi:  [ -9.44006312  25.26992175 -48.01007303  19.75990236 -23.63008584
  62.58992368]
   B: [ -9.44 25.27 -48.01 19.76 -23.63 62.59]

el gráfico de los errores vs iteraciones es:

---

Aplicando Gauss-Seidel

la tabla de aproximaciones sucesivas para el vector X es:

Tabla de iteraciones con AX=B: 
[[ 0.          0.          0.          0.          0.          0.        ]
 [-1.23722149  4.02189753 -5.73196479  2.21089228 -3.51242077  6.91477852]
 [-2.0322951   3.92844105 -6.03202587  1.98957766 -3.97864919  7.00774962]
 [-1.9976784   4.00317297 -6.00049341  1.99807691 -4.00081785  6.99990294]
 [-1.99994191  4.00036665 -5.99978715  2.00000392 -4.00004739  6.99996363]
 [-2.00001274  3.99998231 -5.99999516  2.00000635 -3.99999579  7.00000072]
 [-2.00000008  3.99999833 -6.00000078  1.99999991 -3.99999985  7.00000015]
 [-1.99999994  4.00000007 -6.00000001  1.99999997 -4.00000002  7.        ]]
>>> 

que se obtiene aplicando la función de Gauss-Seidel, tomando como vector inicial el vector de ceros.
Tarea: X=TX+C

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

def gauss_seidel_tabla(A,B,tolera,X = [0], itermax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    if (len(X)==n):
        X = np.zeros(n, dtype=float)
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    tabla = [np.copy(X)]
        
    itera = 0
    while (errado>tolera or itera>itermax):
        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
        tabla.append(np.copy(X))

    tabla = np.array(tabla)
        
    # No converge
    if (itera>itermax):
        X=0
    return(tabla)

# 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

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

# SALIDA
print('Tabla de iteraciones con AX=B: ')
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

Con lo que la fórmula para el objetivo 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

para buscar la raíz:

f(p) = (1-p)^{5} + 5p (1-p)^{4} - 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
>>> 

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

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

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 = 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 Newton-Raphson…

Verificando el polinomio a usar con python:

>>> 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

Tarea: con lo que puede continuar con el procedimiento para encontrar la raiz…

verificar con raiz: 0.3649852264049102

---

Newton-Raphson

# Método de Newton-Raphson
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
import numpy as np

def newtonraphson(funcionx, fxderiva, c, tolera):
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = c - funcionx(c)/fxderiva(c)
        tramo = abs(xnuevo-c)
        c = xnuevo
    return(c)


# PROGRAMA #######################
import matplotlib.pyplot as plt

funcionx = lambda p: 4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
fxderiva = lambda p: 20*p**4 - 60*p**3 + 60*p**2 - 20*p

# INGRESO
c = 0.2
tolera = 0.0000001

# PROCEDIMIENTO
raiz = newtonraphson(funcionx, fxderiva, c, tolera)

# SALIDA
print('raiz en: ', raiz)

el resultado es:

raiz en:  0.3649852264049102
>>> funcionx(0.3649852264049102)
1.4099832412739488e-14
>>>