s1Eva_IT2019_T2 Catenaria cable

Ejercicio: 1Eva_IT2019_T2 Catenaria cable

Las fórmulas con las que se requiere trabajar son:

y = \frac{T_A}{w} cosh \Big( \frac{w}{T_A}x \Big) + y_0 - \frac{T_A}{w}

Donde la altura y del cable está en función de la distancia x.

Además se tiene que:

cosh(z) = \frac{e^z+ e^{-z}}{2}

que sustituyendo la segunda en la primera se convierte en:

y = \frac{T_A}{w} \frac{e^{\frac{w}{T_A}x}+ e^{-\frac{w}{T_A}x}}{2} + y_0 - \frac{T_A}{w}

y usando los valores del enunciado w=12, y0=6 , y=15, x=50 se convierte en:

15 = \frac{T_A}{12} \frac{e^{\frac{12}{T_A}50}+ e^{-\frac{12}{T_A}50}}{2} + 6 - \frac{T_A}{12}

simplificando, para usar el método de búsqueda de raices:

\frac{1}{2}\frac{T_A}{12} e^{\frac{12}{T_A}50} + \frac{1}{2}\frac{T_A}{12} e^{-\frac{12}{T_A}50} - \frac{T_A}{12} - 9 = 0

cambiando la variable \frac{12}{T_A}=x

\frac{1}{2x} e^{50x} + \frac{1}{2x} e^{-50x} - \frac{1}{x}-9=0

la función a usar para la búsqueda de raices es:

f(x)=\frac{1}{2x} e^{50x} + \frac{1}{2x} e^{-50x} - \frac{1}{x}-9

Para el método de Newton-Raphson se tiene que:

x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)}

por lo que se determina:

f'(x)= - \frac{1}{2x^2}e^{50x} + \frac{1}{2x}(50) e^{50x} + - \frac{1}{2x^2} e^{-50x} + \frac{1}{2x}(-50)e^{-50x} + \frac{1}{x^2} f'(x)= -\frac{1}{2x^2}[e^{50x}+e^{-50x}] + + \frac{25}{x}[e^{50x}-e^{-50x}] +\frac{1}{x^2} f'(x)= \Big[\frac{25}{x} -\frac{1}{2x^2}\Big]\Big[e^{50x}+e^{-50x}\Big] +\frac{1}{x^2}

Con lo que se puede inicar las iteraciones.

Por no disponer de valor inicial para TA, considere que el cable colgado no debería tener tensión TA=0 N, pues en la forma x=12/TA se crea una indeterminación. Si no dispone de algún criterio para seleccionar el valor de TA puede iniciar un valor positivo, por ejemplo 120 con lo que el valor de x0=12/120=0.1

Iteración 1

f(0.1)=\frac{1}{2(0.1)} e^{50(0.1)} + \frac{1}{2(0.1)} e^{-50(0.1)} - \frac{1}{0.1}-9 =723.0994 f'(0.1)=\Big[\frac{25}{0.1} - \frac{1}{2(0.1)^2}\Big]\Big[e^{50(0.1)}+e^{-50(0.1)}\Big] + +\frac{1}{(0.1)^2} = 29780.61043 x_{1} = 0.1 -\frac{723.0994}{29780.61043} = 0.07571

error = | x1 – x0| = | 0.07571 – 0.1| = 0.02428

Iteración 2

f(0.07571)=\frac{1}{2(0.07571)} e^{50(0.07571)}+ + \frac{1}{2(0.07571)} e^{-50(0.07571)} - \frac{1}{0.07571}-9 = 269.0042 f'(0.07571)= \Big[\frac{25}{0.07571} -\frac{1}{2(0.07571)^2}\Big]. .\Big[e^{50(0.07571)}+e^{-50(0.07571)}\Big] + +\frac{1}{(0.07571)^2} = 10874.0462 x_{2} = 0.07571 -\frac{269.0042}{10874.0462} = 0.05098

error = | x2 – x1| = |0.05098- 0.02428| = 0.02473

Iteración 3

f(0.05098) = 97.6345 f'(0.05098) = 4144.1544 x_{3} = 0.0274

error = | x3 – x2| = |0.05098- 0.0274| = 0.0236

finalmente después de varias iteraciones, la raiz se encuentra en: 0.007124346154337298

que convitiendo

T_A = \frac{12}{x} = \frac{12}{0.0071243461} = 1684.36 N

Revisión de resultados

Usando como base los algoritmos desarrollados en clase:

['xi', 'xnuevo', 'tramo']
[0.1    0.0757 0.0243]
[0.0757 0.051  0.0247]
[0.051  0.0274 0.0236]
[0.0274 0.0111 0.0163]
[0.0111 0.0072 0.0039]
[7.2176e-03 7.1244e-03 9.3199e-05]
[7.1244e-03 7.1243e-03 3.8351e-08]
raiz en:  0.007124346154337298
TA = 12/x =  1684.365096815854

Algoritmos Python usando el procedimiento de:

http://blog.espol.edu.ec/analisisnumerico/2-3-1-newton-raphson-ejemplo01/

# 1Eva_IT2019_T2 Catenaria cable
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
a = 0.001
b = 0.1
muestras = 51

x0 = 0.1
tolera = 0.00001

fx = lambda x: 0.5*(1/x)*np.exp(50*x) + 0.5*(1/x)*np.exp(-50*x)-1/x -9
dfx = lambda x: -0.5*(1/(x**2))*(np.exp(50*x)+np.exp(-50*x)) + (25/x)*(np.exp(50*x)-np.exp(-50*x)) + 1/(x**2)

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

TA = 12/xnuevo

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

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

