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_IT2010_T2_MN Uso de televisores

Para la función dada:

p(x) =\frac{1}{2.5} \Big(-10 \sin \Big(\frac{12x}{7} \Big) e^{-\frac{24x}{7}} + \frac{48x}{7}e^{-\frac{8x}{7}} + 0.8 \Big)

0≤x≤4

literal a

El enunciado indica encontrar el máximo y luego el mínimo, por lo que la curva bajo análisis es la derivada de la función dp(x)/dx.

Adicionalmente, para encontrar los puntos se requiere usar el método de Newton-Raphson que corresponden a las raíces de dp(x)/dx. La función bajo análisis ahora es la derivada y para el método se usa la derivada: d²p(x)/dx².

Al usar el computador para las fórmulas, se usa la forma simbólica de la función p(x), para obtener dpx y d2px.

primera derivada: 
-3.13469387755102*x*exp(-8*x/7)
 + 2.74285714285714*exp(-8*x/7) 
+ 13.7142857142857*exp(-24*x/7)*sin(12*x/7) 
- 6.85714285714286*exp(-24*x/7)*cos(12*x/7)
segunda derivada: 
3.58250728862974*x*exp(-8*x/7) 
- 6.26938775510204*exp(-8*x/7) 
- 35.265306122449*exp(-24*x/7)*sin(12*x/7) 
+ 47.0204081632653*exp(-24*x/7)*cos(12*x/7)

La gráfica requiere la evaluación las funciones, que por simplicidad de evaluación, su formas simbólicas se convierten a su forma ‘lambda’.

Con la gráfica se verifica que la raíz de dp(x)/dx (en naranja) pasa por el máximo y mínimo de p(x) (en azul).

que se obtienen con las siguientes instrucciones en Python:

# 1ra Evaluación I Término 2010
# tema 2. encendido tv
# Tarea: aplicar el método de Newton-Raphson
# solo se muestra la función y sus derivadas 1 y 2
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# función bajo análisis en forma simbólica
x = sym.Symbol('x')
pxs = (1/2.5)*(-10*sym.sin(12*x/7)*sym.exp(-24*x/7) \
               + (48*x/7)*sym.exp(-8*x/7)+0.8)

# derivadas
dpxs = pxs.diff(x,1)
d2pxs = pxs.diff(x,2)
d2pxs = sym.expand(d2pxs)

# SALIDA
print('primera derivada: ')
print(dpxs)
print('segunda derivada: ')
print(d2pxs)

# conversion a lambda
pxn = sym.utilities.lambdify(x,pxs, 'numpy')
dpxn = sym.utilities.lambdify(x,dpxs, 'numpy')
d2pxn = sym.utilities.lambdify(x,d2pxs, 'numpy')

# observar gráfica
a = 0
b = 4
muestras = 51
tolera = 0.0001

xi = np.linspace(a,b, muestras)
pxi = pxn(xi)
dpxi = dpxn(xi)
d2pxi = d2pxn(xi)

# Gráfica
plt.plot(xi,pxi, label = 'pxi')
plt.plot(xi,dpxi, label = 'dpxi')
plt.plot(xi,d2pxi, label = 'd2pxi')
plt.axhline(0)
plt.legend()
plt.show()

literal b

Usando el método de Newton-Raphson a partir de la primera y segunda derivada según lo planteado, usando x0=0,se tiene:

f(x) = -3.13469387755102 x e^{(-8 x/7)} +2.74285714285714 e^{(-8 x/7)} + 13.7142857142857 e^{(-24 x/7)} \sin\Big(\frac{12x}{7}\Big) - 6.85714285714286 e^{(-24x/7)} \cos\Big(\frac{12x}{7}\Big) f'(x) = 3.58250728862974 x e^{(- 8 x/7)} - 6.26938775510204 e^{(- 8 x/7)} - 35.265306122449 e^{(-24x/7)} \sin \Big(\frac{12x}{7}\Big) + 47.0204081632653 e^{(-24x/7)} \cos\Big(\frac{12x}{7}\Big)

itera = 0, x0 = 0

f(0) = -3.13469387755102 (0) e^{(-8 (0)/7)} +2.74285714285714 e^{(-8 (0)/7)} + 13.7142857142857 e^{(-24 (0)/7)} \sin\Big(\frac{12(0)}{7}\Big) - 6.85714285714286 e^{(-24(0)/7)} \cos\Big(\frac{12(0)}{7}\Big) = -4.1143 f'(0) = 3.58250728862974 (0) e^{(- 8 (0)/7)} - 6.26938775510204 e^{(- 8(0)/7)} - 35.265306122449 e^{(-24(0)/7)} \sin \Big(\frac{12(0)}{7}\Big) + 47.0204081632653 e^{(-24(0)/7)} \cos\Big(\frac{12(0)}{7}\Big) = 40.751 x_1 = 0 - \frac{-4.1143}{40.751} =0.101 error = | 0.101-0 |=0.101

