s1Eva_IT2016_T3_MN Tasa interés anual

Ejercicio: 1Eva_IT2016_T3_MN Tasa interés anual

Propuesta de Solución  empieza con el planteamiento del problema, luego se desarrolla con el método de Bisección y método del Punto Fijo solo con el objetivo de comparar resultados.


Planteamiento del problema

La fórmula del enunciado para el problema es:

A = P \frac{i(1+i)^{n}}{(1+i)^{n} -1}

que con los datos dados se convierte a:

5800 = 35000 \frac{i(1+i)^8}{(1+i)^8 -1} 35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0

que es la forma de f(x) = 0

f(i)=35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800

Intervalo de búsqueda

Como el problema plantea la búsqueda de una tasa de interés, consideramos que:

  • Las tasas de interés no son negativas. ⌉(i<0)
  • Las tasas de interés no son cero en las instituciones bancarias (i≠0)
  • Las tasas muy grandes 1 = 100/100 = 100% tampoco tienen mucho sentido

permite acotar la búsqueda a un intervalo (0,1].
Sin embargo tasas demasiado altas tampoco se consideran en el problema pues el asunto es regulado (superintendencia de bancos), por lo que se podría intentar entre un 1% = 0.01l y la mitad del intervalo 50%= 0.5 quedando

[0.01,0.5]

Tolerancia

si la tolerancia es de tan solo menor a un orden de magnitud que el valor buscado, se tiene que las tasas de interés se representan por dos dígitos después del punto decimal, por lo que la tolerancia debe ser menor a eso.

Por ejemplo: tolerancia < 0.001 o aumentando la precisión

tolera = 0.0001


Método de la Bisección

itera = 1

a = 0.01, b = 0.5

c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 f(0.01) = 35000 \frac{0.01(1+0.01)^8}{(1+0.01)^8-1}-5800 = -1225.83 f(0.5) = 35000 \frac{0.5(1+0.5i)^8}{(1+0.5)^8-1}-5800 = 12410.54

con lo que se verifica que existe cambio de signo al evaluar f(x) en el intervalo y puede existir una raíz.

f(0.255) = 35000 \frac{0.255(1+0.255)^8}{(1+0.255)^8-1}-5800 = 4856.70

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.255]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.255-0.01 = 0.245

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 2

a = 0.01, b = 0.255

c = \frac{0.01+0.255}{2} = 0.1325 f(0.01) = -1225.83 f(0.255) = 4856.70 f(0.1325 ) = 35000 \frac{0.1325(1+0.1325)^8}{(1+0.1325)^8-1}-5800 = 1556.06

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.1325]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.1325-0.01 = 0.1225

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 3

a = 0.01, b = 0.1325

c = \frac{0.01+0.1225}{2} = 0.071 f(0.01) = -1225.83 f(0.1325) = 1556.06 f(0.071 ) = 35000 \frac{0.071(1+0.071)^8}{(1+0.071)^8-1}-5800 = 89.79

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.071]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.071-0.01 = 0.061

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

Para una evaluación del tema en forma escrita es suficiente para mostrar el objetivo de aprendizaje, el valor final se lo encuentra usando el algoritmo.

       raiz en:  0.06724243164062502
error en tramo:  5.981445312500111e-05
iteraciones:  13
>>>

Algoritmo en Python

# 1Eva_IT2016_T3_MN Tasa interés anual
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: 35000*(x*(1+x)**8)/((1+x)**8 -1) -5800
a = 0.01
b = 0.5
tolera = 0.0001

# PROCEDIMIENTO
cuenta = 0
np.set_printoptions(precision=3)
tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = fx(a)
    fb = fx(b)
    fc = fx(c)
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: a = a b = c if cambia > 0:
        a = c
        b = b
    tramo = b-a
    cuenta = cuenta+1

# SALIDA
print('       raiz en: ', c)
print('error en tramo: ', tramo)
print('iteraciones: ',cuenta)


Método del Punto Fijo

El planteamiendo del punto fijo se realiza con x= g(x), por lo que se reordena la ecuación a:

35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0 35000 \frac{i(1+i)^8}{(1+i)^8 -1}=5800 \frac{i(1+i)^8}{(1+i)^8 -1}=\frac{5800}{35000} i=\frac{58}{350} \frac{(1+i)^8 -1} {i(1+i)^8}

con lo que g(x) es:

g(i)=\frac{58}{350} \frac{(1+i)^8 -1} {(1+i)^8}

valor inicial de búsqueda

Para el punto inicial i0 se puede usar uno de los extremos del intervalo propuesto en la sección de planteamiento. Para reducir aun más la búsqueda se pude seleccionar el punto intermedio

 i0 = 0.255

Itera = 1

 i0 = 0.255

g(0.255i)=\frac{58}{350} \frac{(1+0.255)^8 -1} {(1+0.255)^8} = 0.1387

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1387 - 0.255| = 0.1163

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 2

i1 = g(i0 ) = 0.1387

g(0.1387)=\frac{58}{350} \frac{(1+0.1387)^8 -1} {(1+0.1387)^8} = 0.1071

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1071 - 0.1387| = 0,0316

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 3

i2 = g(i1 ) = 0.1071