# Grafica
plt.plot(xp,fp)
plt.xlabel('x=12/TA')
plt.ylabel('f(x)')
plt.axhline(0, color = 'green')
plt.grid()
plt.show()

s1Eva_IT2019_T1 Oxígeno y temperatura en agua

Ejercicio: 1Eva_IT2019_T1 Oxígeno y temperatura en agua

Literal a

Se requiere un polinomio de grado 3 siendo el eje x correspondiente a temperatura. Son necesarios 4 puntos de referencia alrededor de 15 grados, dos a la izquierda y dos a la derecha.

Se observa que los datos en el eje x son equidistantes, h=8, y ordenados en forma ascendente, se cumple con los requisitos para usar diferencias finitas avanzadas. que tiene la forma de:

p_n (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1) + + \frac{\Delta^3 f_0}{3!h^3} (x - x_0)(x - x_1)(x - x_2) + \text{...}

Tabla

xi f[xi] f[x1,x0] f[x2,x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0]
8 11.5 9.9-11.5=
-1.6
-1.5-(-1.6) =
0.1
0.4-0.1=
0.3
16 9.9 8.4-9.9=
-1.5
-1.1-(1.5)=
0.4
24 8.4 7.3-8.4=
-1.1
32 7.3

Con lo que el polinomio buscado es:

p_3 (x) = 11.5 + \frac{-1.6}{8} (x - 8) + + \frac{0.1}{2!8^2} (x - 8)(x - 16) + \frac{0.3}{3!8^3} (x - 8)(x - 16)(x - 24)

Resolviendo y simplificando el polinomio, se puede observar que al aumentar el grado, la constante del término disminuye.

p_3(x)=12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3

para el cálculo del error se puede usar un término adicional del polinomio, añadiendo un punto más a la tabla de diferencia finitas. Se evalúa éste término y se estima el error que dado que el término de grado 3 es del orden de 10-5, el error será menor. (Tarea)

p_3(15)=12.9- 0.15 (15) - 0.00390625 (15)^2 + 0.00009765625 (15)^3

Evaluando el polinomio en temperatura = 15:

p3(15) = 10.1006835937500

literal b

se deriva el polinomio del literal anterior y se evalua en 16:

p'_3(x)=0- 0.15 - 0.00390625 (2) x + 0.00009765625 (3)x^2 p'_3(16)=0- 0.15 - 0.00390625 (2)(16) + 0.00009765625 (3)(16)^2

p’3(16) = -0.20

literal c

El valor de oxígeno usado como referencia es 9, cuyos valores de temperatura se encuentran entre 16 y 24 que se toman como rango inicial de búsqueda [a,b]. Por lo que el polinomio se iguala a 9 y se crea la forma estandarizada del problema:

p_3(x)=9 9 = 12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3 12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3 -9 = 0 f(x) = 3.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3

Para mostrar el procedimiento se realizan solo tres iteraciones,

1ra Iteración
a=16 , b = 24, c = (16+24)/2 = 20
f(a) = 0.9, f(b) = -0.6, f(c) = 0.011
error = |24-16| = 8
como f(c) es positivo, se mueve el extremo f(x) del mismo signo, es decir a

2da Iteración
a=20 , b = 24, c = (20+24)/2 = 22
f(a) = 0.119, f(b) = -0.6, f(c) = -0.251
error = |24-20|= 4
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b

3ra Iteración
a=20 , b = 22, c = (20+22)/2 = 21
f(a) = 0.119, f(b) = -0.251, f(c) = -0.068
error = |22-20| = 2
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b
y así sucesivamente hasta que error< que 10-3

Usando el algoritmo en python se obtendrá la raiz en 20.632 con la tolerancia requerida.


Revisión de resultados

Usando como base los algoritmos desarrollados en clase:

Literal a
[[ 1.   8.  11.5 -1.6  0.1  0.3  0. ]
 [ 2.  16.   9.9 -1.5  0.4  0.   0. ]
 [ 3.  24.   8.4 -1.1  0.   0.   0. ]
 [ 4.  32.   7.3  0.   0.   0.   0. ]]
9.765625e-5*x**3 - 0.00390625*x**2 - 0.15*x + 12.9
p(15) =  10.1006835937500
Literal b
0.00029296875*x**2 - 0.0078125*x - 0.15
dp(16) = -0.200000000000000

literal c
[ i, a, c, b, f(a), f(c), f(b), paso]
1 16.000 20.000 24.000 0.900 0.119 -0.600 8.000 
2 20.000 22.000 24.000 0.119 -0.251 -0.600 4.000 
3 20.000 21.000 22.000 0.119 -0.068 -0.251 2.000 
4 20.000 20.500 21.000 0.119 0.025 -0.068 1.000 
5 20.500 20.750 21.000 0.025 -0.022 -0.068 0.500 
6 20.500 20.625 20.750 0.025 0.001 -0.022 0.250 
7 20.625 20.688 20.750 0.001 -0.010 -0.022 0.125 
8 20.625 20.656 20.688 0.001 -0.004 -0.010 0.062 
9 20.625 20.641 20.656 0.001 -0.002 -0.004 0.031 
10 20.625 20.633 20.641 0.001 -0.000 -0.002 0.016 
11 20.625 20.629 20.633 0.001 0.001 -0.000 0.008 
12 20.629 20.631 20.633 0.001 0.000 -0.000 0.004 
13 20.631 20.632 20.633 0.000 0.000 -0.000 0.002 
14 20.632 20.632 20.633 0.000 0.000 -0.000 0.001 
raiz:  20.63232421875

Algoritmos Python usando la funcion de interpolación y un procedimiento encontrado en:

http://blog.espol.edu.ec/analisisnumerico/5-1-1-diferencias-finitas-avanzadas-polinomio/

