s1Eva_IIT2018_T4 Tasa de interés en hipoteca

literal a

Para adecuar la ecuación, se reemplazan los valores en la fórmula y se iguala a cero.

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

Para evitar inconvenientes con la división para cero se multiplica toda la ecuación por i:

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

literal b

El intervalo de existencia corresponderia a sin interesess e interés máximo de 100%, pero no tendría mucho sentido un preéstamo de 100% anual. por lo que empezar el análisis con 25% puede servir par empezar la búsqueda de la raíz.

Un detalle importante es que las cuotas son mensuales, por lo que las tasas son mensuales, por lo que se busca tasas en el intervalo [ 0/12, 0.25/12]

Se evalua el signo fx(0) = 0, que para que el valor resulte negativo se incrementa al 5% anual 0.05/12 mensual. Encontrado que fx(0.05/12) es negativo, se revisa el valor fx(0.25/12) que resulta positivo.

Encontrado el cambio de signo, se usa el intervalo

[ 0.05/12, 0.25/12]


Literal c

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.25/12 = 0.020833

Se realizan las iteraciones suponiendo que tolera = 1e-4

iteración 1

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

 = 260.80

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

= 69274.06

x_{i+1} = 0.25/12 - \frac{260.80}{69274.06} = 0.017068

error = |0.25/12 – 0.017068| = 0.003764

iteración 2

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

=2.2805

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

= 67792.56

x_{i+1} = 0.25/12 - \frac{2.2805}{67792.56} = 0.0170347

error = |0.017068 – 0.017034| = 3.36e-05

que permite tomar el último valor y convertirlo a tasa anual = 12*0.017034 = 0.204418

Teniendo como resultado una tasa anual de 20.44%


Instrucciones 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(funcionx, fxderiva, 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 - funcionx(xi)/fxderiva(xi)
        tramo = abs(xnuevo-xi)
        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.04/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)

plt.plot(tasa,fi)
plt.title('tasa de interes para Hipoteca')
plt.xlabel('tasa')
plt.ylabel('fx(tasa)')
plt.axhline(0, color='grey')
plt.show()

s1Eva_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

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

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

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:

s1Eva_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

s1Eva_IT2009_T1 Demanda de producto

Desarrollo analítico

– igualar la ecuación al valor buscado, 80

200 t e^{-0.75t} = 80

– forma estándar de la ecuación para el método f(x) = 0:

f(t) = 200 t e^{-0.75t} - 80

– derivada de la ecuación

f'(t) = 200 e^{-0.75t} + 200 t (-0.75) e^{-0.75t} f'(t) = 200 e^{-0.75t}(1-0.75t)

– fórmula del método de Newton-Raphson