g(0.1071)=\frac{58}{350} \frac{(1+0.1071)^8 -1} {(1+0.1071)^8} = 0.0922

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.0922 - 0.1071| = 0,0149

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Observación: Como el error disminuye entre cada iteración, se considera que el método converge, si se realizan suficientes iteraciones se cumplierá con el valor de tolerancia y se habrá llegado a la precisión requerida.

Usando el algoritmo se tiene:

iteraciones: 17
    raiz en: 0.06754389199556779
>>>

Algoritmo en Python

# Algoritmo de Punto Fijo
# x0 valor inicial de búsqueda
# error = tolera

import numpy as np
def puntofijo(gx,antes,tolera, iteramax=50):
    itera = 1 # iteración
    nuevo = gx(antes)
    tramo = abs(nuevo-antes)
    while(tramo>=tolera and itera<=iteramax ):
        antes = nuevo
        nuevo = gx(antes)
        tramo = abs(nuevo-antes)
        itera = itera+1
    respuesta = nuevo
    # Validar respuesta
    if (itera>=iteramax ):
        respuesta = np.nan
    print('iteraciones:',itera)
    return(respuesta)

# PROGRAMA ---------

# INGRESO
gx = lambda i: (58/350)*((1+i)**8-1)/((1+i)**8)

a = 0       # intervalo
b = 0.5

x0 = 0.255
tolera = 0.0001
iteramax  = 50      # itera máximo

# PROCEDIMIENTO
respuesta = puntofijo(gx,0.255,tolera,iteramax)

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

s1Eva_IIT2009_T1 Movimiento de partícula en plano

a. Planteamiento del problema

Las ecuaciones expresan las trayectorias de dos partículas,

x(t) = 3 \sin ^{3}(t)-1 y(t) = 4 \sin (t)\cos (t)

que para que se encuentren o choquen, sus coordenadas deberían ser iguales.

x(t) = y(t) 3 \sin ^{3}(t)-1 = 4 \sin (t)\cos (t)

Se reordena la igualdad, de tal manera que uno de los lados sea cero, de la forma f(t) = 0 que es usada en el algoritmo de búsqueda de raíces.

3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t) = 0 f(t) = 3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t)

b. Intervalo de búsqueda de raíz

Como la variable independiente es tiempo, el evento a buscar se suponte sucede en tiempos positivos t>=0, por lo que el valor inicial a la izquierda del intervalo será a=0

Para el caso de b, a la derecha, se usa lo indicado en el enunciado para la pregunta dell literal b), donde se indica t ∈ [0, π/2], por lo que b = π/2

[0, π/2]

verificando que exista cambio de signo entre f(a) y f(b)

f(0) = 3 \sin ^{3}(0)-1 - 4 \sin (0)\cos (0) = -1 f(\pi /2) = 3 \sin ^{3}(\pi /2)-1 - 4 \sin (\pi /2)\cos (\pi /2) = 3 (1)^3-1 - 4 (1)(0) = 2

con lo que se comprueba que al existir cambio de signo, debe existir una raíz en el intervalo.


c. Método de Falsa Posición


Desarrollo analítico con lápiz y papel

se trata de mostrar los pasos en al menos tres iteraciones del método, usando las siguientes expresiones:

f(t) = 3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t) c = b - f(b) \frac{a-b}{f(a)-f(b)}

[0, π/2]

iteración 1

a = 0 , b = \pi/2

tomando los datos al validar los puntos extremos

f(0) = -1 f(\pi /2) = 2 c = \pi/2 - 2 \frac{0-\pi/2}{-1-2} = \pi/6 f(\pi/6) = 3 \sin ^{3}(\pi/6)-1 - 4 \sin (\pi/6)\cos (\pi/6) = -2.3570

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = |\pi /6 - 0| = \pi/6 a = c = \pi/6 , b = \pi/2

iteración 2

a = \pi/6 , b = \pi/2 f(\pi/6) = -2.3570 f(\pi /2) = 2 c = \pi/2 - 2 \frac{\pi/6-\pi/2}{-1-2} = 1.0901 f(1.0901) = 3 \sin ^{3}(1.0901)-1 - 4 \sin (1.0901)\cos (1.0901) = -0.5486

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = | 1.0901-\pi/6| = 1.0722 a = c = 1.0901 , b = \pi/2

iteración 3

a = 1.0901 , b = \pi/2 f(1.0901) = -0.5486 f(\pi /2) = 2 c = \pi/2 - 2 \frac{-0.5486-\pi/2}{-0.5486-2} = 1.19358 f(1.19358) = 3 \sin ^{3}(1.19358)-1 - 4 \sin (1.19358)\cos (1.19358) = 0.0409

el signo de f(c) es el mismo que f(b), se ajusta el lado derecho

tramo = |b-c| = | \pi/2- 1.19358| = 0.3772 a = 1.0901 , b = 1.19358


Algoritmo con Python

Los parámetros aplicados en el algoritmo son los desarrollados en el planteamiento del problema e intervalo de búsqueda, con lo que se obtiene los siguientes resultados:

 raiz: 1.1864949811467547
error: 9.919955449055884e-05

las instrucciones en python son:

# 1Eva_IIT2009_T1 Movimiento de partícula en plano
import numpy as np
import matplotlib.pyplot as plt