http://blog.espol.edu.ec/analisisnumerico/2-1-1-biseccion-ejemplo01/

# 1Eva_IT2019_T1 Oxígeno y temperatura en mar
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

def interpola_dfinitas(xi,fi):
    '''
    Interpolación de diferencias finitas avanzadas
    resultado: polinomio en forma simbólica
    '''
    # Tabla de diferencias finitas
    titulo = ['i','xi','fi']
    n = len(xi)
    # cambia a forma de columnas
    i = np.arange(1,n+1,1)
    i = np.transpose([i])
    xi = np.transpose([xi])
    fi = np.transpose([fi])
    # Añade matriz de diferencias
    dfinita = np.zeros(shape=(n,n),dtype=float)
    tabla = np.concatenate((i,xi,fi,dfinita), axis=1)
    # Sobre matriz de diferencias, por columnas
    [n,m] = np.shape(tabla)
    c = 3
    diagonal = n-1
    while (c<m):
        # Aumenta el título para cada columna
        titulo.append('df'+str(c-2))
        # calcula cada diferencia por fila
        f = 0
        while (f < diagonal):
            tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
            f = f+1
        
        diagonal = diagonal - 1
        c = c+1

    # POLINOMIO con diferencias finitas
    # caso: puntos en eje x equidistantes
    dfinita = tabla[:,3:]
    n = len(dfinita)
    x = sym.Symbol('x')
    h = xi[1,0]-xi[0,0]
    polinomio = fi[0,0]
    for c in range(1,n,1):
        denominador = np.math.factorial(c)*(h**c)
        factor = dfinita[0,c-1]/denominador
        termino=1
        for f  in range(0,c,1):
            termino = termino*(x-xi[f])
        polinomio = polinomio + termino*factor
    # simplifica polinomio, multiplica los (x-xi)
    polinomio = polinomio.expand()
    
    return(tabla, polinomio)

# INGRESO
tm = [0.,8,16,24,32,40]
ox = [14.6,11.5,9.9,8.4,7.3,6.4]

xi = [8,16,24,32]
fi = [11.5,9.9,8.4,7.3]

# PROCEDIMIENTO
x = sym.Symbol('x')
# literal a
tabla, polinomio = interpola_dfinitas(xi,fi)
p15 = polinomio.subs(x,15)
# literal b
deriva = polinomio.diff(x,1)
dp16 = deriva.subs(x,16)

px =  sym.lambdify(x,polinomio)
xk = np.linspace(np.min(xi),np.max(xi))
pk = px(xk)

# SALIDA
print('Literal a')
print(tabla)
print(polinomio)
print('p(15) = ',p15)
print('Literal b')
print(deriva)
print('dp(16) =', dp16)

# gráfica
plt.plot(tm,ox,'ro')
plt.plot(xk,pk)
plt.axhline(9,color="green")
plt.xlabel('temperatura')
plt.ylabel('concentracion de oxigeno')
plt.grid()
plt.show()

# --------literal c ------------

# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera

# se convierte forma de símbolos a numéricos
buscar = polinomio-9
fx = sym.lambdify(x,buscar)

# INGRESO
a = 16
b = 24
tolera = 0.001

# PROCEDIMIENTO
tabla = []
tramo = b-a

fa = fx(a)
fb = fx(b)
i = 1
while (tramo > tolera):
    c = (a+b)/2
    fc = fx(c)
    tabla.append([i,a,c,b,fa,fc,fb,tramo])
    i = i+1
                 
    cambia = np.sign(fa)*np.sign(fc)
    if (cambia<0):
        b = c
        fb = fc
    else:
        a=c
        fa = fc
    tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)

raiz = c

# SALIDA
print('\n literal c')
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), paso]')
# print(tabla)

# Tabla con formato
n=len(tabla)
for i in range(0,n,1):
    unafila = tabla[i]
    formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
    unafila = formato.format(*unafila)
    print(unafila)
    
print('raiz: ',raiz)

s1Eva_IT2010_T1_MN Demanda y producción sin,log

Ejercicio: 1Eva_IT2010_T1_MN Demanda y producción sin,log

Desarrollo Analítico

Para la demanda, el intervalo de existencia es [0,3]

demanda(t) = sin(t)

Para la oferta, el intervalo de existencia inicia en 1, limitado por la demanda [1,3]

oferta(t) = ln(t)

la oferta satisface la demanda cuando ambas son iguales

demanda(t) = oferta(t) sin(t) = ln(t)

por lo que el tiempo t se encuentra con algun método para determinar la raiz de:

sin(t) - ln(t) = 0 f(t) = sin(t) - ln(t)

Observe que las curvas de oferta y demanda se intersectan en el mismo punto en el eje x que la función f(t).

Use un método para encontrar el valor t que satisface la ecuación.


Algoritmo en Python

instrucciones para la gráfica, no para el algoritmo de búsqueda de raiz que es tarea.

# 1Eva_IT2010_T1_MN Demanda y producción sin,log
# Solo para analizar el problema
# Tarea: Añadir algoritmo de buscar raiz.
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# demanda
ad = 0
bd = 3
muestrasd = 31
# oferta
ao = 1
bo = 3
muestras0  = 21

demanda = lambda t: np.sin(t)
oferta = lambda t: np.log(t)
f = lambda t: demanda(t)-oferta(t)

# PROCEDIMIENTO
tid = np.linspace(ad,bd,muestrasd)
demandai = demanda(tid)

tio = np.linspace(ao,bo,muestras0)
ofertai = oferta(tio)

fi = f(tio)