itera = 1

f(0.101) = -3.13469387755102 (0.101) e^{(-8 (0.101)/7)} +2.74285714285714 e^{(-8 (0.101)/7)} + 13.7142857142857 e^{(-24 (0.101)/7)} \sin\Big(\frac{12(0.101)}{7}\Big) - 6.85714285714286 e^{(-24(0.101)/7)} \cos\Big(\frac{12(0.101)}{7}\Big) = -0.9456 f'(0.101) = 3.58250728862974 (0.101) e^{(- 8 (0.101)/7)} - 6.26938775510204 e^{(- 8 (0.101)/7)} - 35.265306122449 e^{(-24(0.101)/7)} \sin \Big(\frac{12(0.101)}{7}\Big) + 47.0204081632653 e^{(-24(0.101)/7)} \cos\Big(\frac{12(0.101)}{7}\Big) =23.2054 x_2 = 0.101 - \frac{-0.9456}{23.2054} =0.1417 error = | 0.1417-0.101 |=0.0407

itera = 2

f(0.1417) = -3.13469387755102 (0.1417) e^{(-8 (0.1417)/7)} +2.74285714285714 e^{(-8 (0.1417)/7)} + 13.7142857142857 e^{(-24 (0.1417)/7)} \sin\Big(\frac{12(0.1417)}{7}\Big) - 6.85714285714286 e^{(-24(0.1417)/7)} \cos\Big(\frac{12(0.1417)}{7}\Big) = -0.11005 f'(0.1417) = 3.58250728862974 (0.1417) e^{(- 8 (0.1417)/7)} - 6.26938775510204 e^{(- 8 (0.1417)/7)} - 35.265306122449 e^{(-24(0.1417)/7)} \sin \Big(\frac{12(0.1417)}{7}\Big) + 47.0204081632653 e^{(-24(0.1417)/7)} \cos\Big(\frac{12(0.1417)}{7}\Big) = 0.17957 x_3 = 0.1417 - \frac{-0.11005}{0.17957} =0.14784 error = | 0.14784- 0.1417 |=0.0061287

se observa que el error disminuye en cada iteración, por lo que el método converge.

Se continúan las operaciones con el algoritmo obteniendo:

i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 0.     -4.1143 40.751   0.101   0.101 ]
1 [ 0.101  -0.9456 23.2054  0.1417  0.0407]
2 [ 1.4171e-01 -1.1005e-01  1.7957e+01  1.4784e-01  6.1287e-03]
3 [ 1.4784e-01 -2.1916e-03  1.7245e+01  1.4797e-01  1.2708e-04]
4 [ 1.4797e-01 -9.2531e-07  1.7231e+01  1.4797e-01  5.3701e-08]
raíz en:  0.1479664890264113

literal c

la otra raíz se encuentra con x0=1

i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 1.      0.3471 -2.207   1.1573  0.1573]
1 [ 1.1573  0.0539 -1.5338  1.1924  0.0352]
2 [ 1.1924  0.0024 -1.397   1.1941  0.0017]
3 [ 1.1941e+00  5.7065e-06 -1.3904e+00  1.1942e+00  4.1041e-06]
raíz en:  1.1941511721360376

Instrucciones en Python

# 1Eva_IT2010_T2_MN Uso de televisores