#INGRESO
xt = lambda t: 3*(np.sin(t)**3)-1
yt = lambda t: 4*np.sin(t)*np.cos(t)
fx = lambda t: 3*(np.sin(t)**3)-1 - 4*np.sin(t)*np.cos(t) 

a = 0
b = np.pi/2
tolera = 0.001

# intervalo para gráfica
La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Posicion Falsa
tramo = abs(b-a)
fa = fx(a)
fb = fx(b)
while not(tramo<=tolera):
    c = b - fb*(a-b)/(fa-fb)
    fc = fx(c)
    cambio = np.sign(fa)*np.sign(fc)
    if cambio>0:
        # actualiza en izquierda
        tramo = abs(c-a)
        a = c
        b = b
        fa = fc
    else:
        # actualiza en derecha
        tramo = abs(b-c)
        a = a
        b = c
        fb = fc

# para grafica
ti = np.linspace(La,Lb,muestras)
xi = xt(ti)
yi = yt(ti)
fi = fx(ti)

# SALIDA
print(' raiz:', c)
print('error:', tramo)
plt.plot(ti,xi, label='x(t)')
plt.plot(ti,yi, label='y(t)')
plt.plot(ti,fi, label='f(t)')
plt.plot(c,fx(c),'ro')
plt.axhline(0, color='green')
plt.axvline(c, color='magenta')
plt.legend()
plt.xlabel('t')
plt.title('Método de Falsa Posición')
plt.show()

s1Eva_IIT2010_T1 Aproximar con polinomio

Desarrollo Analítico

Ejemplo para Lápiz  y Papel


Tarea 01 Semana 01 Fecha: año/mes/dia
Apellidos Nombres
Referencia: 1Eva_IIT2010_T1 Aproximar con polinomio

Para el ejemplo, supondremos que x0=0

El polinomio de Taylor requerido es de grado 2

P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 +

función f(x) y sus derivadas:

f(x) = e^x \cos (x) +1
f'(x) = e^x \cos (x) - e^x \sin(x) f'(x) = e^x (\cos (x) - \sin(x))
f''(x) = e^x( \cos (x) - \sin(x))+ + e^x (-\sin(x) - \cos(x)) f''(x) = -2 e^x \sin(x))

Punto x0 = 0 (ejemplo), dentro del intervalo.

Observación: escriba las expresiones, reemplazando los valores, asi si en la lección o examen no tuvo tiempo para usar la calculadora, se puede evaluar si realizaba las operaciones con el punto de referencia y expresiones correctas.

f(0) = e^0 \cos (0) +1 = 2 f'(0) = e^0(\cos (0) - \sin(0)) = 1 f''(0) = -2 e^0 \sin(0)) = 0

Sustitución en fórmula de polinomio:

P_{2}(x) = 2+\frac{1}{1} (x-0) + + \frac{0}{2}(x-0)^2 + P_{2}(x) = 2+ x

Desarrollo con Python

Algoritmo desarrollado en clase, usado como taller, modificado para el problema planteado.

Observación: Se reordena el algoritmo para mantener ordenados y separados los bloques de ingreso, procedimiento y salida. Así los bloques pueden ser convertidos fácilmente a funciones algoritmicas def-return.

Observe que la variable n se interprete correctamente como «términos» o «grados» del polinomio de Taylor.

# Aproximación Polinomio de Taylor alrededor de x0
# f(x) en forma simbólica con sympy
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO --------------------
x = sym.Symbol('x')
fx = sym.exp(x)*sym.cos(x) + 1

x0 = 0
n  = 3 # grado de polinomio

# Intervalo para Gráfica
a = 0
b = np.pi
muestras = 21

# PROCEDIMIENTO  -------------
# construye polinomio Taylor
k = 0 # contador de términos
polinomio = 0
while (k <= n):
    derivada   = fx.diff(x,k)
    derivadax0 = derivada.subs(x,x0)
    divisor   = np.math.factorial(k)
    terminok  = (derivadax0/divisor)*(x-x0)**k
    polinomio = polinomio + terminok
    k = k + 1

# forma lambda para evaluación numérica
fxn = sym.lambdify(x,fx,'numpy')
pxn = sym.lambdify(x,polinomio,'numpy')

# evaluar en intervalo para gráfica
xi = np.linspace(a,b,muestras)
fxi = fxn(xi)
pxi = pxn(xi)

# SALIDA  --------------------
print('polinomio p(x)=')
print(polinomio)
print()
sym.pprint(polinomio)

# Gráfica
plt.plot(xi,fxi,label='f(x)')
plt.plot(xi,pxi,label='p(x)')
# franja de error
plt.fill_between(xi,pxi,fxi,color='yellow')
plt.xlabel('xi')
plt.axvline(x0,color='green', label='x0')
plt.axhline(0,color='grey')
plt.title('Polinomio Taylor: f(x) vs p(x)')
plt.legend()
plt.show()

Resultado del algoritmo

Revisar si el polinomio es concordante con lo realizado a lápiz y papel, de no ser así revisar el algoritmo o los pasos realizados en papel, deben ser iguales.
Comprobando que el algoritmo esté correcto y pueda ser usado en otros ejercicios.

 RESTART: D:\MATG1052Ejemplos\Taylot01Tarea01.py 