# SALIDA
plt.plot(tid,demandai, label='demanda')
plt.plot(tio,ofertai, label ='oferta')
plt.plot(tio,fi,label='f(t)= demanda-oferta')
plt.axhline(0,color='black')
plt.axvline(2.2185, color = 'magenta')
plt.xlabel('tiempo')
plt.ylabel('unidades')
plt.legend()
plt.grid()
plt.show()

s1Eva_IT2017_T4 Componentes eléctricos

Ejercicio: 1Eva_IT2017_T4 Componentes eléctricos

Desarrollo Analítico

Solo puede usar las cantidades disponibles de material indicadas, por lo que las cantidades desconocidas de produccion por componente se convieren en las incógnitas x0, x1, x2. Se usa el índice cero por compatibilidad con las instrucciones python.

Material 1 Material 2 Material 3
Componente 1 5 x0 9 x0 3 x0
Componente 2 7 x1 7 x1 16 x1
Componente 3 9 x2 3 x2 4 x2
Total 945 987 1049

Se plantean las fórmulas a resolver:

5 x0 +  7 x1 + 9 x2 = 945
9 x0 +  7 x1 + 3 x2 = 987
3 x0 + 16 x1 + 4 x2 = 1049

Se reescriben en la forma matricial AX=B

\begin{bmatrix} 5 & 7 & 9\\ 9 & 7 & 3 \\ 3& 16 & 4\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 945 \\ 987 \\ 1049 \end{bmatrix}

Se reordena, pivoteando por filas para tener la matriz diagonalmente dominante:

\begin{bmatrix} 9 & 7 & 3\\ 3 & 16 & 4 \\ 5& 7 & 9\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 987 \\ 1049 \\ 945 \end{bmatrix}

Se determina el número de condición de la matriz. Por facilidad con python:

numero de condicion: 4.396316324708121

Obtenido con:

# 1Eva_IT2017_T4 Componentes eléctricos
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
A = np.array([[9., 7.,3.],
              [3.,16.,4.],
              [5., 7.,9.]])

B = np.array([987.,1049.,945.])
# PROCEDIMIENTO

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

# SALIDA
print('numero de condicion:', ncond)

Recordando que la matriz debe ser tipo real, se añade un punto a los dígitos.

El número de condición es cercano a 1, por lo que el sistema si debería converger a la solución.

Desarrollo con Python

La forma AX=B permite usar los algoritmos desarrollados, obteniendo la solución. Se verifica el resultado al realizar la multiplicación de A con el vector respuesta, debe ser el vector B con un error menor al tolerado.

respuesta con Jacobi
[[62.99996585]
 [44.99998037]
 [34.9999594 ]]
verificando:
[[ 986.99943346]
 [1048.99942111]
 [ 944.99932646]]

Si interpreta el resultado, se debe obtener solo la parte entera [63,45,35] pues las unidades producidas son números enteros.

instrucciones:

# 1Eva_IT2017_T4 Componentes eléctricos
import numpy as np

def jacobi(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    xnuevo = np.copy(X)

    itera = 0
    print(itera, X)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        print(itera, xnuevo)
        X = np.copy(xnuevo)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
    return(X)


# INGRESO
A = np.array([[9., 7.,3.],
              [3.,16.,4.],
              [5., 7.,9.]])

B = np.array([987.,1049.,945.])
tolera = 1e-4

X = [0.,0.,0.]

# PROCEDIMIENTO

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

respuesta = jacobi(A,B,tolera,X)

verifica = np.dot(A,respuesta)

# SALIDA
print('numero de condicion:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('verificando:')
print(verifica)

al ejecutar el algoritmo se determina que se requieren 83 iteraciones para cumplir con con el valor de error tolerado.

s1Eva_IIT2018_T4 Tasa de interés en hipoteca

Ejercicio: 1Eva_IIT2018_T4 Tasa de interés en hipoteca

literal a

Siguiendo el desarrollo analítico tradicional, para adecuar la ecuación para los algoritmo de búsquda de raíces de ecuaciones,  se reemplazan los valores en la fórmula.

P = A\Big(\frac{1-(1+i)^{-n}}{i} \Big) 70000 = 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)

Como ambos lados de la ecuación deben ser iguales, si se restan ambos se obtiene una ecuación que tiene como resultado cero, que es la forma ideal para usar en el algoritmo que representa f(x) o en este caso f(i)

70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) = 0

Para evitar inconvenientes con la división para cero en caso que i tome el valor de cero, dado se multiplica toda la ecuación por i:

i \Big[70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) \Big]= i (0) 70000 i - 1200 (1-(1+i)^{-300}) = 0

La ecuación es la utilizada en el algoritmo de búsqueda de raíces pueden ser:

fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) fx(i) = 70000i - 1200(1-(1+i)^{-300})

literal b

El intervalo de existencia correspondería a la tasa de interés mínimo y el interés máximo.

[izquierda, derecha] = [a,b]

Para el intervalo se deben tomar en cuenta algunas consideraciones descritas a continuación:

izquierda:

En el extremo izquierdo, las tasas no son negativas, lo que se interpreta en que un banco paga por que le presten dinero.

Tampoco tiene mucho sentido el valor cero, que son prestamos sin intereses. A menos que sean sus padres quienes le dan el dinero.

Un valor inicial para el interés puede ser por ejemplo 1% ó 0.01, siempre que se cumpla que existe cambio de signo en la función a usar.

derecha:

En el extremo derecho, si se propone por ejemplo i con 100%, o 1.00, no tendría mucho sentido un préstamo con intereses al 100% anual, que resulta en el doble del valor inicial en tan solo un periodo o año.

La tasa de interés de consumo que son de las más alto valor, se encuentran reguladas. En Ecuador es un valor alrededor del 16% anuales o 0.16.