import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx,dfx,xi, tolera, iteramax=100, vertabla=False, precision=4):
    '''
    funciónx y fxderiva en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera and itera<iteramax):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1
    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera, 'No converge,se alcanzó el máximo de iteraciones')
    return(xi)

# INGRESO
fx  = lambda x: -3.13469387755102*x*np.exp(-8*x/7) \
      + 2.74285714285714*np.exp(-8*x/7) \
      + 13.7142857142857*np.exp(-24*x/7)*np.sin(12*x/7) \
      - 6.85714285714286*np.exp(-24*x/7)*np.cos(12*x/7)
dfx = lambda x: (3.58250728862974*x - 6.26938775510204 \
                 - 35.265306122449*np.exp(-16*x/7)*np.sin(12*x/7) \
                 + 47.0204081632653*np.exp(-16*x/7)*np.cos(12*x/7))*np.exp(-8*x/7)

x0 = 0.5
tolera = 0.0001

# PROCEDIMIENTO
respuesta = newton_raphson(fx,dfx,x0, tolera, vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

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

Desarrollo Analítico

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

demanda(t) = sin(t)

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

oferta(t) = ln(t)

la oferta satisface la demanda cuando ambas son iguales

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

por lo que el tiempo t se encuentra con algún método para determinar la raíz de:

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

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

Para encontrar el valor de intersección de f(t) se propone usar el método de la bisección, en el intervalo [1,3]

itera =  0

a = 1, b =3

c=\frac{1+3}{2} = 2 f(1) = sin(1) - ln(1) = 0.8415 f(3) = sin(3) - ln(3) =-0.9575 f(2) = sin(2) - ln(2) =0.2162

cambio de signo a la derecha

a =c = 2 , b = 3 tramo = |3-2| =1

itera =  1

a = 2, b =3

c=\frac{2+3}{2} = 2.5 f(2) = 0.2162 f(3) =-0.9575 f(2.5) = sin(2.5) - ln(2.5) = -0.3178

cambio de signo a la izquierda

a= 2 , b = c = 2.5 tramo = |2.5-2| = 0.5

itera =  2

a = 2, b =2.5

c=\frac{2+2.5}{2} = 2.25 f(2) = 0.2162 f(2.5) = -0.3178 f(2.25) = sin(2.25) - ln(2.25) = -0.3178

cambio de signo a la izquierda

a= 2 , b = c = 2.25 tramo = |2.25-2| = 0.25

El resto de las iteraciones se continúan con el algoritmo,

encontrando la raíz en 2.219 usando tolerancia de 0.001

---

Algoritmo en Python

Desarrollo con el método de la Bisección usando el algoritmo:

método de Bisección
i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [1, 2.0, 3] [ 0.8415 -0.9575  0.2162]
   1.0
1 [2.0, 2.5, 3] [ 0.2162 -0.9575 -0.3178]
   0.5
2 [2.0, 2.25, 2.5] [ 0.2162 -0.3178 -0.0329]
   0.25
3 [2.0, 2.125, 2.25] [ 0.2162 -0.0329  0.0965]
   0.125
4 [2.125, 2.1875, 2.25] [ 0.0965 -0.0329  0.033 ]
   0.0625
5 [2.1875, 2.21875, 2.25] [ 0.033  -0.0329  0.0004]
   0.03125
6 [2.21875, 2.234375, 2.25] [ 0.0004 -0.0329 -0.0162]
   0.015625
7 [2.21875, 2.2265625, 2.234375] [ 0.0004 -0.0162 -0.0079]
   0.0078125
8 [2.21875, 2.22265625, 2.2265625] [ 0.0004 -0.0079 -0.0037]
   0.00390625
9 [2.21875, 2.220703125, 2.22265625] [ 0.0004 -0.0037 -0.0017]
   0.001953125
10 [2.21875, 2.2197265625, 2.220703125] [ 0.0004 -0.0017 -0.0007]
   0.0009765625
raíz en:  2.2197265625

Instrucciones en Python

# 1Eva_IT2010_T1_MN Demanda y producción sin,log
import numpy as np
import matplotlib.pyplot as plt

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)

    itera = 0
    tramo = np.abs(b-a)
    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,fb,fc]))
        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
    return(respuesta)

# INGRESO
fx  = lambda t: np.sin(t) - np.log(t)

a = 1
b = 3
tolera = 0.001

# PROCEDIMIENTO
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

# GRAFICA
ad = 0
bd = 3
# intervalo oferta
# considerando que no existe oferta negativa
ao = 1
bo = 3
muestras = 21

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

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

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

fi = f(tio)

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

Ejercicio: 1Eva_IT2009_T2 Materiales y Productos 3×4

Con los datos de la tabla se plantean las ecuaciones:

	P1	P2	P3	P4
M1	0.2	0.5	0.4	0.2
M2	0.3	0	0.5	0.6
M3	0.4	0.5	0.1	0.2

La cantidad disponible de cada material es: 10, 12, 15 Kg respectivamente, los cuales deben usarse completamente.

1. Plantear el sistema de ecuaciones

0.2 x_0 + 0.5 x_1 + 0.4 x_2 + 0.2 x_3 = 10 0.3 x_0 + 0 x_1 + 0.5 x_2 + 0.6 x_3 = 12 0.4 x_0 + 0.5 x_1 + 0.1 x_2 + 0.2 x_3 = 15

Observe que hay más incógnitas que ecuaciones.

Para equiparar las ecuaciones con el número de incógnitas, podríamos suponer que uno de los productos NO se fabricará. por ejemplo el producto x₃ que podría hacerse igual a cero. Supone que la variable libre es x₃ .

Para mantener la forma de las ecuaciones para otros valores de x₃, se pasa la variable y su coeficiente a la derecha.

\begin{cases} 0.2 x_0 + 0.5 x_1 + 0.4 x_2 = 10 - 0.2 x_3 \\ 0.3 x_0 + 0 x_1 + 0.5 x_2 = 12 - 0.6 x_3 \\ 0.4 x_0 + 0.5 x_1 + 0.1 x_2 = 15 - 0.2 x_3 \end{cases}

Para analizar el ejercicio, se supondrá que el valor de x₃ = 0, lo que permite usar el modelo del problema como A.X=B .En caso de que x₃ sea diferente de cero,  el vector B modifica, y se puede proceder con el sistema de ecuaciones.

2. Convertir a la forma matricial AX = B

Siendo así, suponiendo que x₃ = 0, el ejercicio se puede desarrollar usando:

\begin{pmatrix} 0.2 && 0.5 &&0.4 &&10 \\ 0.3 && 0. && 0.5 && 12.\\ 0.4 && 0.5 && 0.1 && 15. \end{pmatrix}

que de debe pivotear por filas antes de aplicar cualquier método de solución, primero se intercambian la primera y última filas para que el primer valor de la diagonal sea el mayor en la columna.

\begin{pmatrix} 0.4 && 0.5 && 0.1 && 15. \\ 0.3 && 0. && 0.5 && 12. \\0.2 && 0.5 &&0.4 &&10 \end{pmatrix}

luego se repite el proceso a partir de la diagonal en la fila segunda, columna segunda, que al tener un valor de 5  en la tercera fila que es mayor que 0, se intercambia la fila segunda con la tercera

\begin{pmatrix} 0.4 && 0.5 && 0.1 && 15. \\0.2 && 0.5 &&0.4 &&10 \\ 0.3 && 0. && 0.5 && 12. \end{pmatrix}

Para el proceso de eliminación hacia adelante se tiene:

para pivote[0,0] = 0.4,

fila = 0 vs fila 1:
pivote = 0.4, factor = 0.2/0.4 = 0.5

\begin{pmatrix} 0.4 && 0.5 && 0.1 && 15. \\0 && 0.25 &&0.35 &&2.5 \\ 0.3 && 0. && 0.5 && 12. \end{pmatrix}

fila = 0 vs fila 2:
pivote = 0.4, factor = 0.3/0.4 = 0.75

\begin{pmatrix} 0.4 && 0.5 && 0.1 && 15. \\0 && 0.25 &&0.35 &&2.5 \\ 0 && -0.375 && 0.425 && 0.75 \end{pmatrix}

y luego para pivote [1,1]

fila = 1 vs fila 2:

pivote = 0.25, factor =-0.375/0.25 = -1.5

\begin{pmatrix} 0.4 && 0.5 && 0.1 && 15. \\0 && 0.25 &&0.35 &&2.5 \\ 0 && 0. && 0.95 && 4.5 \end{pmatrix}

Para eliminación hacia atrás los factores serán:

fila 2 pivote: 0.95
factor: 0.368421052631579 para fila: 1
factor: 0.10526315789473685 para fila: 0

\begin{pmatrix} 0.4 && 0.5 && 0 && 14.52631579 \\0 && 0.25 &&0 &&0.84210526 \\ 0 && 0. && 0.95 && 4.5 \end{pmatrix}

fila 1 pivote: 0.25
factor: 2.0 para fila: 0

\begin{pmatrix} 0.4 && 0 && 0 && 12.84210526 \\0 && 0.25 &&0 &&0.84210526 \\ 0 && 0. && 0.95 && 4.5 \end{pmatrix}

lo que da como resultado:

\begin{pmatrix} 1 && 0 && 0 && 32.10526316 \\0 && 1 &&0 &&3.36842105 \\ 0 && 0 && 1 && 4.73684211 \end{pmatrix}

que da como solución:

x= [32.10526316, 3.36842105 , 4.73684211]

---

Algoritmo en Python

Para el algoritmo se puede empezar con:

A = np.array([[0.2, 0.5, 0.4],
              [0.3, 0.0, 0.5],
              [0.4, 0.5, 0.1]])
B = np.array([10, 12, 15],dtype=float)

que luego armar el algoritmo y su ejecución, obtendría una solución semejante:

Matriz aumentada
[[ 0.2  0.5  0.4 10. ]
 [ 0.3  0.   0.5 12. ]
 [ 0.4  0.5  0.1 15. ]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
  2 intercambiar filas:  1 y 2
[[ 0.4  0.5  0.1 15. ]
 [ 0.2  0.5  0.4 10. ]
 [ 0.3  0.   0.5 12. ]]
Elimina hacia adelante:
 fila 0 pivote:  0.4
   factor:  0.5  para fila:  1
   factor:  0.7499999999999999  para fila:  2
 fila 1 pivote:  0.25
   factor:  -1.4999999999999998  para fila:  2
 fila 2 pivote:  0.95
[[ 0.4   0.5   0.1  15.  ]
 [ 0.    0.25  0.35  2.5 ]
 [ 0.    0.    0.95  4.5 ]]
Elimina hacia Atras:
 fila 2 pivote:  0.95
   factor:  0.368421052631579  para fila:  1
   factor:  0.10526315789473685  para fila:  0
 fila 1 pivote:  0.25
   factor:  2.0  para fila:  0
 fila 0 pivote:  0.4
[[ 1.          0.          0.         32.10526316]
 [ 0.          1.          0.          3.36842105]
 [ 0.          0.          1.          4.73684211]]
solución X: 
[32.10526316  3.36842105  4.73684211]
>>>

Instrucciones en Python

#1Eva_IT2009_T2 Materiales y Productos 3×4
# Método de Gauss-Jordan
# 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,lu=False,casicero = 1e-15):
    ''' Gauss elimina hacia adelante, a partir de,
    matriz aumentada y pivoteada.
    Para respuesta en forma A=L.U usar lu=True entrega[AB,L,U]
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    L = np.identity(n,dtype=float) # Inicializa L
    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,:]

                L[k,i] = factor # llena L
                
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
                print(AB)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    respuesta = AB
    if vertabla==True:
        print(AB)
    if lu==True:
        U = AB[:,:n-1]
        respuesta = [AB,L,U]
    return(respuesta)