polinomio p(x)=
-x**3/3 + x + 2

   3        
  x         
- -- + x + 2
  3         

Resultados gráficos para x0=0

Continuar con el ejercicio con x0 = π y luego con el siguiente punto x0 = π/2.

Comparar resultados y presentar: Observaciones  y recomendaciones semejantes a las indicadas durante el desarrollo de la clase.

s1Eva_IIT2019_T1 Ecuación Recursiva

la ecuación recursiva es:

x_n = g(x_{n-1}) = \sqrt{3 + x_{n-1}}

literal a y b

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

Verifique la cota de g'(x)

g(x) = \sqrt{3 + x} =(3+x)^{1/2} g'(x) =\frac{1}{2}(3+x)^{-1/2} g'(x) =\frac{1}{2\sqrt{3+x}}

Tarea: verificar que g'(x) es menor que 1 en todo el intervalo.

Literal c

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

Iteración 1: x0=2

g(x_0) = \sqrt{3 + 2} = 2.2361

error = |2.2361 – 2| = 0.2361

Iteración 2: x1 = 2.2361

g(x_1) = \sqrt{3 + 2.2361} = 2.2882

error = |2.2882 – 2.2361| = 0.0522

Iteración 3: x2 = 2.2882

g(x_2) = \sqrt{3 + 2.2882} = 2.2996

error = |2.2996 – 2.28821| = 0.0114

Iteración 4: x3 = 2.2996

g(x_3) = \sqrt{3 + 2.2996} = 2.3021

error = |2.3021- 2.2996| = 0.0025

Iteración 5: x4 = 2.3021

g(x_4) = \sqrt{3 + 2.3021} = 2.3026

error = |2.3021- 2.2996| = 5.3672e-04

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

Se muestra como referencia la tabla resumen.