Considerando las observaciones iniciales del problema, se propone empezar el análisis para la búsqueda de la raíz en el intervalo en un rango más amplio:

[ 0.01, 0.50]

Ser realiza la comprobación que existe cambio de signo en los extremos del intervalo.

fx(0.001) =- 43935.86

fx(0.50) = 67600.0

Para el ejercicio se hace notar que la es tasa nominal anual, pero los pagos son mensuales. Por lo que se debe unificar las tasas de interes a mensuales. Una aproximación es usar las tasas anuales divididas para los 12 meses del año.

Tolerancia al error

La tolerancia se considera en éste ejercicio como el valor de diferencias  (tramo) entre iteraciones con precisión satisfactoria.

Por ejemplo si no negociaremos más con el banco por variaciones de tasas del 0.1% , entonces la tolerancia será de 0.001.

Las publicaciones de tasas en el mercado incluyen dos decimales, por lo que para el ejercicio aumentamos la precisión a : 0.0001

tolera = 1×10-4


Literal c


Se presentan dos formas se solución para el litera c:

– c.1 la requerida en el enunciado con Newton-Raphson

– c.2 una alterna con el método de la Bisección.


c.1. Desarrollo del ejercicio con el método del enunciado Newton-Raphson


Para el método de Newton-Raphson se tiene que:

x_{i+1} = x_{i} - \frac{f(x_0i)}{f'(x_i)}

Se requiere la derivada de la función planteada en el literal a:

fx(i) = 70000i - 1200(1-(1+i)^{-300}) f'x(i) = 70000 + 1200(300)(1+i)^{-301})

tomando como valor inicial xi = 0.16/12 ≈ 0.013

Se realizan las iteraciones suponiendo que tolera = 1×10-4

iteración 1

fx(0.013) = 70000(0.013) - 1200(1-(1+0.013)^{-300})

 = -265.0914

f'x(0.013) = 70000 + 1200(300)(1+0.013)^{-301})

= 62623.3454

x_{2} = 0.013 - \frac{-265.0914}{63318.8456} = 0.01723

error = |0.013 – 0.01723| = 0.004229

iteración 2

fx(0.01723) = 70000i - 1200(1-(1+0.0.01723)^{-300})

= 13.2356

f'x(0.01723) = 70000 + 1200(300)(1+0.01723)^{-301}

= 67895.5656

x_{3} = 0.01723 - \frac{13.2356}{67895.5656} = 0.01703

error = |0.01723 – 0.01703| = 0.0001999

cuyo valr de error está casi dentro del valor de tolerancia,

que permite tomar el último valor como respuesta de tasa mensual

raiz = tasa mensual = 0.01703

Convirtiendo a la tasa tasa anual que es la publicada por las instituciones financieras se tiene que:

tasa anual nominal =  0.01703*12 = 0.2043

Teniendo como resultado una tasa anual de 20.43%


Algoritmo en Python

# 1ra Evaluación II Término 2018
# Tema 4. Tasa de interes para hipoteca
import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx, dfx, xi, tolera):
    '''
    funciónx y fxderiva son de forma numérica
    xi es el punto inicial de búsqueda
    '''
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = xi - fx(xi)/dfx(xi)
        tramo = abs(xnuevo-xi)
        print(xi,xnuevo, tramo)
        xi = xnuevo
    return(xi)

# INGRESO
P = 70000.00
A = 1200.00
n = 25*12
fx  = lambda i: P*i - A*(1-(1+i)**(-n))
dfx = lambda i: P + A*(-n)*(i+1)**(-n-1)
a = 0.01/12
b = 0.25/12
muestras = 51
tolera = 0.0001

# PROCEDIMIENTO
tasa = np.linspace(a,b,muestras)
fi   = fx(tasa)

raiz   = newton_raphson(fx, dfx, b, tolera)
tanual = 12*raiz

# SALIDA
print('raiz encontrada en: ', raiz)
print('tasa anual: ',tanual)

# Gráfica
plt.plot(tasa*12,fi)
plt.title('tasa anual de interes para Hipoteca')
plt.xlabel('tasa')
plt.ylabel('fx(tasa)')
plt.axhline(0, color='green')
plt.show()

c.2. Desarrollo con el método de la Bisección


Desarrollo Analítico con Bisección

Como parte del desarrollo del ejercicio se presenta las iteraciones para el algoritmo, tradicionalmente realizadas con una calculadora.

fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)

iteración 1

a = 0.01, b = 0.5 c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 fx(0.01) = 70000 - 1200\Big(\frac{1-(1+(0.01))^{-300}}{0.01} \Big) = -43935.86 fx(0.255) = 70000 - 1200\Big(\frac{1-(1+(0.255))^{-300}}{0.255} \Big) = 65294.11 fx(0.5) = 70000 - 1200\Big(\frac{1-(1+(0.5))^{-300}}{0.5} \Big) = 67600.0 tramo = 0.5-0.01 =0.49

cambio de signo a la izquierda

a = 0.01, b=0.255

iteración 2

a = 0.01, b = 0.225 c = \frac{a+b}{2} = \frac{0.01+0.225}{2} = 0.1325 fx(0.01) = -43935.86 fx(0.1325) = 70000 - 1200\Big(\frac{1-(1+(0.1325))^{-300}}{0.1325} \Big) = 60943.39 fx(0.225) = 65294.11 tramo = 0.225-0.01 =0.215

cambio de signo a la izquierda

a = 0.01, b=0.1325

iteración 3

a = 0.01, b = 0.1325 c = \frac{a+b}{2} = \frac{0.01+0.1325}{2} = 0.07125 fx(0.01) = -43935.86 fx(0.07125) = 70000 - 1200\Big(\frac{1-(1+(0.07125))^{-300}}{0.07125} \Big) = 53157.89 fx(0.1325) = 60943.39 tramo = 0.1325-0.01 =0.1225

