Ejercicio: 1Eva_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 expresión, de la forma f(t) = 0 para usarla 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 supone 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 del 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)

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

Ejercicio: 1Eva_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

---

Opción 1

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) = f(x_0) + \frac{f'(x_0)}{1!}(x-x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2 P_{2}(x) = 2+\frac{1}{1} (x-0) + + \frac{0}{2}(x-0)^2 + P_{2}(x) = 2+ x

---

Tarea: realizar el ejercicio para x₀ = π/2, verificando que pase por los puntos requeridos. Respuesta usando el algoritmo de Taylor:

p(x) = -4.81047738096535*x**2 + 10.3020830193353*x - 3.31309698247881

---

Opción 2

El polinomio requerido tiene la forma:

p(x) = a + bx + cx^2

por lo que conociendo los pares ordenados por donde debe pasar se puede plantear las ecuaciones y encontrar a,b,c.

f(0) = e^0 \cos (0) +1 = 2 f(\pi/2) = e^{\pi/2} \cos (\pi /2) +1 = 1(0)+1 =1 f(\pi) = e^{\pi} \cos (\pi) +1 = e^{\pi} +1

se encuentra que a = 2 cuando x = 0 y que reemplazando los valores de x =π/2 y x=π se tiene:

2 + (π/2) b + (π/2)2 c = 1
2 +    π  b +   (π)2 c = eπ +1

que se convierte en:

(π/2) b + (π/2)2 c = -1
   π  b +   (π)2 c = -(eπ +1)

al multiplicar la primera ecuación por 2 y restando de la segunda

- π2/2 c = eπ -1
c = (-2/π2)(eπ -1)

y sustituir c en la segunda ecuación:

π b + (π)² (-2/π²)(e^π -1) = -(e^π +1)
π b = -(e^π +1) + 2(e^π -1) = -e^π -1 + 2e^π -2
b = (e^π -3)/π

El polinomio resultante es:

p(x) = 2 + \frac{e^{\pi}-3}{\pi}x + \frac{-1(e^{\pi}-1)}{\pi ^2}x^2

Probando respuesta con los valores en la función y polinomio usando Python, se encuentra que el polinomio pasa por los puntos. Al observar la gráfica observa que se cumple lo requerido pero visualiza el error de aproximación usando el método de la opción 2.

---

Algoritmo 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 algorítmicas 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.

Ejercicio: 1Eva_IIT2019_T1 Ecuación Recursiva

la ecuación recursiva es:

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

literal a y b

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

Verifique la cota de g'(x)

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

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

Literal c

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

Iteración 1: x₀=2

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

error = |2.2361 – 2| = 0.2361

Iteración 2: x₁ = 2.2361

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

error = |2.2882 – 2.2361| = 0.0522

Iteración 3: x₂ = 2.2882

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

error = |2.2996 – 2.28821| = 0.0114

Iteración 4: x₃ = 2.2996

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

error = |2.3021- 2.2996| = 0.0025

Iteración 5: x₄ = 2.3021

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

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

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

Se muestra como referencia la tabla resumen.

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

con el siguiente comportamiento de la funcion:

literal e

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

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

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

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

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

a = 2
b = 2.5
...

Siguiendo las operaciones se obtiene la siguiente tabla:

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

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

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

Ejercicio: 1Eva_IIT2019_T3 Circuito eléctrico

Las ecuaciones del problema son:

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

Planteo del problema en la forma A.X=B

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

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

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

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

B = [-200,-250,100]

conformando la matriz aumentada

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

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

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

Literal a

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

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

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

luego se continúa con la segunda columna:

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

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

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

y planteando se obtiene el último valor

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

con lo que el vector solución es:

[-2.7277672   6.54065815  1.99891216]

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

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

---

literal b

La norma de la matriz infinito se determina como:

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

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

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

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

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

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

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

Tarea: Matriz de transición de Jacobi

---

Literal c

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

las ecuaciones para el método son:

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

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

I_2 =-10

teniendo como resultados de las iteraciones:

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

con lo que el vector resultante es:

respuesta con Gauss-Seidel
[-2.72777 6.5406588 1.9989098]

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

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

---

---

Solución alterna

---

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

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

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

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

Literal a

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

Usando el primer elemento  en la diagonal.

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

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

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

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

18.06879607 x_4 = 36.11793612 x_4 = 1.99891216

para x₃:

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

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

55 x_1 + 25 x_4 = -200

El vector solución obtenido es:

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

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

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

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

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

Y para éste caso, el número de condición se encuentra alejado del valor 1, contrario a la respuesta del la primera forma de solución con la matriz 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:

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

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

Instrucciones en Python usando las funciones creadas en la unidad:

# 1Eva_IIT2019_T3 Circuito eléctrico
# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

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

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

def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss sustituye hacia atras
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    # Sustitución hacia atras
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1):
        suma = 0
        for j in range(i+1,ultcolumna,1):
            suma = suma + AB[i,j]*X[j]
        X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
    return(X)

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