t_{i+1} = t_i - \frac{f(t)}{f'(t)}

– Punto inicial. Como la variable es t, tiempo, el rango de análisis es t>0
El valor inicial de búsqueda se selecciona t0 = 1

iteración 1:

t_{1} = 1 - \frac{200 (1) e^{-0.75(1)} - 80}{200 e^{-0.75(1)}(1-0.75(1))}

error = 0.6128

iteración 2:

t_{2} = 0.3872 - \frac{200 (0.3872) e^{-0.75(0.3872)} - 80}{200 e^{-0.75(0.3872)}(1-0.75(0.3872))}

error = 0.208

iteracion 3:

t_{3} = 0.5952 - \frac{200 (0.5952) e^{-0.75(0.5952)} - 80}{200 e^{-0.75(0.5952)}(1-0.75(0.5952))}

error = 5.3972e-02

tabla de iteraciones

[  xi,       xnuevo,   f(xi),   f'(xi),    tramo]
[  1.        0.3872   14.4733   23.6183   0.6128]
[  0.3872    0.5952  -22.0776  106.1511   0.208 ]
[  5.9518e-01   6.4916e-01  -3.8242e+00   7.0855e+01   5.3972e-02]
[  6.4916e-01   6.5252e-01  -2.1246e-01   6.3069e+01   3.3686e-03]
[  6.5252e-01   6.5254e-01  -7.9031e-04   6.2600e+01   1.2625e-05]
raiz:  0.65253630273

Se obtiene el valor de la raíz con 5 iteraciones, y un error de 1.2625e-05

Desarrollo con Python:

Tarea: desarrollar el literal b

s1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

1. Desarrollo analítico

1.1 Solución iterativa usando el médodo de Jacobi

El sistema de ecuaciones

\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}

cambia a su forma matricial AX=B

\begin{bmatrix} -2 & 5 & 9 \\ 7 & 1 & 1 \\ -3 & 7 &-1 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 6 \\ -26 \end{bmatrix}

para usarla en el algoritmo se intercambian filas, pivoteo, buscando hacerla diagonal dominante:

\begin{bmatrix} 7 & 1 & 1 \\ -3 & 7 &-1\\ -2 & 5 & 9 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 6 \\ -26 \\1 \end{bmatrix}

las ecuaciones para el algoritmo se obtienen despejando una variable diferente en cada ecuación.

\begin{cases}x=(6-y-z)/(7)\\y=(-26+3x+z)/(7)\\z=(1+2x-5y)/(9)\end{cases}

Dado el vector inicial X(0) = [0, 0, 0], y una tolerancia de 0.0001, e desarrollan al menos 3 iteraciones:

Nota. El super índice entre paréntesis denota el número de la iteración. El subíndice denota el elemento del vector. El error se verifica para éste ejercicio como el mayor valor de la diferencia entre iteraciones consecutivas. Analogía al video del acople entre las aeronaves. Si una coordenada aún no es la correcta …. menor a lo tolerado, pues NO hay acople…

Iteración 1

\begin{cases}x=(6-(0)-(0))/(7)\\y=(-26+3(0)+(0))/(7)\\z=(1+2(0)-5(0))/(9)\end{cases}
X(1) = [6/7, -26/7, 1/9] 
    = [0.8571, -3,7142, 0.1111]
diferencia(1) = |X(1) - X(0)|
             = |[0.8571,-3,7142,0.1111] - [0,0,0]|
             = |[0.8571,-3,7142,0.1111]|
errado(1) = maximo|[0.8571,-3,7142,0.1111]|
         = 26/7 ≅ 3.71, 
error es más alto que lo tolerado

Iteración 2

\begin{cases}x=(6-(-26/7)-(1/9))/(7)\\y=(-26+3(6/7)+(1/9))/(7)\\z=(1+2(6/7)-5(-26/7))/(9)\end{cases}
X(2) = [1.3718, -3.3310 ,  2,3650]
diferencia(2) = |X(2) - X(1)|
             = |[1.3718, -3.3310 , 2,3650] - [0.8571, -3,7142, 0.1111]| 
             = |[0.5146, 0.3832, 2.2539]|
errado(2) ≅ 2.2538, 
el error disminuye, pero es más alto que lo tolerado

Iteración 3

\begin{cases}x=(6-(-3.3310)-( 2,3650))/(7)\\y=(-26+3x+z)/(7)\\z=(1+2x-5y)/(9)\end{cases}
X(3) = [0.9951, -2.7884, 2.2665]
diferencia(3) = |[0.9951, -2.7884, 2.2665] - [1.3718, -3.3310 ,  2,3650]|
             = |[-0.3767, 0.5425, -0.0985]|
errado(3) ≅ 0.5425, 
el error disminuye, pero es más alto que lo tolerado

Observación: Si el error disminuye en cada iteración, se puede intuir que se converge a la solución. Se puede continuar con la 4ta iteración…


2. Solución numérica usando Python

Usando el algoritmo desarrollado en clase se obtiene la respuesta con 13 iteraciones:

Xi, errado
[ 0.85714286 -3.71428571  0.11111111] 3.71428571429
[ 1.37188209 -3.33106576  2.36507937] 2.25396825397
[ 0.99514091 -2.78846777  2.26656589] 0.542597991578
[ 0.93170027 -2.96400162  1.8814023 ] 0.385163589245
[ 1.0117999  -3.04621384  1.96482318] 0.0834208883016
[ 1.01162724 -2.99996816  2.02829656] 0.0634733731283
[ 0.99595309 -2.99097453  2.00256614] 0.0257304176216
[ 0.99834406 -3.0013678   1.99408654] 0.0103932672897
[ 1.00104018 -3.00155447  2.0003919 ] 0.00630536414945
[ 1.00016608 -2.99949822  2.00109475] 0.00205624814582
[ 0.99977193 -2.99977243  1.99975814] 0.00133660433338
[ 1.00000204 -3.0001323   1.99982289] 0.000359867509479
[ 1.0000442  -3.00002443  2.00007395] 0.000251063280064
[ 0.99999292 -2.99997049  2.00002339] 5.39347670632e-05
iteraciones:  13

La gráfica muestra las coordenadas de las aproximaciones, observe el espiral que se va cerrando desde el punto inicial en verde al punto final en rojo.

Tarea: convierta el algoritmo a una función, así es más sencillo usar en otros ejercicios.
Se debe controlar el número de iteraciones para verificar la convergencia con iteramax.
NO es necesario almacenar los valores de los puntos, solo el último valor y el contador itera permite determinar si el sistema converge.

# 1ra Evaluación II Término 2011
# Tema 2. Sistema ecuaciones Jacobi
import numpy as np

# INGRESO
A = np.array([[ 7.0, 1, 1],
              [-3.0, 7,-1],
              [-2.0, 5, 9]])

B = np.array([6.0, -26.0, 1.0])

X = np.array([0.0, 0.0,  0.0],dtype=float)
tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
Xi1 = np.zeros(n, dtype=float)
errado = 2*tolera
itera = 0
while not(errado<=tolera or itera>=iteramax):
    i = 0
    while not(i>=n):
        j = 0
        nuevo = B [i]
        while not(j>=m):
            if (i!=j):
                nuevo = nuevo - A[i,j]*X[j]
            j = j+1
        Xi1[i] = nuevo/A[i,i]
        i = i+1
    diferencia = np.abs(X - Xi1)
    errado = np.max(diferencia)
    X = np.copy(Xi1)
    itera = itera +1
    print(Xi1, errado)

# SALIDA
print('iteraciones: ', itera-1)

Tarea: Realice las modificaciones necesarias para usar el algoritmo de Gauss-Seidel. Luego compare resultados.


Revisión de resultados

Si se usan las ecuaciones sin cambiar a diagonal dominante, como fueron presentadas en el enunciado, el algoritmo Jacobi NO converge, el error aumenta, y muy rápido como para observar la espiral hacia afuera en una gráfica.

Xi, errado
[ -0.5   6.   26. ] 26.0
[ 131.5  -16.5   69.5] 132.0
[ 271. -984. -484.] 967.5
[-4638.5 -1407.  -7675. ] 7191.0
[-38055.5  40150.5   4092.5] 41557.5
[ 118792.  262302.  395246.] 391153.5
[ 2434361.5 -1226784.   1479764. ] 2315569.5
[  3591977.5 -18520288.5 -15890546.5] 17370310.5
....

s1Eva_IT2018_T1 Tanque esférico canchas deportivas

a) Para seleccionar el rango para h=[a,b], se observa que el tanque puede estar vacio, a=0 o lleno al máximo, b=2R = 2(3)=6, con lo que se obtiene:

h =[0.0, 6.0]

conociendo la proporción con el valor máximo, se tiene un valor inicial para h0 para las iteraciones.

Vmax = \frac{\pi}{3} (2R)^2 (3R-2R) Vmax = \frac{4\pi }{3}R^{3}= 113.09 h_0 = (6)*30/113.09 = 1.59

b) Usar Newton-Raphson con tolerancia 1e-6

f(h) = \frac{\pi }{3}h^2 (3(3)-h)-30 f(h) = \frac{\pi }{3}(9h^2 -h^3-28.647) f'(h) = \frac{\pi }{3}(18h-3h^2)

el punto siguiente de iteración es:

h_{i+1} = h_{i} -\frac{f(h)}{f'(h)} = h_{i}-\frac{ \frac{\pi }{3}(9h^2 -h^3-28.647)}{ \frac{\pi }{3}(18h-3h^2)} h_{i+1} = h_{i} -\frac{(9h^2 -h^3-28.647)}{(18h-3h^2)}

con lo que se realiza la tabla de iteraciones:

hi hi+1 error orden
1.590 2.061 0.47 10-1
2.061 2.027 -0.034 10-2
2.027 2.02686 -0.00014 10-4
2.02686 2.0268689 -2.32E-09 10-9

En valor práctico 2.028 m usando flexómetro, a menos que use medidor laser con precisión 10-6 usará más dígitos con un tanque de 6 metros de altura ganará una precisión de una gota de agua para usar en duchas o regar el césped .

c) El orden de convergencia del error observando las tres primeras iteraciones es cuadrático