[[ x ,   g(x), 	 tramo  ] 
 [1.      2.      1.    ]
 [2.      2.2361  0.2361]
 [2.2361  2.2882  0.0522]
 [2.2882  2.2996  0.0114]
 [2.2996  2.3021  0.0025]
 [2.3021  2.3026  5.3672e-04]
 [2.3026  2.3027  1.1654e-04]
 [2.3027  2.3028  2.5305e-05]
raiz:  2.3027686193257098

con el siguiente comportamiento de la funcion:

literal e

Realizando el mismo ejercicio para el método de la bisección, se requiere cambiar a la forma f(x)=0

x = \sqrt{3 + x} 0 = \sqrt{3 + x} -x f(x) = \sqrt{3 + x} -x

tomando como intervalo el mismo que el inicio del problema [1,3], al realizar las operaciones se tiene que:

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

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

a = 2
b = 2.5
...

Siguiendo las operaciones se obtiene la siguiente tabla:

[ i, a,   c,   b,    f(a),  f(c),  f(b), paso]
 1 1.000 2.000 3.000 1.000  0.236 -0.551 2.000 
 2 2.000 2.500 3.000 0.236 -0.155 -0.551 1.000 
 3 2.000 2.250 2.500 0.236  0.041 -0.155 0.500 
 4 2.250 2.375 2.500 0.041 -0.057 -0.155 0.250 
 5 2.250 2.312 2.375 0.041 -0.008 -0.057 0.125 
 6 2.250 2.281 2.312 0.041  0.017 -0.008 0.062 
 7 2.281 2.297 2.312 0.017  0.005 -0.008 0.031 
 8 2.297 2.305 2.312 0.005 -0.001 -0.008 0.016 
 9 2.297 2.301 2.305 0.005  0.002 -0.001 0.008 
10 2.301 2.303 2.305 0.002  0.000 -0.001 0.004 
11 2.303 2.304 2.305 0.000 -0.001 -0.001 0.002 
12 2.303 2.303 2.304 0.000 -0.000 -0.001 0.001 
13 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
14 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
15 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
16 2.303 2.303 2.303 0.000  0.000 -0.000 0.000 
raiz:  2.302764892578125

Donde se observa que para la misma tolerancia de 10-4, se incrementan las iteraciones a 16. Mientra que con punto fijo eran solo 8.

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

s1Eva_IIT2019_T3 Circuito eléctrico

Las ecuaciones del problema  son:

55 I_1 - 25 I_4 =-200 -37 I_3 - 4 I_4 =-250 -25 I_1 - 4 I_3 +29 I_4 =100 I_2 =-10

Planteo del problema en la forma A.X=B

A = np.array([[ 55, 0,   0,-25],
              [  0, 0, -37, -4],
              [-25, 0,  -4, 29],
              [  0, 1,   0,  0]],dtype=float)

B = np.array([[-200],
              [-250],
              [ 100],
              [ -10]],dtype=float)

El ejercicio se puede simplificar con una matriz de 3×3 dado que una de las corrientes I2 es conocida con valor -10, queda resolver el problema para
[I1 ,I3 ,I4 ]

A = np.array([[ 55,   0,-25],
              [  0, -37, -4],
              [-25,  -4, 29]],dtype=float)

B = np.array([[-200],
              [-250],
              [ 100]],dtype=float)

conformando la matriz aumentada

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

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

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

Literal a

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

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

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

luego se continúa con la segunda columna:

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

y para el método de Gauss se emplea sustitución hacia atras
se determina el valor de I4

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

y planteando se obtiene el último valor

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

con lo que el vector solución es:

[[-2.7277672]
 [ 6.54065815]
 [ 1.99891216]]

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

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

literal b
La norma de la matriz infinito se determina como:

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

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

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

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

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

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

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

Tarea: Matriz de transición de Jacobi


Literal c

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

las ecuaciones para el método son:

I_1 =\frac{-200 + 25 I_4}{55} I_3 = \frac{-250+ 4 I_4}{-37} I_4 =\frac{100 +25 I_1 + 4 I_3}{29}

Como I2 es constante, no se usa en las iteraciones

I_2 =-10

teniendo como resultados de las iteraciones:

iteración:  1
 Xnuevo:      [-3.63636364  6.75675676  1.99891216]
 diferencia:  [3.63636364   6.75675676  1.99891216]
 error:  6.756756756756757
 iteración:  2
 Xnuevo:      [-2.7277672   6.54065815  1.99891216]
 diferencia:  [0.90859643   0.21609861  0.        ]
 error:  0.9085964348406557
 iteración:  3
 Xnuevo:      [-2.7277672   6.54065815  1.99891216]
 diferencia:  [0. 0. 0.]
 error:  0.0

con lo que el vector resultante es:

 Método de Gauss-Seidel
respuesta de A.X=B : 
[[-2.7277672 ]
 [ 6.54065815]
 [ 1.99891216]]

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

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


Solución alterna


Usando la matriz de 4×4, los resultados son iguales para las corrientes
[I1 ,I2 , I3 ,I4 ]. Realizando la matriz aumentada,

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

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

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

Literal a

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

Usando el primer elemento  en la diagonal.

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

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

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

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

18.06879607 x_4 = 36.11793612 x_4 = 1.99891216

para x3:

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

como ejercicio, continuar con x1, dado que x2=-10

55 x_1 + 25 x_4 = -200

El vector solución obtenido es:

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

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

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

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

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

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


Algoritmo en Python

Presentado por partes para revisión:

# Método de Gauss,
# Recibe las matrices A,B
# presenta solucion X que cumple: A.X=B
import numpy as np

# INGRESO
A = np.array([[ 55,   0,-25],
              [  0, -37, -4],
              [-25,  -4, 29]],dtype=float)

B = np.array([[-200],
              [-250],
              [ 100]],dtype=float)

# PROCEDIMIENTO
# Matriz Aumentada
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
print('aumentada: ')
print(AB)

# pivotea por fila AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = np.abs(AB[i:,i])
    dondemax = np.argmax(columna)
    # revisa dondemax no está en la diagonal
    if (dondemax != 0):
        # intercambia fila
        temporal = np.copy(AB[i,:])
        AB[i,:] = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal
print('pivoteada: ')
print(AB)

# Gauss elimina hacia adelante
# tarea: verificar términos cero
print('Elimina hacia adelante')
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1
    print('i: ', i)
    for k in range(adelante,n,1):
        print('k: ',k)
        if (np.abs(pivote)>=casicero):
            factor = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - factor*AB[i,:]
        print(AB)
print('Elimina hacia adelante')
print(AB)

# Sustitución hacia atras
X = np.zeros(n,dtype=float) 
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    suma = 0
    for j in range(i+1,ultcolumna,1):
        suma = suma + AB[i,j]*X[j]
    X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
X = np.transpose([X])

# Verifica resultado
verifica = np.dot(A,X)

# SALIDA
print('por sustitución hacia atras')
print('el vector solución X es:')
print(X)

print('verificar que A.X = B')
print(verifica)

# -------------------------------------------------
# # Algoritmo Gauss-Seidel,
# matrices, métodos iterativos
# Referencia: Chapra 11.2, p.310, pdf.334
#      Rodriguez 5.2 p.162
# ingresar iteramax si requiere más iteraciones

import numpy as np

def gauss_seidel(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)
    
    itera = 0
    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]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)
        print(' iteración: ', itera)
        print(' Xnuevo: ', X)
        print(' diferencia: ', diferencia)
        print(' error: ' ,errado)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=0
    return(X)

# Programa de prueba #######
tolera = 0.0001

# PROCEDIMIENTO
# usando la pivoteada por filas
n = len(B)
Xgs = np.zeros(n, dtype=float)
respuesta = gauss_seidel(AB[:,:n],AB[:,n],tolera,Xgs)
verificags = np.dot(A,respuesta)

# SALIDA
print(' Método de Gauss-Seidel')
print('respuesta de A.X=B : ')
print(respuesta)
print('verificar A.X: ')
print(verificags)

s1Eva_IIT2019_T2 Proceso Termodinámico

la ecuación para el problema se describe como:

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

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

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

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

literal a

Con los puntos seleccionados se escriben las ecuaciones del polinomio:

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

usando los valores de la tabla:

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

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

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

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

matriz aumentada

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

matriz pivoteada

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

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

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

[[ 0.04994705]
 [-0.41595438]
 [ 1.        ]]

con lo que el polinomio de interpolación es:

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

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