cambio de signo a la izquierda

a = 0.01, b=0.07125

y se continuaría con las iteraciones, hasta cumplir que tramo<=tolera

Tabla de datos obtenidos

tabla para Bisección
i a c b f(a) f(c) f(b) tramo
1 0.01 0.255 0.5 -43935.86 65294.11 67600.0 0.49
2 0.01 0.1325 0.255 -43935.86 60943.39 65294.11 0.215
3 0.01 0.07125 0.1325 -43935.86 53157.89 60943.39 0.1225

hasta lo calculado la raiz se encontraría en el intervalo [0.01,0.07125] con error esitmado de 0.06125, aún por mejorar con más iteraciones.

Algoritmo en Python para Bisección

  • El algoritmo bisección usa las variables a y b, por lo que los limites en el intervalo usados son [La,Lb]
  • para el problema la variable ‘i’ se usa en el eje x.
  • La selección de cambio de rango [a,b] se hace usando solo el signo del valor.
  • El algoritmo presentado es tal como se explica en la parte conceptual

Se deja como tarea convertir el algoritmo a funcion def-return de Python.

# 1Eva_IIT2018_T4 Tasa de interés en hipoteca
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
P = 70000.00
A = 1200.00
n = 25*12
fi = lambda i: P - A*(1-((1+i)**-n))/i

# Intervalo de observación
# e inicio de Bisección
La = 0.01
Lb = 0.50

tolera = 0.0001 #grafica

muestras = 21

# PROCEDIMIENTO

# Método de Bisección
a = La
b = Lb
c = (a+b)/2
tramo = np.abs(b-a)
while (tramo>tolera):
    fa = fi(a)
    fb = fi(b)
    fc = fi(c)
    cambio = np.sign(fc)*np.sign(fa)
    if (cambio>0):
        a = c
        b = b
    else:   
        b = c
        a = a
    c = (a+b)/2
    tramo = np.abs(b-a)

# Para la gráfica
tasa = np.linspace(La,Lb,muestras)
fr = fi(tasa)

# SALIDA
print('a, f(a):', a,fa)
print('c, f(c):', c,fc)
print('b, f(b):', b,fb)
print('la raiz esta entre: \n',a,b)
print('con un error de: ', tramo)
print('raiz es tasa buscada: ', c)
print('tasas anual buscada: ',c*12)

# Gráfica
plt.plot(tasa,fr)
plt.axhline(0, color='green')
plt.title('tasa de interes mensual')
plt.show()

la ejecución del algoritmo da como resultado

>>> 
 RESTART: D:/MATG1052Ejemplos/HipotecaInteres.py 
a, f(a): 0.016998291015625 -385.52828922150366
c, f(c): 0.0170281982421875 -145.85350695741363
b, f(b): 0.01705810546875 92.28034212642524
la raiz esta entre: 
 0.016998291015625 0.01705810546875
con un error de:  5.981445312500111e-05
raiz es tasa buscada:  0.0170281982421875
tasas anual buscada:  0.20433837890625

y la gráfica obtenida es:

s1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

Ejercicio: 1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

El ejercicio tiene dos partes: la interpolación y el integral.

Literal a

No se especifica el método a seguir, por lo que se puede seleccionar el de mayor preferencia.

Usando Lagrange, con los puntos primero, medio y último:

p_2(t) = 0\frac{(t-4)(t-8)}{(0-4)(0-8)} + + 27.77\frac{(t-0)(t-8)}{(4-0)(4-8)} + + 41.10\frac{(t-0)(t-4)}{(8-0)(8-4)} p_2(t) = 0 + 27.77\frac{t(t-8)}{-16}) + + 41.10\frac{t(t-4)}{32} p_2(t) = -1.73(t^2-8t) + 1.28(t^2-4t) p_2(t) = -0.45 t^2 + 8.72t


Literal b

El tema de integración para primera evaluación se realiza de forma analítica.

Como introducción a la Unidad 7 de la 2da Evaluación, se usa el método del trapecio para el polinomio encontrado en el literal anterior:

v(t) =  -0.45125*t**2 + 8.7475*t
distancia recorrida:  202.8912640000001

Realice la integración analítica del polinomio encontrado y determine el error entre los métodos.


Desarrollo con Python

las instrucciones en Python para el ejercicio son:

# 1ra Evaluación II Término 2018
# Tema 1. Interpolar velocidad del paracaidista

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# Literal a)
def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.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*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

# INGRESO
t = [0.0, 2, 4, 6, 8]
v = [0.0, 16.40, 27.77, 35.64, 41.10]

xi = [t[0],t[2],t[4]]
yi = [v[0],v[2],v[4]]

muestras = 51

# PROCEDIMIENTO
polinomio = interpola_lagrange(xi,yi)
velocidad = polinomio.subs('x','t')

# Para graficar
vt = sym.lambdify('t',velocidad)
a = t[0]
b = t[-1]
ti = np.linspace(a, b, muestras)
vi = vt(ti)

# SALIDA
print('v(t) = ', velocidad)
# Grafica
plt.plot(t,v,'ro')
plt.plot(ti,vi)
plt.title('Interpolar velocidad de paracidista')
plt.xlabel('t')
plt.ylabel('v')
plt.show()


# Literal b
def integratrapecio(funcionx,a,b,tramos):
    h = (b-a)/tramos
    x = a
    suma = funcionx(x)
    for i in range(0,tramos-1,1):
        x = x+h
        suma = suma + 2*funcionx(x)
    suma = suma + funcionx(b)
    area = h*(suma/2)
    return(area)