Tarea: Realizar los cálculos con Python, luego aplique otro método. Añada sus observaciones y conclusiones.

s1Eva_IT2018_T3 Temperatura en nodos de placa

a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior:

a=\frac{50+c+100+b}{4} b=\frac{a+30+50+d}{4} c=\frac{a+60+100+d}{4} d=\frac{b+60+c+30}{4}

que reordenando se convierte en:

4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90

simplificando:

4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90

que a forma matricial se convierte en:

A = [[ 4, -1, -1, 0.0],
     [ 1, -4,  0, 1.0],
     [ 1,  0, -4, 1.0],
     [ 0,  1,  1,-4.0]]
B = [[ 150.0],
     [ -80.0],
     [-160.0],
     [ -90.0]]

Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas.  Se aumentó el punto decimal a los valores de la matriz A y el vector B  para que sean considerados como números reales.

El número de condición es: np.linalg.cond(A) = 3.0

que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es «estable» y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.

b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial

 
X(0) = [[60.0],
        [40], 
        [70],
        [50]] 

reemplazando los valores iniciales en cada ecuación sin cambios.

iteración 1
a=\frac{50+70+100+40}{4} = 65

b=\frac{60+30+50+50}{4} = 47.5 c=\frac{60+60+100+50}{4} = 67.5 d=\frac{40+60+70+30}{4} = 50
X(1) = [[65],
        [47.5],
        [67.5],
        [50]]

vector de error = 
     [|65-60|,
      |47.5-40|,
      |67.5-70|,
      |50-50|]
  = [|5|,
     |7.5|,
     |-2.5|,
     |0|]
errormax = 7.5

iteración 2
a=\frac{50+67.5+100+47.5}{4} = 66.25

b=\frac{65+30+50+50}{4} = 48.75 c=\frac{65+60+100+50}{4} = 68.75 d=\frac{47.5+60+67.7+30}{4} = 51.3
X(2) = [[66.25],
        [48.75],
        [68.75],
        [51.3]]

vector de error = 
       [|66.25-65|,
        |48.75-47.5|,
        |68.75-67.5|,
           |51.3-50|] 
  = [|1.25|,
     |1.25|,
     |1.25|,
     |1.3|]
errormax = 1.3

iteración 3
a=\frac{50+68.75+100+48.75}{4} = 66.875

b=\frac{66.25+30+50+51.3}{4} = 49.38 c=\frac{66.25+60+100+51.3}{4} = 69.3875 d=\frac{48.75+60+68.75+30}{4} = 51.875
X(2) = [[66.875],
        [49.38],
        [69.3875],
        [51.875]]

vector de error = 
      [|66.875-66.25|,
       |49.38-48.75|,
       |69.3875-68.75|,
       |51.875-51.3|]
    = [|0.655|,
       |0,63|,
       |0.6375|,
        |0.575|]
errormax = 0.655
con error relativo de:
100*0.655/66.875 = 0.97%

siguiendo las iteraciones se debería llegar a:

>>> np.linalg.solve(A,B)
array([[ 67.5],
       [ 50. ],
       [ 70. ],
       [ 52.5]])