p_2(2.4) = 0.04994705 (2.4)^2 - 0.41595438 (2.4) + 1.0 f(2.4)=e^{-0.5(2.4)} error = |f(2.4)-p_2(2.4)|
Evaluando p(2.4):  0.2894044975129779
Error en 2.4:      0.011789714399224216
Error relativo:    0.039143230291095066

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


literal b

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


Algoritmo en Python

Presntado por secciones, semejante a lo desarrollado en clases

# Interpolación por polinomio
import numpy as np
import matplotlib.pyplot as plt

# INGRESO, Datos de prueba
fx = lambda x: np.exp(-0.5*x)

xi = np.array([0,2,4])

# determina vector
n = len(xi)
fi = np.zeros(n,dtype=float)
for i in range(0,n,1):
    fi[i]= fx(xi[i])

# PROCEDIMIENTO
# Arreglos numpy 
xi = np.array(xi)
fi = np.array(fi)
n = len(xi)

# Vector B en columna
B = np.array(fi)
B = fi[:,np.newaxis]

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for f in range(0,n,1):
    for c in range(0,n,1):
        potencia = (n-1)-c
        D[f,c] = xi[f]**potencia

# Resolver matriz aumentada
coeficiente =  np.linalg.solve(D,B)

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

# Convierte polinomio a funcion 
# para evaluación más rápida
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio
k = 100
inicio = np.min(xi)
fin = np.max(xi)
puntosx = np.linspace(inicio,fin,k)
puntosy = px(puntosx)

puntosf = fx(puntosx)

# Puntos evaluados
x1 = 2.4
px1 = px(x1)
fx1 = fx(x1)
errorx1 = np.abs(px1-fx1)
errorx1rel = errorx1/fx1

x2 = 1.7
px2 = px(x1)
fx2 = fx(x1)
errorx2 = np.abs(px1-fx1)
errorx2rel = errorx1/fx1
    
# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print()
print('Evaluando en X1: ',x1)
print('Evaluando p(x1): ',px1)
print('Error en x1:     ',errorx1)
print(' Error relativo: ', errorx1rel)
print()
print('Evaluando en X2: ',x2)
print('Evaluando p(x2): ',px2)
print('Error en x2:     ',errorx2)
print(' Error relativo: ', errorx2rel)

# Grafica
plt.plot(xi,fi,'o',label='muestras' )
plt.plot(puntosx,puntosy,label='p(x)')
plt.plot(puntosx,puntosf,label='f(x)')
plt.xlabel('x')
plt.legend()
plt.show()

s1Eva_IIT2019_T4 Concentración de químico

formula a usar:

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

Se sustituyen los valores dados con:
C0 = 4, Cent = 10, C = 0.93 Cent.

0.93(10) = 10 ( 1 - e^{-0.04t}) + 4 e^{-0.03t}

igualando a cero para forma estandarizada del algoritmo,

10( 1 - e^{-0.04t}) + 4 e^{-0.03t} - 9.3 = 0 0.7 - 10 e^{-0.04t} + 4 e^{-0.03t} = 0

se usas las funciones f(t) y f'(t) para el método de Newton-Raphson,

f(t) = 0.7 - 10 e^{-0.04t} + 4 e^{-0.03t} f'(t) = - 10(-0.04) e^{-0.04t} + 4(-0.03) e^{-0.03t} f'(t) = 0.4 e^{-0.04t} - 0.12 e^{-0.03t}

con lo que se pueden realizar los cálculos de forma iterativa.

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

De no disponer de la gráfica de f(t), y se desconoce el valor inicial para t0 se usa 0. Como no se indica la tolerancia, se estima en 10-4

Iteración 1

t0 = 0

f(0) = 0.7 - 10 e^{-0.04(0)} + 4 e^{-0.03(0)} = 5.3 f'(0) = 0.4 e^{-0.04(0)} - 0.12 e^{-0.03(0)} = -0.28 t_{1} = 0 -\frac{5.3}{-0.28} = 18.92 error = |18.92-0| = 18.92

Iteración 2

f(18.92) = 0.7 - 10 e^{-0.04(18.92)} + 4 e^{-0.03(18.92)} = -1.72308 f'(18.92) = 0.4 e^{-0.04(18.92)} - 0.12 e^{-0.03(18.92)} = 0.119593 t_{2} = 18.92 -\frac{-1.723087}{0.119593} = 33.3365 error = |33.3365 - 18.92| = 14.4079

Iteración 3

f(33.3365) = 0.7 - 10 e^{-0.04(33.3365)} + 4 e^{-0.03(33.3365)} = -0.4642 f'(33.3365) = 0.4 e^{-0.04(33.3365)} - 0.12 e^{-0.03(33.3365)} = 0.06128 t_{3} = 33.3365 -\frac{-0.46427}{-5.8013} = 40.912 error = |40.912 - 33.3365| = 7.5755

Observando que los errores disminuyen entre cada iteración, se encuentra que el método converge.

y se forma la siguiente tabla:

['xi' ,  'xnuevo', 'error']
[ 0.      18.9286  18.9286]
[18.9286  33.3365  14.4079]
[33.3365  40.912    7.5755]
[40.912   42.654     1.742]
[42.654   42.7316   0.0776]
[4.2732e+01 4.2732e+01 1.4632e-04]
raiz en:  42.731721341402796