B = [-200,-250,100]

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

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_sustituyeAtras(AB,vertabla=True)

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

---

literal c

Resultados con el algoritmo de Gauss Seidel

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

Instrucciones en Python

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

def gauss_seidel(A,B,X0,tolera, iteramax=100,
                 vertabla=False, precision=4):
    ''' Método de Gauss Seidel, tolerancia, vector inicial X0
        para mostrar iteraciones: vertabla=True
    '''
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float)
    X0 = np.array(X0, dtype=float)
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = 2*tolera*np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('Iteraciones Gauss-Seidel')
        print('itera,[X]')
        print('      diferencia,errado')
        print(itera, X, errado)
        np.set_printoptions(precision)
    while (errado>tolera and itera<iteramax):
        for i in range(0,n,1):
            xi = B[i]
            for j in range(0,m,1):
                if (i!=j):
                    xi = xi-A[i,j]*X[j]
            xi = xi/A[i,i]
            diferencia[i] = np.abs(xi-X[i])
            X[i] = xi
        errado = np.max(diferencia)
        itera = itera + 1
        if vertabla==True:
            print(itera, X)
            print('  ', diferencia, errado)
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas, entrega matriz aumentada AB
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

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

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

B = [-200,-250,100]

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

tolera = 0.00001
iteramax = 100
verdecimal = 7

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

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

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

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

.

Ejercicio: 1Eva_IIT2019_T2 Proceso Termodinámico

la ecuación para el problema se describe como:

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

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

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

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

literal a

Con los puntos seleccionados se escriben las ecuaciones del polinomio:

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

usando los valores de la tabla:

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

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

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

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

matriz aumentada

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

matriz pivoteada

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

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

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

[ 0.04994705 -0.41595438  1.        ]

con lo que el polinomio de interpolación es:

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

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

p_2(2.4) = 0.04994705 (2.4)^2 - 0.41595438 (2.4) + 1.0 f(2.4)=e^{-0.5(2.4)} error = |f(2.4)-p_2(2.4)|

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

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

---

literal b

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

---

Algoritmo en Python

Resultado con el algoritmo

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

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

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

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

Presentado por secciones, semejante a lo desarrollado en clases

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

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

# determina vector
fi= fx(xi)

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

B = fi
n = len(xi)

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

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

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

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

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


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


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

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

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

 

Ejercicio: 1Eva_IIT2019_T4 Concentración de químico

formula a usar:

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

Se sustituyen los valores dados con:
C₀ = 4, C_ent = 10, C = 0.93 C_ent.

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

igualando a cero para forma estandarizada del algoritmo,

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

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

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

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

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

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

Iteración 1

t₀ = 0

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

Iteración 2

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

Iteración 3

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

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

y se forma la siguiente tabla:

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

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

---

Algoritmo en Python

# 1Eva_IIT2019_T4
# Método de Newton-Raphson

import numpy as np
import matplotlib.pyplot as plt

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

x0 = 0
tolera = 0.001

a = 0
b = 60
muestras = 21

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

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

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

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

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

Ejercicio: 1Eva_IT2019_T3 Vector perpendicular a plano

Literal a

Para que un vector sea perpendicular a otro, se debe cumplir que
V₁.V₂ =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:

v₁ = (2,-3,0.525)
v₂ = (2.55,1,-4)
v₃ = (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]

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⁽¹⁾ = [-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⁽²⁾ = [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⁽³⁾ = [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.

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, y₀=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 raíces:

\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 raíces 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 T_A, considere que el cable colgado no debería tener tensión T_A=0 N, pues en la forma x=12/T_A se crea una indeterminación. Si no dispone de algún criterio para seleccionar el valor de T_A puede iniciar un valor positivo, por ejemplo 120 con lo que el valor de x₀=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 = | x₁ – x₀| = | 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 = | x₂ – x₁| = |0.05098- 0.02428| = 0.02473

Iteración 3

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

error = | x₃ – x₂| = |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:

https://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()

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:

p₃(15) = 10.1006835937500

literal b

se deriva el polinomio del literal anterior y se evalúa 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’₃(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:

tabla de diferencias finitas
['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']
[[ 0.   8.  11.5 -1.6  0.1  0.3  0. ]
 [ 1.  16.   9.9 -1.5  0.4  0.   0. ]
 [ 2.  24.   8.4 -1.1  0.   0.   0. ]
 [ 3.  32.   7.3  0.   0.   0.   0. ]]
dfinita:  [-1.6  0.1  0.3  0. ]
11.5 +
+( -1.6 / 8.0 )* ( x - 8.0 )
+( 0.1 / 128.0 )*  (x - 16.0)*(x - 8.0) 
+( 0.3 / 3072.0 )*  (x - 24.0)*(x - 16.0)*(x - 8.0) 
polinomio simplificado
9.8e-5*x**3 - 0.003923*x**2 - 0.149752*x + 12.898912
Literal a
9.8e-5*x**3 - 0.003923*x**2 - 0.149752*x + 12.898912
p(15) =  10.1007070000000

Literal b
0.000294*x**2 - 0.007846*x - 0.149752
dp(16) = -0.200024000000000
método de Bisección
i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [16, 20.0, 24] [ 0.9     0.1187 -0.6   ]
   4.0
1 [20.0, 22.0, 24] [ 0.1187 -0.2509 -0.6   ]
   2.0
2 [20.0, 21.0, 22.0] [ 0.1187 -0.0683 -0.2509]
   1.0
3 [20.0, 20.5, 21.0] [ 0.1187  0.0246 -0.0683]
   0.5
4 [20.5, 20.75, 21.0] [ 0.0246 -0.022  -0.0683]
   0.25
5 [20.5, 20.625, 20.75] [ 0.0246  0.0013 -0.022 ]
   0.125
6 [20.625, 20.6875, 20.75] [ 0.0013 -0.0104 -0.022 ]
   0.0625
7 [20.625, 20.65625, 20.6875] [ 0.0013 -0.0045 -0.0104]
   0.03125
8 [20.625, 20.640625, 20.65625] [ 0.0013 -0.0016 -0.0045]
   0.015625
9 [20.625, 20.6328125, 20.640625] [ 0.0013 -0.0002 -0.0016]
   0.0078125
10 [20.625, 20.62890625, 20.6328125] [ 0.0013  0.0006 -0.0002]
   0.00390625
11 [20.62890625, 20.630859375, 20.6328125] [ 0.0006  0.0002 -0.0002]
   0.001953125
12 [20.630859375, 20.6318359375, 20.6328125] [ 1.9762e-04  1.5502e-05 -1.6661e-04]
   0.0009765625
raíz en:  20.6318359375

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

Interpolación por Diferencias finitas avanzadas

Método de la Bisección – Ejemplo con Python

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

def interpola_dfinitasAvz(xi,fi, vertabla=False,
                       precision=6, casicero = 1e-15):
    '''Interpolación de diferencias finitas
    resultado: polinomio en forma simbólica,
    redondear a cero si es menor que casicero 
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    n = len(xi)
    # revisa tamaños de paso equidistantes
    h_iguales = pasosEquidistantes(xi, casicero)
    if vertabla==True:
        np.set_printoptions(precision)
    # POLINOMIO con diferencias Finitas avanzadas
    x = sym.Symbol('x')
    polisimple = sym.S.Zero # expresión del polinomio con Sympy
    if h_iguales==True:
        tabla,titulo = dif_finitas(xi,fi,vertabla)
        h = xi[1] - xi[0]
        dfinita = tabla[0,3:]
        if vertabla==True:
            print('dfinita: ',dfinita)
            print(fi[0],'+')
        n = len(dfinita)
        polinomio = fi[0]
        for j in range(1,n,1):
            denominador = np.math.factorial(j)*(h**j)
            factor = np.around(dfinita[j-1]/denominador,precision)
            termino = 1
            for k in range(0,j,1):
                termino = termino*(x-xi[k])
            if vertabla==True:
                txt1='';txt2=''
                if n<=2 or j<=1:
                    txt1 = '('; txt2 = ')'
                print('+(',np.around(dfinita[j-1],precision),
                      '/',np.around(denominador,precision),
                      ')*',txt1,termino,txt2)
            polinomio = polinomio + termino*factor
        # simplifica multiplicando entre (x-xi)
        polisimple = polinomio.expand() 
    if vertabla==True:
        print('polinomio simplificado')
        print(polisimple)
    return(polisimple)

def dif_finitas(xi,fi, vertabla=False):
    '''Genera la tabla de diferencias finitas
    resultado en: [título,tabla]
    Tarea: verificar tamaño de vectores
    '''
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    # Tabla de Diferencias Finitas
    titulo = ['i','xi','fi']
    n = len(xi)
    ki = np.arange(0,n,1)
    tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
    tabla = np.transpose(tabla)
    # diferencias finitas vacia
    dfinita = np.zeros(shape=(n,n),dtype=float)
    tabla = np.concatenate((tabla,dfinita), axis=1)
    # Calcula tabla, inicia en columna 3
    [n,m] = np.shape(tabla)
    diagonal = n-1
    j = 3
    while (j < m):
        # Añade título para cada columna
        titulo.append('df'+str(j-2))
        # cada fila de columna
        i = 0
        while (i < diagonal):
            tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
            i = i+1
        diagonal = diagonal - 1
        j = j+1
    if vertabla==True:
        print('tabla de diferencias finitas')
        print(titulo)
        print(tabla)
    return(tabla, titulo)

def pasosEquidistantes(xi, casicero = 1e-15):
    ''' Revisa tamaños de paso h en vector xi.
    True:  h son equidistantes,
    False: h tiene tamaño de paso diferentes y dónde.
    '''
    xi = np.array(xi,dtype=float)
    n = len(xi)
    # revisa tamaños de paso equidistantes
    h_iguales = True
    if n>3: 
        dx = np.zeros(n,dtype=float)
        for i in range(0,n-1,1): # calcula hi como dx
            dx[i] = xi[i+1]-xi[i]
        for i in range(0,n-2,1): # revisa diferencias
            dx[i] = dx[i+1]-dx[i]
            if dx[i]<=casicero: # redondea cero
                dx[i]=0
            if abs(dx[i])>0:
                h_iguales=False
                print('tamaños de paso diferentes en i:',i+1,',',i+2)
        dx[n-2]=0
    return(h_iguales)

# PROGRAMA ----------------

# 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
polinomio = interpola_dfinitasAvz(xi,fi, vertabla=True)
p15 = polinomio.subs(x,15)
# literal b
derivap = polinomio.diff(x,1)
dp16 = derivap.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(polinomio)
print('p(15) = ',p15)
print('Literal b')
print(derivap)
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 ------------

def biseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
    '''
    Algoritmo de Bisección
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    error = tolera
    '''
    fa = fx(a)
    fb = fx(b)
    tramo = np.abs(b-a)
    itera = 0
    cambia = np.sign(fa)*np.sign(fb)
    if cambia<0: # existe cambio de signo f(a) vs f(b)
        if vertabla==True:
            print('método de Bisección')
            print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
            print('  ','tramo')
            np.set_printoptions(precision)
            
        while (tramo>=tolera and itera<=iteramax):
            c = (a+b)/2
            fc = fx(c)
            cambia = np.sign(fa)*np.sign(fc)
            if vertabla==True:
                print(itera,[a,c,b],np.array([fa,fc,fb]))
            if (cambia<0):
                b = c
                fb = fc
            else:
                a = c
                fa = fc
            tramo = np.abs(b-a)
            if vertabla==True:
                print('  ',tramo)
            itera = itera + 1
        respuesta = c
        # Valida respuesta
        if (itera>=iteramax):
            respuesta = np.nan

    else: 
        print(' No existe cambio de signo entre f(a) y f(b)')
        print(' f(a) =',fa,',  f(b) =',fb) 
        respuesta=np.nan
    return(respuesta)
# 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
respuesta = biseccion(fx,a,b,tolera,vertabla=True)

# SALIDA
print('raíz en: ', respuesta)