Categoría: Sol_1Eva 2007-2010

  • s1Eva_IIT2009_T1 Movimiento de partícula en plano

    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
    

    s1EIIT20019T1 Mov Punto Plano 01

    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()
    
  • s1Eva_IIT2010_T1 Aproximar con polinomio

    Ejercicio: 1Eva_IIT2010_T1 Aproximar con polinomio

    Desarrollo Analítico

    Encabezado ejemplo para el desarrollo con Papel y lápiz:


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

    Opción 1. Usando el polinomio de Taylor

    Usando Unidad 1 Polinomio de Taylor, 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 +

    Se escribe la función f(x) y sus derivadas para el polinomio:

    f(x) = e^x \cos (x) +1

    primera derivada

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

    segunda derivada

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

    Punto de referencia x0 = 0, tomado como ejemplo dentro del intervalo.

    Observación: escriba las expresiones, reemplazando los valores. Un criterio de evaluación es que, 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 el polinomio de Taylor planteado:

    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

    El error referenciado a los otros puntos será:

    errado = |f(x)-p(x)| errado = \Big|\Big( e^x \cos (x) +1 \Big) - (2+ x)\Big|

    usando un punto diferente a x0=0, como x = π/2

    errado = \Big| \Big( e^{\frac{\pi}{2}}\cos \Big( \frac{\pi}{2} \Big) +1 \Big) - \Big(2+ \frac{\pi}{2}\Big) \Big| errado = | 1 - 3.5707| = 2.5707

    Tarea: calcular el error,  para x = π, verificando que pase por los puntos requeridos.



    Opción 2. Usando el polinomio de interpolación

    Siguiendo la unidad 4 para interpolación, 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 (-2/π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.

    polinomio Taylor aproximacion


    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 math
    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   = 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.

    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.

    polinomio Taylor error de aproximacion

  • s1Eva_IT2010_T2_MN Uso de televisores

    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: d2p(x)/dx2.

    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)
    
  • s1Eva_IT2010_T1_MN Demanda y producción sin,log

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

    Desarrollo Analítico

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

    demanda(t) = sin(t)

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

    oferta(t) = ln(t)

    la oferta satisface la demanda cuando ambas son iguales

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

    por lo que el tiempo t se encuentra con 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()
    
  • s1Eva_IT2009_T2 Materiales y Productos 3×4

    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 x3 que podría hacerse igual a cero. Supone que la variable libre es x3 .

    Para mantener la forma de las ecuaciones para otros valores de x3, 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 x3 = 0, lo que permite usar el modelo del problema como A.X=B .En caso de que x3 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 x3 = 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 x3, pregunte el valor al inicio, y vuelva a calcular.

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

    Observe el rango para x3, 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?

     

  • s1Eva_IT2009_T1 Demanda de un producto alcanza la producción

    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 t0 = 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 fb(t) cuando f'(t)= 0.

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

    y fb'(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

     

  • s1Eva_IIT2008_T3_MN Ganancia en inversión

    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 = [3.2 , 3.8 , 4.2 , 4.5 ]
    fi = [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}
    Pivoteo parcial por filas AB: ----
    [[91.125 20.25   4.5    1.     6.85 ]
     [74.088 17.64   4.2    1.     7.25 ]
     [54.872 14.44   3.8    1.     6.42 ]
     [32.768 10.24   3.2    1.     5.12 ]]
    Eliminación hacia adelante: ------
    i: 0 k: 1 factor: 0.8130370370370371
    i: 0 k: 2 factor: 0.6021618655692729
    i: 0 k: 3 factor: 0.35959396433470514
    [[91.125  20.25    4.5     1.      6.85  ]
     [ 0.      1.176   0.5413  0.187   1.6807]
     [ 0.      2.2462  1.0903  0.3978  2.2952]
     [ 0.      2.9582  1.5818  0.6404  2.6568]]
    i: 1 k: 2 factor: 1.9100529100529142
    i: 1 k: 3 factor: 2.5154950869236625
    [[ 9.1125e+01  2.0250e+01  4.5000e+00  1.0000e+00  6.8500e+00]
     [ 0.0000e+00  1.1760e+00  5.4133e-01  1.8696e-01  1.6807e+00]
     [ 0.0000e+00  0.0000e+00  5.6296e-02  4.0729e-02 -9.1503e-01]
     [ 0.0000e+00  0.0000e+00  2.2011e-01  1.7010e-01 -1.5710e+00]]
    i: 2 k: 3 factor: 3.9097744360903106
    [[ 9.1125e+01  2.0250e+01  4.5000e+00  1.0000e+00  6.8500e+00]
     [ 0.0000e+00  1.1760e+00  5.4133e-01  1.8696e-01  1.6807e+00]
     [ 0.0000e+00  0.0000e+00  5.6296e-02  4.0729e-02 -9.1503e-01]
     [ 0.0000e+00  0.0000e+00  0.0000e+00  1.0860e-02  2.0065e+00]]
    Eliminación hacia atras: ------
    i: 3  pivote: 0.010860484544693044
     k: 2 factor: 3.7501994302000807
     k: 1 factor: 17.214974358977564
     k: 0 factor: 92.07692307694026
    [[ 9.1125e+01  2.0250e+01  4.5000e+00  0.0000e+00 -1.7791e+02]
     [ 0.0000e+00  1.1760e+00  5.4133e-01  0.0000e+00 -3.2862e+01]
     [ 0.0000e+00  0.0000e+00  5.6296e-02  0.0000e+00 -8.4400e+00]
     [ 0.0000e+00  0.0000e+00  0.0000e+00  1.0860e-02  2.0065e+00]]
    i: 2  pivote: 0.05629629629629451
     k: 1 factor: 9.615789473684513
     k: 0 factor: 79.93421052631832
    [[ 9.1125e+01  2.0250e+01  0.0000e+00  0.0000e+00  4.9674e+02]
     [ 0.0000e+00  1.1760e+00  0.0000e+00  0.0000e+00  4.8295e+01]
     [ 0.0000e+00  0.0000e+00  5.6296e-02  0.0000e+00 -8.4400e+00]
     [ 0.0000e+00  0.0000e+00  0.0000e+00  1.0000e+00  1.8476e+02]]
    i: 1  pivote: 1.1759999999999984
     k: 0 factor: 17.219387755102066
    [[  91.125     0.        0.        0.     -334.876 ]
     [   0.        1.176     0.        0.       48.2952]
     [   0.        0.        1.        0.     -149.9209]
     [   0.        0.        0.        1.      184.7569]]
    i: 0  pivote: 91.125
    [[  91.125     0.        0.        0.     -334.876 ]
     [   0.        1.        0.        0.       41.0673]
     [   0.        0.        1.        0.     -149.9209]
     [   0.        0.        0.        1.      184.7569]]
    AB:
    [[   1.        0.        0.        0.       -3.6749]
     [   0.        1.        0.        0.       41.0673]
     [   0.        0.        1.        0.     -149.9209]
     [   0.        0.        0.        1.      184.7569]]
    Método de Gauss-Jordan
    solución X: 
    [  -3.6749   41.0673 -149.9209  184.7569]

    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
    # Sistemas de Ecuaciones A.X=B
    import numpy as np
    
    # INGRESO
    xi = [3.2 , 3.8 , 4.2 , 4.5 ]
    fi = [5.12, 6.42, 7.25, 6.85]
    
    A = [[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 = [6.85,7.25,6.42,5.12]
    
    # PROCEDIMIENTO
    np.set_printoptions(precision=4)
    casicero = 1e-15
    # Matrices como arreglo, numeros reales
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    
    # Matriz aumentada AB
    B_columna = np.transpose([B])
    AB  = np.concatenate((A,B_columna),axis=1)
    
    # Pivoteo parcial por filas
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        if (dondemax !=0): # NO en diagonal
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:]  = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal
    print('Pivoteo parcial por filas AB: ----')
    print(AB)
    
    # Eliminación hacia adelante
    print('Eliminación hacia adelante: ------')
    for i in range(0,n-1,1): # una fila
        pivote   = AB[i,i]
        adelante = i + 1
        for k in range(adelante,n,1): # diagonal adelante
            factor  = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - AB[i,:]*factor
    
            print('i:',i,'k:',k, 'factor:',factor)
        print(AB)
    
    # Eliminación hacia atras
    print('Eliminación hacia atras: ------')
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1): # una fila
        pivote   = AB[i,i]
        print('i:',i,' pivote:',pivote)
        
        atras = i - 1
        for k in range(atras,0-1,-1): # diagonal adelante
            factor  = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - AB[i,:]*factor
            for j in range(0,m,1): # casicero revisa
                if abs(AB[k,j])<casicero:
                    AB[k,j]=0
    
            print(' k:',k, 'factor:',factor)
        print(AB)
        
        AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
        
    print('AB:')
    print(AB)
    X = np.copy(AB[:,ultcolumna])
    
    # SALIDA
    print('Método de Gauss-Jordan')
    print('solución X: ')
    print(X)
    
    # Revisa polinomio
    import matplotlib.pyplot as plt
    a = np.min(xi)
    b = np.max(xi)
    muestras = 21
    
    px = lambda x: X[0]*x**3 + X[1]*x**2 +X[2]*x + X[3]
    xk = np.linspace(a,b,muestras)
    pxk = px(xk)
    
    # gráfica
    plt.plot(xk,pxk)
    plt.plot(xi,fi,'ro')
    plt.show()
    
  • s1Eva_IIT2008_T1 Distribuidores de productos

    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 X1 para A = precio en nodoA - costo de transporteA

    siguiendo el mismo procedimiento,

    precio en X1 para A: (3.1-0.2)
    precio de X1 para B: (2.8-0.3)
    precio de X1 para C: (2.7-0.4)
    precio de X1 para X2: (X2-0.1)
    precio de X1 para X3: (X3-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 X3

    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)
    
  • s1Eva_IIT2008_T1_MN Bacterias contaminantes

    Ejercicio: 1Eva_IIT2008_T1_MN Bacterias contaminantes

    a.  Planteamiento del problema

    estandarizar la fórmula a la forma f(x) = 0

    c = 70 e^{-1.5t} + 25 e^{-0.075t}

    el valor que debe tomar c = 9, por lo que la función a desarrollar se convierte en:

    9 = 70 e^{-1.5t} + 25 e^{-0.075t}

    y la que se usará en el algoritmo de búsqueda de raices es:

    70 e^{-1.5t} + 25 e^{-0.075t} -9 = 0 f(t) = 70 e^{-1.5t} + 25 e^{-0.075t} -9 = 0

    b. Intervalo de búsqueda [a,b]

    Como la variable t representa el tiempo, se inicia el análisis con cero por la izquierda o un poco mayor. a=0

    Para verificar que existe raíz se evalúa f(a) y f(b) buscando valores donde se presenten cambios de signo.

    La función depende de tiempo, por lo que a = 0, y b seleccionar un valor mayor.

    Para el caso de b,  si fuesen minutos, se pueden probar con valores considerando el tiempo que duraría un "experimento" en el laboratorio durante una clase. Por ejemplo 60 minutos, 30 minutos, etc, para obtener un punto "b" donde se garantice un cambio de signo f(a) y f(b)

    Para el ejemplo, se prueba con b = 15

    a = 0
    b = 15
    f(a) = f(0) = 70*e0 + 25 e0 -9 = 86
    f(b) = f(15) = 70*e-1.5*15. + 25 e-0.075*15 -9 = - 0.0036
    

    c. Número de iteraciones

    Calcular el número de iteraciones para llegar a la raíz con el error tolerado

    error = 0.001

    \frac{|15-0|}{2^n} = 0.001 15/0.001 = 2^n log(15/0.001) = nlog(2) n = \frac{log(15/0.001)}{log(2)} = 13.87 n = 14

    Verificando la selección usando una gráfica, usando 50 tramos entre [a,b] o 51 muestras en el intervalo.

    bacteriasLago02


    d. iteraciones

    Se encuentra la derivada f'(t) y se aplica el algoritmo Newton-Raphson con valor inicial cero.

    f(t) = 70 e^{-1.5t} + 25 e^{-0.075t} -9 = 0 f'(t) = 70(-1.5) e^{-1.5t} + 25(-0.075) e^{-0.075t}

    itera = 0 , t0=0

    f(0) = 70 e^{-1.5(0)} + 25 e^{-0.075(0)} -9 = 86 f'(0) = 70(-1.5) e^{-1.5(0)} + 25(-0.075) e^{-0.075(0)} =-106.875 t_1 = 0 - \frac{86}{-106.875} = 0.8047 errado = |0.8047- 0| = 0.8047

    itera = 1

    f(0.8047) = 70 e^{-1.5(0.8047)} + 25 e^{-0.075(0.8047)} -9 = 35.472 f'(0.8047) = 70(-1.5) e^{-1.5(0.8047)} + 25(-0.075) e^{-0.075(0.8047)} = -33.1694 x_2 = 0.8047 - \frac{35.472}{-33.1694} = 1.8741 errado = |1.8741-0.8047| = 1.0694

    itera = 2

    f(1.8741) = 70 e^{-1.5(1.8741)} + 25 e^{-0.075(1.8741)} -9 = 16.93146 f'(1.8741) = 70(-1.5) e^{-1.5(1.8741)} + 25(-0.075) e^{-0.075(1.8741)} = -7.9434 t_1 = 1.8741 - \frac{16.9314}{-7.9434} = 4.0056 errado = |4.0056 -1.8741| = 2.1315

    Se observa que el valor del error crece en cada iteración, por lo que si solo se hacen tres iteraciones, se concluye que el método No converge.
    Sin embargo con el algoritmo se observa que el método luego comienza a disminuir, por lo que se calcula hasta la iteración 7, obteniendo un error de 10-6 y una raíz en 13.62.

    El resultado final se lo encuentra al usar el algoritmo:

    método de Newton-Raphson
    i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
    0 [   0.       86.     -106.875     0.8047    0.8047]
    1 [  0.8047  35.472  -33.1694   1.8741   1.0694]
    2 [ 1.8741 16.9314 -7.9434  4.0056  2.1315]
    3 [ 4.0056  9.6848 -1.6465  9.8875  5.8819]
    4 [ 9.8875  2.9093 -0.8932 13.1445  3.257 ]
    5 [13.1445  0.3282 -0.6996 13.6136  0.4691]
    6 [ 1.3614e+01  5.7056e-03 -6.7543e-01  1.3622e+01  8.4474e-03]
    7 [ 1.3622e+01  1.8070e-06 -6.7500e-01  1.3622e+01  2.6771e-06]
    raiz en:  13.622016772385326
    

    Se usa el algoritmo en Python para encontrar el valor. El algoritmo Newton-Raphson mostrado es más pequeño que por ejemplo la bisección, pero requiere realizar un trabajo previo para encontrar la derivada de la función.

    Algoritmo en Python

    # 1Eva_IIT2008_T1_MN Bacterias contaminantes
    # Método de Newton-Raphson
    import numpy as np
    
    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)
    
    
    # PROGRAMA ---------------------------------
    fx = lambda t: 70*(np.e**(-1.5*t)) +25*(np.e**(-0.075*t))-9
    dfx = lambda t: -70*1.5*(np.e**(-1.5*t)) -25*0.075*(np.e**(-0.075*t))
    
    # INGRESO
    x0 = 0
    tolera = 0.001
    
    # PROCEDIMIENTO
    raiz = newton_raphson(fx, dfx, x0, tolera, vertabla=True)
    
    # SALIDA
    print('raiz en: ', raiz)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    a = 0
    b = 15
    muestras = 21
    
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    plt.plot(xi,fi, label='f(x)')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.grid()
    plt.legend()
    plt.show()
    

    Tarea: Para el problema, realice varios métodos y compare el número de iteraciones y el trabajo realizado al plantear el problema al implementar cada uno.

  • s1Eva_IT2008_T1 Raíz de funcion(f)

    Ejercicio: 1Eva_IT2008_T1 Raíz de funcion(f)

    Pasos a Seguir usando: BISECCION

    1. Plantear la fórmula estandarizada f(x) = 0
    2. Seleccionar el rango de análisis [a,b] donde exista cambio de signo.
    3. Calcular el número de iteraciones para llegar a la raíz con el error tolerado
    4. Calcular la raíz:
      4.1 Solución manual en papel: Realizar la tabla que muestre las iteraciones del método:
      4.2 Solución usando el algoritmo: Usar el algoritmo para encontrar la raiz.

    1. Plantear la fórmula estandarizada f(x) = 0

    \sqrt{f} . ln \Big( R\frac{\sqrt{f}}{2.51} \Big) - 1.1513 = 0

    La función depende de la variable f, por lo que por facilidad de trabajo se cambiará a x para no ser confundida con una función f(x).
    El valor de R también se reemplaza con R=5000 como indica el problema.
    El error tolerado para el problema es de 0.00001

    \sqrt{x} . ln \Big( 5000\frac{\sqrt{x}}{2.51} \Big) - 1.1513 = 0

    2. Seleccionar el rango de análisis [a,b] donde exista cambio de signo

    La función tiene un logaritmo, por lo que no será posible iniciar con cero, sin o con  valor un poco mayor a=0.01. Para el límite superior se escoge para prueba b=2. y se valida el cambio de signo en la función.

    >>> import numpy as np
    >>> fx = lambda x: np.sqrt(x)*np.log((5000/2.51)*np.sqrt(x))-1.1513
    >>> fx(0.01)
    -0.62186746547214999
    >>> fx(2)
    10.082482845673042
    

    validando el rango por cambio de signo en la función [0.01 ,2]

    3. Calcular el número de iteraciones para llegar a la raíz con el error tolerado

    error = 0.00001

    \frac{|2-0.01|}{2^n} = 0.001 1.99/0.00001 = 2^n log(1.99/0.00001) = nlog(2) n = \frac{log(1.99/0.00001)}{log(2)} = 17.6 n = 18

    4. Calcular la raíz:

    Usando el algoritmo se encuentra que la raiz está en:

    raiz:  0.0373930168152
    f(raiz) =  -2.25294254252e-06
    

    Primeras iteraciones de la tabla resultante:

     
    a,b,c, fa, fb, fc, error
    [[  1.00000000e-02   2.00000000e+00   1.00500000e+00  -6.21867465e-01
        6.46707903e+00   6.46707903e+00   1.99000000e+00]
     [  1.00000000e-02   1.00500000e+00   5.07500000e-01  -6.21867465e-01
        4.01907320e+00   4.01907320e+00   9.95000000e-01]
     [  1.00000000e-02   5.07500000e-01   2.58750000e-01  -6.21867465e-01
        2.36921961e+00   2.36921961e+00   4.97500000e-01]
     [  1.00000000e-02   2.58750000e-01   1.34375000e-01  -6.21867465e-01
        1.26563722e+00   1.26563722e+00   2.48750000e-01]
     [  1.00000000e-02   1.34375000e-01   7.21875000e-02  -6.21867465e-01
        5.36709904e-01   5.36709904e-01   1.24375000e-01]
     [  1.00000000e-02   7.21875000e-02   4.10937500e-02  -6.21867465e-01
        6.51903790e-02   6.51903790e-02   6.21875000e-02]
    ...
    

    el número de iteraciones es filas-1 de la tabla

    >>> len(tabla)
    19
    

    con lo que se comprueba los resultados.