Observando la gráfica de la función puede observar el resultado:


Algoritmo en Python

# 1Eva_IIT2019_T4
# Método de Newton-Raphson

import numpy as np
import matplotlib.pyplot as plt

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

x0 = 0
tolera = 0.001

a = 0
b = 60
muestras = 21

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

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

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

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

# grafica
plt.plot(xk,fk)
plt.axhline(0, color='black')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.show()

s1Eva_IT2019_T3 Vector perpendicular a plano

Literal a

Para que un vector sea perpendicular a otro, se debe cumplir que
V1.V2 =0.

\begin{bmatrix} 2 \\ -3 \\ a \end{bmatrix} . \begin{bmatrix} b \\ 1 \\ -4 \end{bmatrix} = 0

se obtiene entonces la ecuación:

(2)(b)+(-3)(1)+(a)(-4)=0
2b -3 -4a =0

procediendo de la misma forma con los siguientes pares de vectores:

\begin{bmatrix} 2 \\ -3 \\ a \end{bmatrix} . \begin{bmatrix} 3 \\ c \\ 2 \end{bmatrix} = 0

se obtiene entonces la ecuación:

(2)(3)+(-3)(c)+(a)(2)=0
6 -3c +2a = 0
\begin{bmatrix} b \\ 1 \\ -4 \end{bmatrix} . \begin{bmatrix} 3 \\ c \\ 2 \end{bmatrix} = 2

se obtiene entonces la ecuación:

(b)(3)+(1)(c)+(-4)(2)=2
3b +c -8 =2

se obtiene el sistema de ecuaciones:

\begin{cases}-4a+2b=3 \\ 2a-3c=-6 \\3b+c=10 \end{cases}

Literal b

Se convierte a la forma matricial Ax=B

\begin{bmatrix} -4 && 2 && 0 \\ 2 && 0 && -3 \\ 0 && 3 && 1\end{bmatrix}.\begin{bmatrix} a \\b\\c \end{bmatrix} = \begin{bmatrix} 3 \\ -6 \\ 10 \end{bmatrix}

se crea la matriz aumentada:

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 2 && 0 && -3 &&-6 \\ 0 && 3 && 1 && 10 \end{bmatrix}

se pivotea por filas buscando hacerla diagonal dominante:

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 2 && 0 && -3 &&-6 \end{bmatrix}

se aplica el algoritmo de eliminación hacia adelante:
1ra Iteración

la eliminación del primer término columna es necesario solo para la tercera fila:

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 2 -(-4)\Big( \frac{2}{-4}\Big) && 0-2\Big(\frac{2}{-4}\Big) && -3 -0\Big(\frac{2}{-4}\Big) &&-6 -3\Big(\frac{2}{-4}\Big) \end{bmatrix} \begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 0 && 1 && -3 && -4.5 \end{bmatrix}

2da Iteración

\begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 0 && 1 -3\Big(\frac{1}{3}\Big) && -3-(1)\Big(\frac{1}{3}\Big) && -4.5-10\Big(\frac{1}{3}\Big) \end{bmatrix} \begin{bmatrix} -4 && 2 && 0 && 3 \\ 0 && 3 && 1 && 10 \\ 0 && 0 && -\frac{10}{3} && -7.8333 \end{bmatrix}

Aplicando eliminación hacia atras

(-10/3)c = -7.8333
c = -7.8333(-3/10) = 2.35

3b +c = 10
b= (10-c)/3 = (10-2.35)/3 = 2.55

-4a +2b =3
a= (3-2b)/(-4) = (3-2(2.55))/(-4) = 0.525

como resultado, los vectores buscados:

v1 = (2,-3,0.525)
v2 = (2.55,1,-4)
v3 = (3,2.35,2)

comprobando resultados:

>>> np.dot((2,-3,0.525),(2.55,1,-4))
-4.440892098500626e-16
>>> np.dot((2,-3,0.525),(3,2.35,2))
-6.661338147750939e-16
>>> np.dot((2.55,1,-4),(3,2.35,2))
2.0

Los dos primeros resultados son muy cercanos a cero, por lo que se los considera válidos

literal c

Para usar el método de Jacobi, se despeja una de las variables de cada ecuación:

\begin{cases} a = \frac{2b -3}{4} \\b = \frac{10-c}{3} \\c = \frac{2a+6}{3} \end{cases}

usando el vector x(0) = [0,0,0]

1ra iteración

a = \frac{2b -3}{4} = \frac{2(0) -3}{4} = -\frac{3}{4} b = \frac{10-c}{3} = \frac{10-0}{3} = \frac{10}{3} c = \frac{2a+6}{3 }= \frac{2(0)+6}{3} = 2

x(1) = [-3/4,10/3,2]
diferencias = [-3/4-0,10/3-0,2-0] = [-3/4,10/3,2]
error = max|diferencias| = 10/3 = 3.3333

2da iteración

a = \frac{2b -3}{4} = \frac{2(10/3) -3}{4} = \frac{11}{12} b = \frac{10-c}{3} = \frac{10-2}{3} = \frac{8}{3} c = \frac{2a+6}{3} = \frac{2(-3/4)+6}{3} = \frac{3}{2}