# INGRESO
# El ingreso es el polinomio en forma lambda
# se mantienen las muestras
tramos = muestras-1
# PROCEDIMIENTO
distancia = integratrapecio(vt,a,b,tramos)

# SALIDA
print('distancia recorrida: ', distancia)

s1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

Ejercicio: 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

El tema es semejante al tema 1, cambiando el método de interpolación.


Literal a

Se usan los puntos de las posiciones 0, 3 y 5.
en la fórmula:

p_2(x) = b_0 + b_1x + b_2 x^2

en la fórmula:

punto x[0] = 1, y[0]= 1.84

1.84 = b_0 + b_1(1) + b_2 (1)^2 1.84 = b_0 + b_1 + b_2

punto x[3] = 1.5, y[3]= 2.28

2.28 = b_0 + b_1(1.5) + b_2 (1.5)^2 2.28 = b_0 + 1.5 b_1 + 2.25 b_2

punto x[5] = 2.1, y[5]= 3.28

3.28= b_0 + b_1(2.1) + b_2 (2.1)^2 3.28= b_0 + 2.1 b_1 + 4.41 b_2

se obtiene el sistema de ecuaciones:

b_0 + b_1 + b_2 = 1.84 b_0 + 1.5 b_1 + 2.25 b_2 = 2.28 b_0 + 2.1 b_1 + 4.41 b_2 = 3.28

Con lo que se plantea la forma Ax=B:

A = \begin{bmatrix} 1 & 1 & 1\\ 1 & 1.5 & 2.25 \\1 & 2.1 & 4.41 \end{bmatrix} B = \begin{bmatrix} 1.84\\ 2.28 \\ 3.28 \end{bmatrix}

y se obtiene el resultado de la interpolación.


Literal b

Se requiere calcular una norma de suma de filas. es suficiente para demostrar el conocimiento del concepto el usar A.

Se adjunta el cálculo del número de condición y la solución al sistema de ecuaciones:

suma de columnas:  [3.   4.75 7.51]
norma A:  7.51
numero de condicion:  97.03737354737129
solucion: 
[ 2.03272727 -0.90787879  0.71515152]

El comentario importante corresponde al número de condición, que es un número muy alto para usar un método iterativo, por lo que la solución debe ser un método directo.
Se puede estimar será un número mucho mayor que 1, pues la matriz no es diagonal dominante.


Instrucciones en Python

# 1ra Evaluación II Término 2018
# Tema 3. Interpolar con sistema de ecuaciones

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# forma matricial para interpolar
A = np.array([[1, 1. , 1.  ],
              [1, 1.5, 2.25],
              [1, 2.1, 4.41]])

B = np.array([1.84, 2.28, 3.28])

# literal b
sumacolumnas = np.sum(A, axis =1)
norma = np.max(sumacolumnas)
print('suma de columnas: ', sumacolumnas)
print('norma A: ', norma)

numerocondicion = np.linalg.cond(A)
print('numero de condicion: ', numerocondicion)

solucion = np.linalg.solve(A,B)
print('solucion: ')
print(solucion)

s1Eva_IIT2018_T2 Distancia mínima a un punto

Ejercicio: 1Eva_IIT2018_T2 Distancia mínima a un punto

Literal a

Se requiere analizar la distancias entre una trayectoria y el punto = [1,1]

Al analizar las distancias de ex y el punto [1,1] se trazan lineas paralelas a los ejes desde el punto [1,1], por lo que se determina que el rango de x = [a,b] para distancias se encuentra en:

a > 0, a = 0.1
b < 1, b = 0.7

El ejercicio usa la fórmula de distancia entre dos puntos:

d = \sqrt{(x_2-x_1)^2+(y_2- y_1)^2}

en los cuales:

[x1,y1] = [1,1]
[x2,y2] = [x, ex]

que al sustituir en la fórmula se convierte en:

d = \sqrt{(x-1)^2+(e^x- 1)^2}

que es lo requerido en el literal a


Literal b

Para encontrar el punto más cercano, se debe encontrar el mínimo de la distancia, se podría derivar la función y encontrar la raiz en cero.

Considere simplificar la función a un polinomio, donde tiene dos opciones:

b.1 Polinomio de Taylor, que también requiere derivadas (descartado)

b.2 evaluar la función en varios puntos, interpolar y obtener un polinomio.

Al reutilizar el algoritmo del tema 1 se obtiene lo planteado en b.2, usando un polinomio de grado 3, con muestras de 4 puntos equidistantes en el eje x.

polinomio
0.867192074184622*x**3 + 1.22015957396232*x**2 - 1.21861610672236*x + 1.01491694350023
derivada polinomio:
2.60157622255387*x**2 + 2.44031914792464*x - 1.21861610672236

Al aplicar un método para encontrar raíces se tiene que:

en distancia mínima en x= 0.3644244280922699
con y = 1.4396851273165785
distancia =  0.7728384816953889


Instrucciones en Python

# 1ra Evaluación II Término 2018
# Tema 2. Distancia mínima a un punto
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# PRESENTA PROBLEMA
# INGRESO
punto = [1,1]
trayectoria = lambda x: np.exp(x)

a = 0
b = 1
muestras = 51

# PROCEDIMIENTO
xif = np.linspace(a,b,muestras)
trayecto = trayectoria(xif)

# SALIDA
plt.plot(xif,trayecto, label = 'trayectoria = e^x')
plt.plot(punto[0],punto[1],'ro')
plt.axhline(punto[0], color='grey')
plt.axvline(punto[1], color='grey')
plt.title('distancia a un punto')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# ------------------------
# PARA ANALIZAR DISTANCIAS

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.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*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