def gauss_eliminaAtras(AB, vertabla=False, precision=5, casicero = 1e-15):
    ''' Gauss-Jordan elimina hacia atras
    Requiere la matriz triangular inferior
    Tarea: Verificar que sea triangular inferior
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    ultfila = n-1
    ultcolumna = m-1
    if vertabla==True:
        print('Elimina hacia Atras:')
        
    for i in range(ultfila,0-1,-1):
        pivote = AB[i,i]
        atras = i-1  # arriba de la fila i
        if vertabla==True:
            print(' fila',i,'pivote: ', pivote)
            
        for k in range(atras,0-1,-1):
            if (np.abs(AB[k,i])>=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')
 
        AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
    X = np.copy(AB[:,ultcolumna])
    
    if vertabla==True:
        print(AB)
    return(X)

# PROGRAMA ------------------------
# INGRESO
A = np.array([[0.2, 0.5, 0.4],
              [0.3, 0.0, 0.5],
              [0.4, 0.5, 0.1]])
B = np.array([10, 12, 15],dtype=float)

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

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_eliminaAtras(AB,vertabla=True)

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

Encontrada la solución, modifique el algoritmo para calcular B en función de x₃, pregunte el valor al inicio, y vuelva a calcular.

x3 = float(input('cantidad a producir de cuarto producto: ')

Observe el rango para x₃, por ejemplo que debe ser mayor o igual que cero, pues no hay producción negativa. De producir solamente ése un producto, el valor máximo de unidades a obtener no superaría lo posible con la cantidad de materiales disponible.

Ejemplo:

x3 = 1
B3 = [0.2,0.6,0.2]
Bnuevo = B - x3*B3

>>> Bnuevo
array([  9.8, 11.4,  14.8])

Y se vuelve a generar el sistema A.X=Bnuevo

Observación: Algunos productos se fabrican por unidades, no necesariamente por peso o volumen. ¿comentarios al respecto?

 

Ejercicio: 1Eva_IT2009_T1 Demanda de un producto alcanza la producción

Desarrollo analítico

– igualar la ecuación al valor buscado, 80

200 t e^{-0.75t} = 80

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

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

– derivada de la ecuación

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

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

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

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

iteración 0:

f(1) = 200 (1) e^{-0.75(1)} - 80 =14.4733 f'(1) = 200 e^{-0.75(1)}(1-0.75(1)) = 23.6183 t_{1} = 1 - \frac{14.4733}{23.6183} =0.3872 error = |0.3872 - 1| = 0.6128

iteración 1:

f(0.3872)= 200 (0.3872) e^{-0.75(0.3872)} - 80 = -22.0776 f'(0.3872 )= 200 e^{-0.75(0.3872)}(1-0.75(0.3872)) = 106.1511 t_{2} = 0.3872 - \frac{-22.0776}{106.1511} = 0.5952 error = |0.5952- 0.3872| = 0.208

iteración 2:

f(0.5952)=200 (0.5952) e^{-0.75(0.5952)} - 80 = -3.8242 f'(0.5952) = 200 e^{-0.75(0.5952)}(1-0.75(0.5952)) = 70.855 t_{3} = 0.5952 - \frac{-3.8242}{70.855} = 0.64916

error = |0.64916-0.5952| = 0.053972
…
tabla de iteraciones

i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 1.     14.4733 23.6183  0.3872  0.6128]
1 [  0.3872 -22.0776 106.1511   0.5952   0.208 ]
2 [ 5.9518e-01 -3.8242e+00  7.0855e+01  6.4916e-01  5.3972e-02]
3 [ 6.4916e-01 -2.1246e-01  6.3069e+01  6.5252e-01  3.3686e-03]
4 [ 6.5252e-01 -7.9031e-04  6.2600e+01  6.5254e-01  1.2625e-05]
5 [ 6.5254e-01 -1.1069e-08  6.2599e+01  6.5254e-01  1.7683e-10]
raíz en:  0.6525363029069534

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

---

Desarrollo con Python

# 1Eva_IT2009_T1 Demanda de un producto alcanza la producción

import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(fx,dfx,xi, tolera, iteramax=100, vertabla=False, precision=4):
    '''
    fx y dfx en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('método de Newton-Raphson')
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1

    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera, 'No converge,se alcanzó el máximo de iteraciones')

    return(xi)