x(2) = [11/12,8/3,3/2]
diferencias = [11/12-(-3/4),8/3-10/3,3/2-2] = [20/12, -2/3, -1/2]
error = max|diferencias| = 5/3= 1.666

3ra iteración

a = \frac{2b -3}{4} = \frac{2(8/3)-3}{4} = \frac{7}{12} b = \frac{10-c}{3} = \frac{10-3/2}{3} = \frac{17}{6} c = \frac{2a+6}{3} = \frac{2(11/12)+6}{3} = 2.6111

x(3) = [7/12, 17/6, 2.6111]
diferencias = [7/12-11/12, 17/6-8/3, 2.6111-3/2] = [-1/3, 1/6, 1.1111]
error = max|diferencias| = 1.1111

Los errores disminuyen en cada iteración, por lo que el método converge,
si se analiza en número de condición de la matriz A es 2, lo que es muy cercano a 1, por lo tanto el método converge.


Revisión de resultados

Resultados usando algoritmos desarrollados en clase:

matriz aumentada: 
[[-4.   2.   0.   3. ]
 [ 0.   3.   1.  10. ]
 [ 2.   0.  -3.  -6. ]]
Elimina hacia adelante
[[-4.   2.   0.   3. ]
 [ 0.   3.   1.  10. ]
 [ 0.   1.  -3.  -4.5]]
Elimina hacia adelante
[[-4.   2.   0.        3.      ]
 [ 0.   3.   1.       10.      ]
 [ 0.   0.  -3.333333 -7.833333]]
Elimina hacia adelante
[[-4.   2.   0.        3.      ]
 [ 0.   3.   1.       10.      ]
 [ 0.   0.  -3.333333 -7.833333]]
Elimina hacia atras
[[ 1.  -0.  -0.   0.525]
 [ 0.   1.   0.   2.55 ]
 [-0.  -0.   1.   2.35 ]]
el vector solución X es:
[[0.525]
 [2.55 ]
 [2.35 ]]
verificar que A.X = B
[[ 3.]
 [10.]
 [-6.]]

Número de condición de A: 2.005894

los resultados para [a,b,c]:
[0.525 2.55  2.35 ]

producto punto entre vectores:
v12:  0.0
v13:  1.3322676295501878e-15
v23:  2.0

Algoritmos en Python:

# 1Eva_IT2019_T3 Vector perpendicular a plano
import numpy as np

# INGRESO
A = np.array([[-4.,2,0],
              [ 0., 3,1],
              [ 2.,0,-3]])
B = np.array([3.,10,-6])

# PROCEDIMIENTO
B = np.transpose([B])

casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

print('matriz aumentada: ')
print(AB)
# Gauss elimina hacia adelante
# tarea: verificar términos cero
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(pivote)>=casicero):
            factor = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - factor*AB[i,:]
    print('Elimina hacia adelante')
    print(AB)

# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            factor = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*factor - AB[i,:]
        else:
            factor= 'division para cero'
print('Elimina hacia atras')
print(AB)

X = AB[:,ultcolumna]

# Verifica resultado
verifica = np.dot(A,X)

# SALIDA
print('el vector solución X es:')
print(np.transpose([X]))

print('verificar que A.X = B')
print(np.transpose([verifica]))

numcond = np.linalg.cond(A)

# para comprobar respuestas

v1 = [2,-3,X[0]]
v2 = [X[1],1,-4]
v3 = [3,X[2],2]

v12 = np.dot(v1,v2)
v13 = np.dot(v1,v3)
v23 = np.dot(v2,v3)
      
# SALIDA
print('\n Número de condición de A: ', numcond)

print('\n los resultados para [a,b,c]:')
print(X)

print('\n productos punto entre vectores:')
print('v12: ',v12)
print('v13: ',v13)
print('v23: ',v23)

Tarea, usar el algoritmo de Jacobi usado en el taller correspondiente.

s1Eva_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_T2 Catenaria cable2

Desarrollo con error en fórmula ‘y’, en argumento Cosh(), revisar fórmulas en referencia del problema por error de tipografía:
————————-
Las fórmulas con las que se requiere trabajar son:

y = \frac{T_A}{w} cosh \Big( \frac{T_A}{w}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^{\Big( \frac{T_A}{w}x \Big)} + e^{-\Big( \frac{T_A}{w}x \Big)}}{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^{\Big( \frac{T_A}{12}50 \Big)} + e^{-\Big( \frac{T_A}{12}50 \Big)}}{2} + 6 - \frac{T_A}{12}

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

\frac{T_A}{24} e^{\Big( \frac{50}{12} T_A\Big)} + \frac{T_A}{24} e^{-\Big( \frac{50}{12} T_A\Big)} - \frac{T_A}{12} -9 =0

estandarizando a TA/12 =x

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

la función a usar es:

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

Para el méodo 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{x}{2}(50)e^{50x} + \frac{1}{2}e^{50x} + + \frac{x}{2}(-50) e^{-50x} + \frac{1}{2} e^{-50x} - 1 f'(x) = 25x(e^{50x}-e^{-50x}) + + \frac{1}{2}(e^{50x} +e^{-50x}) - 1