def posicionfalsa(fx,a,b,tolera):
    fa = fx(a)
    fb = fx(b)
    c = b - fb*(a-b)/(fa-fb)
    
    tramo = abs(c-a)
    while (tramo > tolera):
        fc = fx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if (cambia > 0):
            tramo = abs(c-a)
            a=c
            fa=fc
        else:
            tramo = abs(c-b)
            b=c
            fb=fc
        c = b - fb*(a-b)/(fa-fb)
        print(tramo)
    respuesta = c
    
    # Valida respuesta
    fa = fx(a)
    fb = fx(b)
    cambio = np.sign(fa)*np.sign(fb)
    if (cambio>0):
        respuesta = np.nan
        
    return(respuesta)

# INGRESO
# Trayectoria y punto tomados de sección PRESENTA PROBLEMA
y = lambda x: np.sqrt((x-punto[0])**2+(trayectoria(x)-punto[1])**2)

a = 0.1
b = 0.7
muestras = 51

tolera = 0.0001
muestrap = 4 # Para el polinomio

# PROCEDIMIENTO
yif = y(xif)
xip = np.linspace(a,b,muestrap)
yip = y(xip)

x = sym.Symbol('x')

polinomio = interpola_lagrange(xip,yip)

px = sym.lambdify('x',polinomio)
pxi = px(xif)

dpx = polinomio.diff('x',1)
dpxn = sym.lambdify('x',dpx)
fxi = dpxn(xif)

raiz = posicionfalsa(dpxn,a,b,tolera)

# SALIDA
print('polinomio')
print(polinomio)
print('derivada polinomio:')
print(dpx)
print('en distancia mínima en x=',raiz)
print('con y =', trayectoria(raiz))
print('distancia = ', y(raiz))

# GRAFICA
plt.plot(xif,yif, label = 'distancia')
plt.plot(xif,pxi, label = 'p(x)')
plt.plot(xif,fxi, label = 'dpx(x)')
plt.axhline(0, color = 'grey')
plt.axvline(raiz,  color = 'grey')
plt.legend()
plt.title('Distancia mínima de (1,1) a e^x')
plt.xlabel('x')
plt.ylabel('distancia')
plt.show()

s1Eva_IIT2007_T2 Aplicar Gauss-Seidel

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

s1Eva_IT2011_T1 Encontrar α en integral

Ejercicio: 1Eva_IT2011_T1 Encontrar α en integral

Desarrollo Analítico

Se iguala la ecuación al valor buscado = 10, y se resuelve

\int_{\alpha}^{2\alpha} x e^{x}dx = 10

siendo: μ = x , δv = ex, δu = δx , v = ex

\int u dv = uv - \int v \delta u xe^x \Big|_{\alpha}^{2 \alpha} - \int_{\alpha}^{2\alpha} e^{x}dx - 10 = 0 2\alpha e^{2 \alpha} -\alpha e^{\alpha} - (e^{2\alpha} - e^{\alpha}) - 10 = 0 (2\alpha-1)e^{2 \alpha}+ (1-\alpha) e^{\alpha} - 10 = 0

la función a usar en el método es

f(\alpha) = (2\alpha-1)e^{2 \alpha}+ (1-\alpha)e^{\alpha} -10

Se obtiene la derivada para el método de Newton Raphson

f'(\alpha) = 2e^{2 \alpha} + 2(2\alpha-1)e^{2 \alpha} - e^{\alpha} + (1-\alpha) e^{\alpha} f'(\alpha) = (2 + 2(2\alpha-1))e^{2 \alpha} +(-1 + (1-\alpha)) e^{\alpha} f'(\alpha) = 4\alpha e^{2 \alpha} -\alpha e^{\alpha}

la fórmula para el método de Newton-Raphson

\alpha_{i+1} = \alpha_i - \frac{f(\alpha)}{f'(\alpha)}

se prueba con α0 = 1, se evita el valor de cero por la indeterminación que se da por f'(0) = 0

iteración 1:

f(1) = (2(1)-1)e^{2(1)}+ (1-(1))e^{(1)} -10 f'(1) = 4(1) e^{2 (1)} -(1) e^{(1)} \alpha_{1} = 1- \frac{f(1)}{f'(1)}

\alpha_{1} = 1.0973
error = 0.0973

iteración 2:
f(2) = (2(1.0973)-1)e^{2(1.0973)}+ (1-(1.0973))e^{(1.0973)} -10

f'(2) = 4(1.0973) e^{2 (1.0973)} -(1.0973) e^{(1.0973)} \alpha_{2} = 1.0973 - \frac{f(1.0973)}{f'(1.0973)}

\alpha_{2} = 1.0853
error = 0.011941

iteración 3:
f(3) = (2(1.0853)-1)e^{2(1.0853)}+ (1-(1.0853))e^{(1.0853)} -10

f'(3) = 4(1.0853) e^{2 (1.0853)} -(1.0853) e^{(1.0853)} \alpha_{3} = 1.0853- \frac{f(1.0853)}{f'(1.0853)}

\alpha_{3} = 1.0851
error = 0.00021951

[  xi,          xnuevo,      f(xi),      f'(xi),       tramo ]
[  1.           1.0973      -2.6109      26.8379       0.0973]
[  1.0973e+00   1.0853e+00   4.3118e-01   3.6110e+01   1.1941e-02]
[  1.0853e+00   1.0851e+00   7.6468e-03   3.4836e+01   2.1951e-04]
[  1.0851e+00   1.0851e+00   2.5287e-06   3.4813e+01   7.2637e-08]
raiz:  1.08512526549

se obtiene el valor de la raíz con 4 iteraciones, con error de aproximación de 7.2637e-08

Desarrollo con Python