# INGRESO
fx  = lambda t: 200*t*np.exp(-0.75*t) - 80
dfx = lambda t: 200*np.exp(-0.75*t) + 200*t*(-0.75)*np.exp(-0.75*t)

#fx  = lambda t: 200*np.exp(-0.75*t) + 200*t*(-0.75)*np.exp(-0.75*t)
#dfx = lambda t: (112.5*t - 300.0)*np.exp(-0.75*t)

x0 = 1
tolera = 0.00001

# PROCEDIMIENTO
respuesta = newton_raphson(fx,dfx,x0, tolera, vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

Para la gráfica se añade:

# grafica
a = 0
b = 2
muestras = 21

xi = np.linspace(a,b,muestras)
fi = fx(xi)
plt.plot(xi,fi)
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()

---

Literal b

Para el caso de encontrar el máximo se usaría la expresión f_b(t) cuando f'(t)= 0.

Con lo que para el algoritmo la expresión nueva f(t) es :

y f_b‘(t) es f''(t):

f_b(t) = f'(t) = 200 e^{-0.75t} + 200 t (-0.75) e^{-0.75t} f_b'(t) = f''(t) = (112.5t-300)e^{-0.75t}

se desarrolla las expresiones con Sympy :

>>> import sympy as sym
>>> t = sym.Symbol('t')
>>> f = 200*t*sym.exp(-0.75*t)-80
>>> sym.diff(f,t,1)
-150.0*t*exp(-0.75*t) + 200*exp(-0.75*t)
>>> sym.diff(f,t,2)
(112.5*t - 300.0)*exp(-0.75*t)
>>> 

Se actualiza el algoritmo con las funciones:

fx  = lambda t: 200*np.exp(-0.75*t) + 200*t*(-0.75)*np.exp(-0.75*t)
dfx = lambda t: (112.5*t - 300.0)*np.exp(-0.75*t)

y se obtiene como resultado:

método de Newton-Raphson
i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [  1.      23.6183 -88.5687   1.2667   0.2667]
1 [  1.2667   3.8674 -60.9117   1.3302   0.0635]
2 [ 1.3302e+00  1.7560e-01 -5.5445e+01  1.3333e+00  3.1671e-03]
3 [ 1.3333e+00  4.1611e-04 -5.5183e+01  1.3333e+00  7.5406e-06]
raíz en:  1.3333333332906878

 

Ejercicio: 1Eva_IIT2008_T3_MN Ganancia en inversión

Se dispone de los datos (x, f(x)), en donde x es un valor de inversión y f(x) es un valor de ganancia, ambos en miles de dólares.

xi = np.array([3.2 , 3.8 , 4.2 , 4.5 ]) 
fi = np.array([5.12, 6.42, 7.25, 6.85])

Los datos se usan junto al modelo propuesto:

f(x) = a_1 x^3 + a_2 x^2 + a_3 x + a_4

generando las expresiones del sistema de ecuaciones:

a_1 (3.2)^3 + a_2 (3.2)^2 + a_3 (3.2) + a_4 = 5.12 a_1 (3.8)^3 + a_2 (3.8)^2 + a_3 (3.8) + a_4 = 6.42 a_1 (4.2)^3 + a_2 (4.2)^2 + a_3 (4.2) + a_4 = 7.25 a_1 (4.5)^3 + a_2 (4.5)^2 + a_3 (4.5) + a_4 = 6.85

Se convierte a la forma Ax=B
\begin{bmatrix} (3.2)^3 && (3.2)^2 && (3.2) && 1 \\ (3.8)^3 && (3.8)^2 && (3.8) && 1 \\ (4.2)^3 && (4.2)^2 && (4.2) && 1 \\ (4.5)^3 && (4.5)^2 && (4.5) && 1 \end{bmatrix} . \begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ a_4 \end{bmatrix} = \begin{bmatrix} 5.12 \\ 6.42 \\ 7.25 \\6.85 \end{bmatrix}

Se crea la matriz aumentada

\begin{bmatrix} (3.2)^3 && (3.2)^2 && (3.2) && 1 && 5.12\\ (3.8)^3 && (3.8)^2 && (3.8) && 1 && 6.42 \\ (4.2)^3 && (4.2)^2 && (4.2) && 1 && 7.25 \\ (4.5)^3 && (4.5)^2 && (4.5) && 1 && 6.85 \end{bmatrix}

Se pivotea por filas:

\begin{bmatrix} (4.5)^3 && (4.5)^2 && (4.5) && 1 && 6.85 \\ (4.2)^3 && (4.2)^2 && (4.2) && 1 && 7.25 \\ (3.8)^3 && (3.8)^2 && (3.8) && 1 && 6.42 \\ (3.2)^3 && (3.2)^2 && (3.2) && 1 && 5.12 \end{bmatrix}

Como se pide un método directo, se inicia con el algoritmo de eliminación hacia adelante con factor para cada fila a partir de la diagonal.

\begin{bmatrix} (4.5)^3 && (4.5)^2 && (4.5) \\ (4.2)^3 - (4.5)^3\frac{(4.2)^3}{(4.5)^3} && (4.2)^2 - (4.5)^2\frac{(4.2)^3}{(4.5)^3} && (4.2) - (4.5)\frac{(4.2)^3}{(4.5)^3} \\ (3.8)^3 - (4.5)^3\frac{(3.8)^3}{(4.5)^3} && (3.8)^2 - (4.5)^2\frac{(3.8)^3}{(4.5)^3} && (3.8) -(4.5)\frac{(3.8)^3}{(4.5)^3} \\ (3.2)^3 - (4.5)^3\frac{(3.2)^3}{(4.5)^3} && (3.2)^2 -(4.5)^2\frac{(3.2)^3}{(4.5)^3} && (3.2) - (4.5)\frac{(3.2)^3}{(4.5)^3} \end{bmatrix}

continua a la derecha de la matriz:

\begin{bmatrix} 1 && 6.85 \\ 1 - 1\frac{(4.2)^3}{(4.5)^3} && 7.25 - 6.85\frac{(4.2)^3}{(4.5)^3} \\ 1 -1\frac{(3.8)^3}{(4.5)^3} && 6.42 - 6.85\frac{(3.8)^3}{(4.5)^3}\\1 -1\frac{(3.2)^3}{(4.5)^3} && 5.12-6.85\frac{(3.2)^3}{(4.5)^3} \end{bmatrix}

[[91.125  20.25   4.5    1.     6.85 ]
 [ 0.      1.176  0.541  0.186  1.680]
 [ 0.      2.246  1.090  0.397  2.295]
 [ 0.      2.958  1.581  0.640  2.656]]
Elimina hacia adelante
[[91.125  20.250  4.500  1.000  6.850]
 [ 0.      1.176  0.541  0.186  1.680]
 [ 0.      0.     0.056  0.040 -0.915]
 [ 0.      0.     0.220  0.170 -1.571]]
Elimina hacia adelante
[[91.125  20.250  4.500  1.000  6.850]
 [ 0.      1.176  0.541  0.186  1.680]
 [ 0.      0.000  0.056  0.040 -0.915]
 [ 0.      0.000  0.000  0.010  2.006]]
Elimina hacia adelante
[[91.125  20.250  4.500  1.000  6.850]
 [ 0.      1.176  0.541  0.186  1.680]
 [ 0.      0.     0.056  0.040 -0.915]
 [ 0.      0.     0.     0.010  2.006]]
Elimina hacia atras
[[1. 0. 0. 0.   -3.67490842]
 [0. 1. 0. 0.   41.06730769]
 [0. 0. 1. 0. -149.92086081]
 [0. 0. 0. 1.  184.75692308]]
el vector solución X es:
[[  -3.67490842]
 [  41.06730769]
 [-149.92086081]
 [ 184.75692308]]

verificar que A.X = B
[[6.85]
 [7.25]
 [6.42]
 [5.12]]

con lo que el polinomio buscado es:

f(x) = -3.67490842 x^3 + 41.06730769 x^2 + -149.92086081 x + 184.75692308

que genera la siguiente gráfica:

para encontrar la cantidad necesaria a invertir y obtener 6.0 de ganancia en f(x):

6.0 = -3.67490842 x^3 + 41.06730769 x^2 + -149.92086081 x + 184.75692308

que para usar en el algoritmo se realiza se reordena como g(x) = 0

-3.67490842 x^3 + 41.06730769 x^2 +

-149.92086081 x + 184.75692308 - 6.0 = 0
y se aplica la búsqueda de raices en el rango [a, b] que de la gráfica se estima en [3.2, 3.8]

La ejecución del algoritmo de búsqueda queda como tarea.

---

Algoritmo en python para obtener la gráfica y respuesta a la matriz:

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

# INGRESO
puntos = np.array([[3.2, 5.12],
                   [3.8, 6.42],
                   [4.2, 7.25],
                   [4.5, 6.85]])

A = np.array([[4.5**3 , 4.5**2 , 4.5 , 1],
              [4.2**3 , 4.2**2 , 4.2 , 1],
              [3.8**3 , 3.8**2 , 3.8 , 1],
              [3.2**3 , 3.2**2 , 3.2, 1]])

B = np.array([[6.85],
              [7.25],
              [6.42],
              [5.12]])

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

# Revisa polinomio
a = np.min(puntos[:,0])
b = np.max(puntos[:,0])
muestras = 51

px = lambda x: X[0]*x**3 + X[1]*x**2 +X[2]*x + X[3]
xi = np.linspace(a,b,muestras)
pxi = px(xi)

# gráfica
plt.plot(xi,pxi)
plt.plot(puntos[:,0],puntos[:,1],'ro')
plt.show()

Ejercicio: 1Eva_IIT2008_T1 Distribuidores de productos

literal a

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

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

siguiendo el mismo procedimiento,

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

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

Se continua con el mismo proceso para los siguientes nodos:

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

Para X₃

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

El sistema de ecuaciones se convierte en:

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

Para resolver se cambia a la forma Ax=B

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

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

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

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

Eliminación hacia adelante

i=0, f=0

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

Eliminación hacia adelante para i = 1, j=1

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

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

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

---

literal b

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

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

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

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

---

literal c

Se plantean las ecuaciones para el método de Jacobi a partir del sistema de ecuaciones, a partir del pivoteo por filas:

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

Si consideramos que el costo mínimo podría ser 2, el precio debería ser mayor x0 = [2,2,2]

itera =  0

x_1 = \frac{7.1 +(2) + (2)}{5} = 2.22 x_2 = \frac{ 5.5 +(2) + (2)}{4} = 2.375 x_3 = \frac{5 +(2) + (2)}{4} =2.25 errado = max|[2.22-2, 2.375-2, 2.25-2]| = max |[0.22, 0.375, 0.25]| = 0.375 X_1 = [2.22, 2.375, 2.25]|

itera = 1

x_1 = \frac{7.1 +(2.375) + (2.25)}{5} = 2.345 x_2 = \frac{ 5.5 +(2.22) + (2.25)}{4} = 2.4925 x_3 = \frac{5 +(2.22) + (2.375)}{4} =2.39875 errado = max|[2.345-2.22, 2.4925-2.375, 2.39875 - 2.25]| = 0.14874999999999972 X_2 = [2.345, 2.4925, 2.39875]

itera = 2

x_1 = \frac{7.1 +(2.4925) + (2.39875)}{5} = 2.39825 x_2 = \frac{ 5.5 +(2.345) + (2.39875)}{4} = 2.5609375 x_3 = \frac{5 +(2.345) + (2.4925)}{4} = 2.459375 errado = max|[2.39825 - 2.345, 2.5609375- 2.4925, 2.459375 - 2.39875]| = 0.06843749999999993 X_3 = [2.39825, 2.5609375, 2.459375 ]

El error disminuye entre iteraciones, por lo que el método converge.

Los datos de las iteraciones usando el algoritmo son:

Matriz aumentada
[[ 5.  -1.  -1.   7.1]
 [-1.   4.  -1.   5.5]
 [-1.  -1.   4.   5. ]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
itera,[X],errado
0 [2. 2. 2.] 1.0
1 [2.22  2.375 2.25 ] 0.375
2 [2.345   2.4925  2.39875] 0.14874999999999972
3 [2.39825   2.5609375 2.459375 ] 0.06843749999999993
4 [2.4240625  2.58940625 2.48979687] 0.030421875000000043
5 [2.43584062 2.60346484 2.50336719] 0.014058593750000181
6 [2.44136641 2.60980195 2.50982637] 0.006459179687499983
7 [2.44392566 2.61279819 2.51279209] 0.0029962402343750583
8 [2.44511806 2.61417944 2.51418096] 0.001388874511718985
9 [2.44567208 2.61482476 2.51482437] 0.0006453167724607134
10 [2.44592983 2.61512411 2.51512421] 0.00029983517456066977
11 [2.44604966 2.61526351 2.51526348] 0.0001393951034542873
12 [2.4461054  2.61532829 2.51532829] 6.480845146183967e-05
numero de condición: 2.5274158815808474
respuesta con Jacobi
[2.4461054 2.61532829 2.51532829]

Algoritmo en Python

# 1Eva_IIT2008_T1 Distribuidores de productos
# Método de Jacobi
import numpy as np

def jacobi(A,B,tolera,X0,iteramax=100, vertabla=False):
    ''' Método de Jacobi, 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 = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)
    xnuevo = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('itera,[X],errado')
        print(itera, xnuevo, errado)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        X = np.copy(xnuevo)
        itera = itera + 1
        if vertabla==True:
            print(itera, xnuevo, 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)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print('AB')
    return(AB)


# INGRESO
A = [[  5, -1, -1],
     [ -1,  4, -1],
     [ -1, -1,  4]]
B = [7.1, 5.5,5]
X0  = [2,2,2]

tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO
# numero de condicion
ncond = np.linalg.cond(A)
# Pivoteo parcial por filas
AB = pivoteafila(A,B,vertabla=True)
n = len(A)
A = AB[:,:n]
B = AB[:,n]

respuesta = jacobi(A,B,tolera,X0,vertabla=True)

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