Categoría: Soluciones

Ejercicios resueltos de examen

  • s1Eva_IT2017_T2 Tanque esférico-volumen

    Ejercicio: 1Eva_IT2017_T2 Tanque esférico-volumen

    a. Planteamiento del problema

    V = \frac{\pi h^{2} (3r-h)}{3}

    Si r=1 m y V=0.75 m3,

    0.75 = \frac{\pi h^{2} (3(1)-h)}{3} 0.75 -\frac{\pi h^{2} (3(1)-h)}{3} = 0 f(h) = 0.75 -\frac{\pi h^{2} (3-h)}{3} = 0.75 -\frac{\pi}{3}(h^{2} (3)-h^{2}h) = 0.75 -\frac{\pi}{3}(3 h^{2} - h^{3}) f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3}

    b. Intervalo de búsqueda de raíz

    El tanque vacio tiene h=0 y completamente lleno h= 2r = 2(1) = 2, por lo que el intevalo tiene como extremos:

    [0,2]

    tanque esferico llenado altura h

    Verificando que exista cambio de signo en la función:

    f(0) = 0.75 -\pi (0)^{2} + \frac{\pi}{3}(0)^{3} = 0.75 f(2) = 0.75 -\pi (2)^{2} + \frac{\pi}{3}(2)^{3}= -3.4387

    y verificando al usar la gráfica dentro del intervalo:

    grafica tanque esferico llenado V vs h

    Tolerancia

    Se indica en el enunciado como 0.01 que es la medición mínima a observar con un flexómetro.

    tolera = 0.01

    cinta métrica


    c. Método de Newton-Raphson
    d. Método de Punto Fijo


    c. Método de Newton-Raphson

    El método de Newton-Raphson requiere la derivada de la función:

    x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3} f'(h) = -2\pi h + \pi h^{2}

    Tomando como punto inicial de búsqueda el extremo izquierdo del intervalo, genera una división para cero. Por lo que se mueve un poco a la derecha, algo más cercano a la raiz, viendo la gráfica por ejemplo 0.1

    x0 = 0.1

    iteración 1

    i =0

    f(0.1) = 0.75 -\pi (0.1)^{2} + \frac{\pi}{3}(0.1)^{3} =0.7196 f'(0.1) = -2\pi (0.1) + \pi (0.1)^{2} = -0.5969 x_{1} = x_0 -\frac{0.7496 }{-0.0625} = 1.3056 tramo = |x_{1}-x_{0}| = |0.1-1.3056 | = 1.2056

    iteración 2

    i =1

    f(1.3056) = 0.75 -\pi (1.3056)^{2} + \frac{\pi}{3}(1.3056)^{3} = -2.2746 f'(1.3056) = -2\pi (1.3056) + \pi (1.3056)^{2} =-2.8481 x_{1} = x_0 -\frac{-2.2746}{-2.8481} = 0.5069 tramo = |x_{2}-x_{1}|=|0.5069-1.3056|=0.7987

    iteración 3

    i =2

    f(0.5069) = 0.75 -\pi (0.5069)^{2} + \frac{\pi}{3}(0.5069)^{3} = 0.0789 f'(0.5069) = -2\pi (0.5069) + \pi (0.5069)^{2} =-2.3780 x_{1} = x_0 -\frac{0.0789}{-2.3780} = 0.5401 tramo = |x_{3}-x_{2}| =|0.5401 - 0.5069| = 0.0332

    Observe que el tramo disminuye en cada iteración , por lo que el método converge, si se siguen haciendo las operaciones se tiene que:

     [ xi, xnuevo, tramo]
    [[1.00000000e-01 1.30560920e+00 1.20560920e+00]
     [1.30560920e+00 5.06991599e-01 7.98617601e-01]
     [5.06991599e-01 5.40192334e-01 3.32007350e-02]
     [5.40192334e-01 5.39518667e-01 6.73667593e-04]]
    raiz 0.5395186666699257
    

    Instrucciones en Python

    para Método de Newton-Raphson

    # 1Eva_IT2017_T2 Tanque esférico-volumen
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    fx = lambda h: 0.75 - np.pi*(h**2)+(np.pi/3)*h**3
    dfx = lambda h: -2*np.pi*h+np.pi*(h**2)
    
    # Parámetros de método
    a = 0
    b = 2
    tolera = 0.01
    x0 = 0.1
    iteramax = 100
    
    # parametros de gráfica
    La = a
    Lb = b
    muestras = 21
    
    # PROCEDIMIENTO
    # Newton-Raphson
    tabla = []
    tramo = abs(2*tolera)
    xi = x0
    while (tramo>=tolera):
        xnuevo = xi - fx(xi)/dfx(xi)
        tramo = abs(xnuevo-xi)
        tabla.append([xi,xnuevo,tramo])
        xi = xnuevo
    tabla = np.array(tabla)
    
    # calcula para grafica
    hi = np.linspace(La,Lb,muestras)
    fi = fx(hi)
    gi = dfx(hi)
    
    # SALIDA
    print(' [ xi, xnuevo, tramo]')
    print(tabla)
    print('raiz', xnuevo)
    plt.plot(hi,fi)
    plt.plot(xi,fx(xi),'ro')
    plt.axhline(0, color='green')
    plt.xlabel('h')
    plt.ylabel('V')
    plt.title('Método Newton-Raphson')
    plt.show()
    

    Planteo con Punto Fijo


    d. Método de Punto Fijo

    Del planteamiento del problema en el literal a, se tiene que:

    0.75 = \frac{\pi h^{2} (3(1)-h)}{3}

    de donde se despeja una h:

    \frac{3(0.75)}{\pi (3(1)-h) } = h^{2} h = \sqrt{\frac{3*0.75}{\pi (3-h) }} h = \sqrt{\frac{2.25}{\pi (3-h) }}

    con lo que se obtienen las expresiones a usar en el método

    identidad = h g(h) = \sqrt{\frac{2.25}{\pi (3-h) }}

    El punto inicial de búsqueda debe encontrarse en el intervalo, se toma el mismo valor que x0 en el método de Newton-Raphson

    x0 = 0.10

    método punto fijo en tanque esférico

    Iteracion 1

    x_0= 0.10 g(0.10) = \sqrt{\frac{2.25}{\pi (3-(0.10) }}= 0.4969 tramo= |0.4969-0.10| = 0.3869

    Iteracion 2

    x_1= 0.4969 g(0.4969) = \sqrt{\frac{2.25}{\pi (3-(0.4969 ) }}= 0.5349 tramo= |0.5349- 0.4969| = 0.038

    Iteracion 3

    x_2 =0.5349 g(0.5349) = \sqrt{\frac{2.25}{\pi (3-(0.5349) }}= 0.5390 tramo= |0.5390 - 0.5349| = 0.0041

    con lo que se cumple el criterio de tolerancia, y se obtiene la raiz de:

    raiz = 0.5390

    Tabla de resultados, donde se observa que el tramo o error en cada iteración disminuye, por lo que el método converge.

     [i,xi,xi+1,tramo]
    [[1.         0.1        0.4969553  0.3969553 ]
     [2.         0.4969553  0.5349116  0.03795631]
     [3.         0.5349116  0.53901404 0.00410243]]
    raiz 0.5390140355891347
    >>> 
    

    Instrucciones en Python

    para Método de Punto-Fijo, recordamos que el método puede diverger, por lo que se añade el parámetro iteramax

    # 1Eva_IT2017_T2 Tanque esférico-volumen
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    fx = lambda h: h
    gx = lambda h: np.sqrt(2.25/(np.pi*(3-h)))
    
    a = 0.1
    b = 2
    tolera = 0.01
    iteramax = 100
    
    La = a
    Lb = b
    muestras = 21
    
    # PROCEDIMIENTO
    # Punto Fijo
    tabla = []
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    while not(tramo<=tolera or i>=iteramax):
        tabla.append([i,a,b,tramo])
        a = b
        b = gx(a)
        tramo = abs(b-a)
        i = i+1
    tabla.append([i,a,b,tramo])
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax):
        respuesta = np.nan
    tabla = np.array(tabla)
    
    # calcula para grafica
    hi = np.linspace(La,Lb,muestras)
    fi = fx(hi)
    gi = gx(hi)
    
    # SALIDA
    print(' [i, xi, xi+1, tramo]')
    print(tabla)
    print('raiz', respuesta)
    plt.plot(hi,fi, label = 'identidad', color='green')
    plt.plot(hi,gi, label = 'g(h)', color = 'orange')
    plt.plot(b,gx(b),'ro')
    plt.axhline(0, color='green')
    plt.xlabel('h')
    plt.ylabel('V')
    plt.title('Método Punto Fijo')
    plt.legend()
    plt.axhline(0, color='green')
    plt.show()
    

    Otra forma de probar la convergencia es que |g'(x)|<1 que se observa en la una gráfica adicional, lo que limita aún más el intervalo de búsqueda.

    Desarrollo en la siguiente clase.


  • s1Eva_IT2016_T3_MN Tasa interés anual

    Ejercicio: 1Eva_IT2016_T3_MN Tasa interés anual

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


    Planteamiento del problema

    La fórmula del enunciado para el problema es:

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

    que con los datos dados se convierte a:

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

    que es la forma de f(x) = 0

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

    Intervalo de búsqueda

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

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

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

    [0.01,0.5]

    Tolerancia

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

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

    tolera = 0.0001


    Método de la Bisección

    itera = 1

    a = 0.01, b = 0.5

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

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

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

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

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

    tramo = b-a = 0.255-0.01 = 0.245

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

    itera = 2

    a = 0.01, b = 0.255

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

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

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

    tramo = b-a = 0.1325-0.01 = 0.1225

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

    itera = 3

    a = 0.01, b = 0.1325

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

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

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

    tramo = b-a = 0.071-0.01 = 0.061

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

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

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

    Algoritmo en Python

    # 1Eva_IT2016_T3_MN Tasa interés anual
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    fx = lambda x: 35000*(x*(1+x)**8)/((1+x)**8 -1) -5800
    a = 0.01
    b = 0.5
    tolera = 0.0001
    
    # PROCEDIMIENTO
    cuenta = 0
    np.set_printoptions(precision=3)
    tramo = b-a
    while not(tramo<tolera):
        c = (a+b)/2
        fa = fx(a)
        fb = fx(b)
        fc = fx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if cambia < 0: a = a b = c if cambia > 0:
            a = c
            b = b
        tramo = b-a
        cuenta = cuenta+1
    
    # SALIDA
    print('       raiz en: ', c)
    print('error en tramo: ', tramo)
    print('iteraciones: ',cuenta)
    


    Método del Punto Fijo

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

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

    con lo que g(x) es:

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

    valor inicial de búsqueda

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

     i0 = 0.255

    Itera = 1

     i0 = 0.255

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

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

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

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

    Itera = 2

    i1 = g(i0 ) = 0.1387

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

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

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

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

    Itera = 3

    i2 = g(i1 ) = 0.1071

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

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

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

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

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

    Usando el algoritmo se tiene:

    iteraciones: 17
        raiz en: 0.06754389199556779
    >>>

    Algoritmo en Python

    # Algoritmo de Punto Fijo
    # x0 valor inicial de búsqueda
    # error = tolera
    
    import numpy as np
    def puntofijo(gx,antes,tolera, iteramax=50):
        itera = 1 # iteración
        nuevo = gx(antes)
        tramo = abs(nuevo-antes)
        while(tramo>=tolera and itera<=iteramax ):
            antes = nuevo
            nuevo = gx(antes)
            tramo = abs(nuevo-antes)
            itera = itera+1
        respuesta = nuevo
        # Validar respuesta
        if (itera>=iteramax ):
            respuesta = np.nan
        print('iteraciones:',itera)
        return(respuesta)
    
    # PROGRAMA ---------
    
    # INGRESO
    gx = lambda i: (58/350)*((1+i)**8-1)/((1+i)**8)
    
    a = 0       # intervalo
    b = 0.5
    
    x0 = 0.255
    tolera = 0.0001
    iteramax  = 50      # itera máximo
    
    # PROCEDIMIENTO
    respuesta = puntofijo(gx,0.255,tolera,iteramax)
    
    # SALIDA
    print('    raiz en:',respuesta)
    
  • s2Eva_IT2015_T2 EDO Deflexión de mástil

    Ejercicio: 2Eva_IT2015_T2 EDO Deflexión de mástil

    \frac{\delta ^2 y}{\delta x^2} = \frac{f}{2EI} (L-x)^2 x= 0, y = 0, \frac{\delta y}{\delta x} = 0

    Ecuación en forma estandarizada:

    y'' = \frac{f}{2EI} (L-x)^2

    Sistema de ecuaciones

    y' = z = f(x,y,z) z' = \frac{f}{2EI} (L-x)^2 = g(x,y,z)

    siendo:

    \frac{f}{2EI} =\frac{60}{2(1.25 \times 10 ^{-8}) (0.05)} = 4.8 \times 10^{-6}

    El mástil mide L=30 m y se requieren 30 intervalos, entonces h = 30/30 = 1.

    Usando muestras = tramos +1 = 31

    Se plantea una solución usando Runge Kutta de 2do Orden

    f(x,y,z) = z g(x,y,z) = (4.8 \times 10^{-6})(30-x)^2

    Desarrollo analítico

    Las iteraciones se realizan para llenar los datos de la tabla,

    tabla de resultados
    i xi yi zi
    0 0 0 0
    1 1 0.00216 0.00417
    2 2 0.00835 0.00807
    3 3 tarea ...

    iteración 1

    i = 0 ; x0= 0; y0 = 0; z0 = 0

    K_{1y} = h f(x_0,y_0, z_0)= 1(0) = 0

     

    K_{1z} = h g(x_0,y_0, z_0) = = 1(4.8 \times 10^{-6})(30-0)^2 = 0.00432 K_{2y} = h f(x_0+h,y_0+K_{1y}, z_0+K_{1z}) = = 1(0+0.00432) = 0.00432 K_{2z} = h g(x_0+h,y_0+K_{1y}, z_0+K_{1z}) = = 1(4.8 \times 10^{-6})(30-(0+1))^2 = 0.00403 y_1 = y_0 + \frac{K_{1y}+K_{2y}}{2} = = 0 + \frac{0+0.00432}{2} = 0.00216 z_1 = z_0 + \frac{K_{1z}+K_{2z}}{2} = = 0 + \frac{0.00432+0.00403}{2} = 0.00417

    iteración 2

    i = 2 ; x1= 1; y1 = 0.00216; z1 = 0.00417

    K_{1y} = h f(x_1,y_1, z_1)= 1(0.00417) = 0.00417 K_{1z} = h g(x_1,y_1, z_1) = = 1(4.8 \times 10^{-6})(30-1)^2 = 0.00403 K_{2y} = h f(x_1+h,y_1+K_{1y}, z_1+K_{1z}) = = 1(0.00417+0.00403) = 0.00821 K_{2z} = h g(x_1+h,y_1+K_{1y}, z_1+K_{1z}) = = 1(4.8 \times 10^{-6})(30-(1+1))^2 = 0.00376 y_2 = y_1 + \frac{K_{1y}+K_{2y}}{2} = = 0.00216 + \frac{0.00417+0.00821}{2} = 0.00835 z_2 = z_2 + \frac{K_{1z}+K_{2z}}{2} = = 0.00417 + \frac{0.00403+0.00376}{2} = 0.00807

    iteración 3
    i = 2; x2= 2; y2 = 0.00835; z2 = 0.00807
    tarea ...

    Para facilitar los cálculos se propone usa el algoritmo en Python, como se describe en la siguiente sección.


    Algoritmo en Python

    Al usar el algoritmo se puede comparar los resultados entre Runge-Kutta de 2do orden y de 4to Orden.

    De los resultados se presenta la siguiente gráfica

    EDO2EIT2015T2 Deflexion Mastil 01

    Observe en la gráfica la diferencia de escalas entre los ejes.

    Runge Kutta 2do Orden
     [x, y, z]
    [[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
     [  1.00000000e+00   2.16000000e-03   4.17840000e-03]
     [  2.00000000e+00   8.35680000e-03   8.07840000e-03]
     [  3.00000000e+00   1.83168000e-02   1.17096000e-02]
     [  4.00000000e+00   3.17760000e-02   1.50816000e-02]
     [  5.00000000e+00   4.84800000e-02   1.82040000e-02]
     ...
     [  2.90000000e+01   9.29856000e-01   4.32216000e-02]
     [  3.00000000e+01   9.73080000e-01   4.32240000e-02]
     [  3.10000000e+01   1.01630400e+00   4.32264000e-02]]
    
    Runge Kutta 4do Orden
     [x, y, z]
    [[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
     [  1.00000000e+00   2.11240000e-03   4.17760000e-03]
     [  2.00000000e+00   8.26240000e-03   8.07680000e-03]
     [  3.00000000e+00   1.81764000e-02   1.17072000e-02]
     [  4.00000000e+00   3.15904000e-02   1.50784000e-02]
     [  5.00000000e+00   4.82500000e-02   1.82000000e-02]
     ...
     [  2.90000000e+01   9.28800400e-01   4.31984000e-02]
     [  3.00000000e+01   9.72000000e-01   4.32000000e-02]
     [  3.10000000e+01   1.01520040e+00   4.32016000e-02]]
    >>> 
    

    Las instrucciones en Python obtener los resultados:

    # 2Eva_IT2015_T2 Deflexión de mástil
    # solución propuesta: edelros@espol.edu.ec
    import numpy as np
    
    def rungekutta2fg(fx,gx,x0,y0,z0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,3),dtype=float)
        # incluye el punto [x0,y0]
        estimado[0] = [x0,y0,z0]
        xi = x0
        yi = y0
        zi = z0
        for i in range(1,tamano,1):
            K1y = h * fx(xi,yi,zi)
            K1z = h * gx(xi,yi,zi)
            
            K2y = h * fx(xi+h, yi + K1y, zi + K1z)
            K2z = h * gx(xi+h, yi + K1y, zi + K1z)
    
            yi = yi + (K1y+K2y)/2
            zi = zi + (K1z+K2z)/2
            xi = xi + h
            
            estimado[i] = [xi,yi,zi]
        return(estimado)
    
    def rungekutta4fg(fx,gx,x0,y0,z0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,3),dtype=float)
        # incluye el punto [x0,y0]
        estimado[0] = [x0,y0,z0]
        xi = x0
        yi = y0
        zi = z0
        for i in range(1,tamano,1):
            K1y = h * fx(xi,yi,zi)
            K1z = h * gx(xi,yi,zi)
            
            K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
            K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
            
            K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
            K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)
    
            K4y = h * fx(xi+h, yi + K3y, zi + K3z)
            K4z = h * gx(xi+h, yi + K3y, zi + K3z)
    
            yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
            zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
            xi = xi + h
            
            estimado[i] = [xi,yi,zi]
        return(estimado)
    
    # INGRESO
    f = 60
    L = 30
    E = 1.25e8
    I = 0.05
    x0 = 0
    y0 = 0
    z0 = 0
    tramos = 30
    
    fx = lambda x,y,z: z
    gx = lambda x,y,z: (f/(2*E*I))*(L-x)**2
    
    # PROCEDIMIENTO
    muestras = tramos + 1
    h = L/tramos
    tabla2 = rungekutta2fg(fx,gx,x0,y0,z0,h,muestras)
    xi2 = tabla2[:,0]
    yi2 = tabla2[:,1]
    zi2 = tabla2[:,2]
    
    tabla4 = rungekutta4fg(fx,gx,x0,y0,z0,h,muestras)
    xi4 = tabla4[:,0]
    yi4 = tabla4[:,1]
    zi4 = tabla4[:,2]
    
    # SALIDA
    print('Runge Kutta 2do Orden')
    print(' [x, y, z]')
    print(tabla2)
    print('Runge Kutta 4do Orden')
    print(' [x, y, z]')
    print(tabla4)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    plt.plot(xi2,yi2, label='Runge-Kutta 2do Orden')
    plt.plot(xi4,yi4, label='Runge-Kutta 4do Orden')
    plt.title('Deflexión de mástil')
    plt.xlabel('x')
    plt.ylabel('y: deflexión')
    plt.legend()
    plt.grid()
    plt.show()
    
  • s1Eva_IT2015_T4 Lingotes metales

    Ejercicio: 1Eva_IT2015_T4 Lingotes metales

    Se plantea que cada lingote debe aportar una proporción xi al lingote nuevo a ser fundido.

    Se dispone de los compuestos de cada lingote por filas:

    compuesto = np.array([[ 20, 50, 20, 10],
                          [ 30, 40, 10, 20],
                          [ 20, 40, 10, 30],
                          [ 50, 20, 20, 10]])
    proporcion = np.array([ 27, 39.5, 14, 19.5])

    El contenido final de cada componente basado en los aportes xi de cada lingote para cada componente.

    Ejemplo para los 27 gramos de oro

    20x_1 + 30x_2+ 20x_3 + 50x_4 = 27

    se realiza el mismo procedimiento para los otros tipos de metal.

    50x_1 + 40x_2+ 40x_3 + 20x_4 = 39.5 20x_1 + 10x_2+ 10x_3 + 20x_4 = 14 10x_1 + 20x_2+ 30x_3 + 10x_4 = 19.5

    Las ecuaciones se escriben en la forma matricial Ax=B

    \begin{bmatrix} 20 && 30&& 20 &&50 \\ 50 && 40 && 40 && 20 \\ 20 && 10 && 10 && 20 \\ 10 && 20 && 30 && 10 \end{bmatrix} = \begin{bmatrix} x_1 \\ x_2 \\x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 27 \\ 39.5 \\ 14 \\ 19.5 \end{bmatrix}

    Para resolver se plantea la matriz aumentada

    \begin{bmatrix} 20 && 30&& 20 &&50 && 27\\ 50 && 40 && 40 && 20 && 39.5\\ 20 && 10 && 10 && 20 && 14 \\ 10 && 20 && 30 && 10 && 19.5 \end{bmatrix}

    se pivotea por filas la matriz:

    \begin{bmatrix} 50 && 40 && 40 && 20 && 39.5\\ 20 && 30&& 20 &&50 && 27\\ 10 && 20 && 30 && 10 && 19.5 \\ 20 && 10 && 10 && 20 && 14 \end{bmatrix}

    Para eliminar hacia adelante:

    \begin{bmatrix} 50 && 40 && 40 && 20 && 39.5 \\ 20 - 50\frac{20}{50} && 30-40\frac{20}{50} && 20-40\frac{20}{50} && 50-20\frac{20}{50} && 27-39.5\frac{20}{50}\\ 10 - 50\frac{10}{50} && 20-40\frac{10}{50} && 30-40\frac{10}{50} && 10-20\frac{10}{50} && 19.5-39.5\frac{10}{50} \\ 20 - 50\frac{20}{50} && 10-40\frac{20}{50} && 10-40\frac{20}{50} && 20-20\frac{20}{50} && 14-39.5\frac{20}{50} \end{bmatrix}

    continuando con el desarrollo:

    Elimina hacia adelante
    [[50.  40.  40.  20.  39.5]
     [ 0.  14.   4.  42.  11.2]
     [ 0.  12.  22.   6.  11.6]
     [ 0.  -6.  -6.  12.  -1.8]]
    Elimina hacia adelante
    [[ 50.  40.  40.      20.  39.5]
     [  0.  14.   4.      42.  11.2]
     [  0.   0.  18.5714 -30.   2. ]
     [  0.   0.  -4.2857  30.   3. ]]
    Elimina hacia adelante
    [[ 50.  40.  40.      20.     39.5   ]
     [  0.  14.   4.      42.     11.2   ]
     [  0.   0.  18.5714 -30.      2.    ]
     [  0.   0.   0.      23.0769  3.4615]]
    Elimina hacia adelante
    [[ 50.  40.   40.      20.     39.5  ]
     [  0.  14.    4.      42.     11.2  ]
     [  0.   0.   18.5714 -30.      2.   ]
     [  0.   0.    0.      23.0769  3.4615]]
    Elimina hacia atras
    [[ 1.    0.    0.    0.    0.25]
     [ 0.    1.    0.    0.    0.25]
     [ 0.    0.    1.   -0.    0.35]
     [ 0.    0.    0.    1.    0.15]]
    el vector solución X es:
    [[0.25]
     [0.25]
     [0.35]
     [0.15]]
    verificar que A.X = B
    [[39.5]
     [27. ]
     [19.5]
     [14. ]]
    

    Las proporciones de cada lingote a usar para el nuevo lingote que cumple con lo solicitado son:

    [0.25, 0.25, 0.35, 0.15]


    Algoritmo en python

    usado para la solución es:

    # 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
    A = np.array([[ 50., 40, 40, 20],
                  [ 20., 30, 20, 50],
                  [ 10., 20, 30, 10],
                  [ 20., 10, 10, 20]])
    B1 = np.array([ 39.5, 27, 19.5, 14])
    
    B = np.transpose([B1])
    
    # 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]))
    

    Tarea: Revisar sobre la última pregunta.

  • s1Eva_IT2015_T2 Salida cardiaca

    Ejercicio: 1Eva_IT2015_T2 Salida cardiaca

    Solución presentada como introducción al tema de interpolación y solución de sistemas de ecuaciones.
    No realiza el literal c, no se desarrolla el tema de integrales.

    Note que el desarrollo del tema permite aumentar el grado del polinomio de interpolación, lo que se afecta al tamaño del sistema de ecuaciones (matriz).

    Los valores obtenidos con la solución propuesta son:

    solución para X por Gauss-Seidel
    [[-0.42175867]
     [ 0.15610893]
     [ 0.02736763]]
    verifica que A.X = B:
    [[ -7.02831455e-05]
     [  1.50012970e+00]
     [  3.20000000e+00]]
    polinomio interpolación, con puntos:  3
    0.0273676337623498*x**2 + 0.156108926206362*x - 0.421758670607596
    >>>
    

    La gráfica para observar los datos experimentales es:
    s1EIT2015 salida cardiaca 01

    La gráfica con polinomio de interpolación de grado 2, con tres puntos:

    s1EIT2015 salida cardiaca 02

    instrucciones del problema para la solución por partes en python:

    # 1ra Evaluación I Término 2015
    # Tema 2. Flujo de sangre en corazón
    # Tarea: parte c), no se ha realizado el áre bajo curva
    #        falta calcular salida cardiaca.
    import numpy as np
    import matplotlib.pyplot as plt
    
    # Gráfica de datos experimentales:
    t = np.array([2,6,9,12,15,18])
    y = np.array([0,1.5,3.2,4.1,3.4,2.0])
    
    # SALIDA
    plt.plot(t,y)
    plt.title('datos del experimento: t vs concentración ')
    plt.show()
    
    # Sistema de ecuaciones para aproximar a polinomio grado 2
    # para grado dos usa tres puntos,
    # por ejemplo usando el punto [ 2, 0] del experimento
    # a + b*(2) + c*(2**2) =  0
    A = np.array([[1,2,2**2],
                  [1,6,6**2],
                  [1,9,9**2]])
    B = np.array([[0],
                  [1.5],
                  [3.2]])
    
    tolera = 0.0001
    X = np.zeros(len(B), dtype = float)
    
    # usando numpy para solucion de matrices
    # Xnp = np.linalg.solve(A,B)
    # print('solución para A.X=B con numpy')
    # print(Xnp)
    
    # algoritmo Gauss-Seidel
    iteramax=100
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    
    itera = 0
    while (errado>tolera or itera>iteramax):
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)
        itera = itera + 1
    # Vector en columna
    X =  np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=0
    
    Xgs = X
    
    # Metodo numérico Gauss_Seidel
    verifica = np.dot(A,Xgs)
    print('solución para X por Gauss-Seidel')
    print(Xgs)
    
    # verificar resultado
    print('verifica que A.X = B: ')
    print(verifica)
    
    # Observar interpolacion con polinomio creado
    pt = lambda t: Xgs[0,0]+ Xgs[1,0]*t + Xgs[2,0]*t**2
    
    ti = np.linspace(2,18,501)
    pti = pt(ti)
    
    plt.plot(ti,pti, label = 'interpolacion')
    plt.plot(t,y,'*', label = 'datos experimento')
    plt.title('interpolación con polinomio')
    plt.legend()
    plt.show()
    
    # polinomio en sympy
    import sympy as sp
    x = sp.Symbol('x')
    polinomio = 0
    k = len(Xgs)
    for i in range(0,k,1):
        polinomio = polinomio + Xgs[i,0]*x**i
    print('polinomio interpolación, con puntos: ', k) 
    print(polinomio)
    
  • s2Eva_IT2012_T3_MN EDO Taylor 2 Contaminación de estanque

    Ejercicio: 2Eva_IT2012_T3_MN EDO Taylor 2 Contaminación de estanque

    La ecuación a resolver con Taylor es:

    s'- \frac{26s}{200-t} - \frac{5}{2} = 0

    Para lo que se plantea usar la primera derivada:

    s'= \frac{26s}{200-t}+\frac{5}{2}

    con valores iniciales de s(0) = 0, h=0.1

    La fórmula de Taylor para tres términos es:

    s_{i+1}= s_{i}+s'_{i}h + \frac{s''_{i}}{2}h^2 + error

    Para el desarrollo se compara la solución con dos términos, tres términos y Runge Kutta.

    1. Solución con dos términos de Taylor

    Iteraciones

    i = 0, t0 = 0, s(0)=0

    s'_{0}= \frac{26s_{0}}{200-t_{0}}+\frac{5}{2} = \frac{26(0)}{200-0}+\frac{5}{2} = \frac{5}{2} s_{1}= s_{0}+s'_{0}h = 0+ \frac{5}{2}*0.1= 0.25

    t1 =  t0+h = 0+0.1 = 0.1

    i=1


    s'_{1}= \frac{26s_{1}}{200-t_{1}}+\frac{5}{2} = \frac{26(0.25)}{200-0.1}+\frac{5}{2} = 2.5325 s_{2}= s_{1}+s'_{1}h = 0.25 + (2.5325)*0.1 = 0.5032

    t2 =  t1+h = 0.1+0.1 = 0.2

    i=2,

    resolver como tarea


    2. Resolviendo con Python

    estimado
     [xi,yi Taylor,yi Runge-Kutta, diferencias]
    [[ 0.0  0.0000e+00  0.0000e+00  0.0000e+00]
     [ 0.1  2.5000e-01  2.5163e-01 -1.6258e-03]
     [ 0.2  5.0325e-01  5.0655e-01 -3.2957e-03]
     [ 0.3  7.5980e-01  7.6481e-01 -5.0106e-03]
     [ 0.4  1.0197e+00  1.0265e+00 -6.7714e-03]
     [ 0.5  1.2830e+00  1.2916e+00 -8.5792e-03]
     [ 0.6  1.5497e+00  1.5601e+00 -1.0435e-02]
     [ 0.7  1.8199e+00  1.8322e+00 -1.2339e-02]
     [ 0.8  2.0936e+00  2.1079e+00 -1.4294e-02]
     [ 0.9  2.3710e+00  2.3873e+00 -1.6299e-02]
     [ 1.0  2.6519e+00  2.6703e+00 -1.8357e-02]
     [ 1.1  2.9366e+00  2.9570e+00 -2.0467e-02]
     [ 1.2  3.2250e+00  3.2476e+00 -2.2632e-02]
     [ 1.3  3.5171e+00  3.5420e+00 -2.4853e-02]
     [ 1.4  3.8132e+00  3.8403e+00 -2.7129e-02]
     [ 1.5  4.1131e+00  4.1426e+00 -2.9464e-02]
     [ 1.6  4.4170e+00  4.4488e+00 -3.1857e-02]
     [ 1.7  4.7248e+00  4.7592e+00 -3.4310e-02]
     [ 1.8  5.0368e+00  5.0736e+00 -3.6825e-02]
     [ 1.9  5.3529e+00  5.3923e+00 -3.9402e-02]
     [ 2.0  5.6731e+00  5.7152e+00 -4.2043e-02]]
    error en rango:  0.04204310894163932

    contamina Estanque 02


    2. Algoritmo en Python

    # EDO. Método de Taylor 3 términos 
    # estima la solucion para muestras espaciadas h en eje x
    # valores iniciales x0,y0
    # entrega arreglo [[x,y]]
    import numpy as np
    
    def edo_taylor2t(d1y,x0,y0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,2),dtype=float)
        # incluye el punto [x0,y0]
        estimado[0] = [x0,y0]
        x = x0
        y = y0
        for i in range(1,tamano,1):
            y = y + h*d1y(x,y) # + ((h**2)/2)*d2y(x,y)
            x = x+h
            estimado[i] = [x,y]
        return(estimado)
    
    def rungekutta2(d1y,x0,y0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,2),dtype=float)
        # incluye el punto [x0,y0]
        estimado[0] = [x0,y0]
        xi = x0
        yi = y0
        for i in range(1,tamano,1):
            K1 = h * d1y(xi,yi)
            K2 = h * d1y(xi+h, yi + K1)
    
            yi = yi + (K1+K2)/2
            xi = xi + h
            
            estimado[i] = [xi,yi]
        return(estimado)
    
    # PROGRAMA PRUEBA
    # 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua
    
    # INGRESO.
    # d1y = y' = f, d2y = y'' = f'
    d1y = lambda x,y: 26*y/(200-x)+5/2
    x0 = 0
    y0 = 0
    h = 0.1
    muestras = 20
    
    # PROCEDIMIENTO
    puntos = edo_taylor2t(d1y,x0,y0,h,muestras)
    xi = puntos[:,0]
    yi = puntos[:,1]
    
    # Con Runge Kutta
    puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
    xiRK2 = puntosRK2[:,0]
    yiRK2 = puntosRK2[:,1]
    
    # diferencias
    diferencias = yi-yiRK2
    error = np.max(np.abs(diferencias))
    tabla = np.copy(puntos)
    tabla = np.concatenate((puntos,np.transpose([yiRK2]),
                            np.transpose([diferencias])),
                           axis = 1)
    
    # SALIDA
    np.set_printoptions(precision=4)
    print('estimado[xi,yi Taylor,yi Runge-Kutta, diferencias]')
    print(tabla)
    print('error en rango: ', error)
    
    # Gráfica
    import matplotlib.pyplot as plt
    plt.plot(xi[0],yi[0],'o',
             color='r', label ='[x0,y0]')
    plt.plot(xi[0:],yi[0:],'-',
             color='g',
             label ='y Taylor 2 términos')
    plt.plot(xiRK2[0:],yiRK2[0:],'-',
             color='blue',
             label ='y Runge-Kutta 2Orden')
    plt.axhline(y0/2)
    plt.title('EDO: Taylor 2T vs Runge=Kutta 2Orden')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()
    

    Usando Taylor con 3 términos

    estimado
     [xi,        yi,        d1yi,      d2yi      ]
    [[0.         0.         2.5        0.325     ]
     [0.1        0.251625   2.53272761 0.32958302]
     [0.2        0.50654568 2.56591685 0.33423301]
     [0.3        0.76480853 2.59957447 0.33895098]
     [0.4        1.02646073 2.63370731 0.34373796]
     [0.5        1.29155015 2.66832233 0.348595  ]
     [0.6        1.56012536 2.70342658 0.35352316]
     [0.7        1.83223563 2.73902723 0.35852351]
     [0.8        2.10793097 2.77513155 0.36359715]
     [0.9        2.38726211 2.81174694 0.36874519]
     [1.         2.67028053 2.84888087 0.37396876]
     [1.1        2.95703846 2.88654098 0.37926901]
     [1.2        3.24758891 2.92473497 0.3846471 ]
     [1.3        3.54198564 2.96347069 0.39010422]
     [1.4        3.84028323 3.00275611 0.39564157]
     [1.5        4.14253705 3.04259931 0.40126036]
     [1.6        4.44880328 3.08300849 0.40696184]
     [1.7        4.75913894 3.12399199 0.41274727]
     [1.8        5.07360187 3.16555827 0.41861793]
     [1.9        5.39225079 3.2077159  0.42457511]
     [2.         5.71514526 0.         0.        ]]

     

  • s1Eva_IT2012_T2_MN Modelo Leontief

    Ejercicio: 1Eva_IT2012_T2_MN Modelo Leontief

    Planteamiento

    X - TX = D

    A(I-T) = D

    (I-T)X = D

    para el algoritmo:

    A = I - T

    B = D


    Algoritmo en Python

    Resultados del algoritmo

    respuesta X: 
    [[158.56573701]
     [288.73225044]
     [323.87373581]]
    verificar A.X=B: 
    [[ 79.99999997]
     [139.99999998]
     [200.        ]]
    >>> itera
    8
    >>>

    Instrucciones en Python

    # Método de Gauss-Seidel
    # solución de sistemas de ecuaciones
    # por métodos iterativos
    
    import numpy as np
    
    # INGRESO
    T = np.array([[0.40, 0.03, 0.02],
                  [0.06, 0.37, 0.10],
                  [0.12, 0.15, 0.19]])
    A = np.identity(3) - T
    
    B = np.array([80.0, 140.0, 200.0],dtype=float)
    
    X0 = np.array([200.0,200.0,200.0])
    
    tolera = 0.00001
    iteramax = 100
    
    # PROCEDIMIENTO
    # Gauss-Seidel
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    #  valores iniciales
    X = np.copy(X0)
    diferencia = np.ones(n, dtype=float)
    errado = 2*tolera
    
    itera = 0
    while not(errado<=tolera or itera>iteramax):
        # por fila
        for i in range(0,n,1):
            # por columna
            suma = 0 
            for j in range(0,m,1):
                # excepto diagonal de A
                if (i!=j): 
                    suma = suma-A[i,j]*X[j]
            
            nuevo = (B[i]+suma)/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            
            X[i] = nuevo
        errado = np.max(diferencia)
        itera = itera + 1
    
    # Respuesta X en columna
    X = np.transpose([X])
    
    # revisa si NO converge
    if (itera>iteramax):
        X=0
    # revisa respuesta
    verifica = np.dot(A,X)
    
    # SALIDA
    print('respuesta X: ')
    print(X)
    print('verificar A.X=B: ')
    print(verifica)
    
  • s1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

    Ejercicio: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

    1. Desarrollo analítico

    Solución iterativa usando el método de Jacobi

    El sistema de ecuaciones

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

    cambia a su forma matricial AX=B

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

    y luego como matriz aumentada y realizar el pivoteo parcial por filas, buscando hacerla diagonal dominante:

    \begin{pmatrix} 7 & 1 & 1 &6\\ -3 & 7 &-1 & -26\\ -2 & 5 & 9 &1 \end{pmatrix}

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

    x_{i+1}=\frac{6-y_i-z_i}{7} y_{i+1}=\frac{-26+3x_i+z_i}{7} z_{i+1}=\frac{1+2x_i-5y_i}{9}

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

    El error se verifica para éste ejercicio como el mayor valor de la diferencia entre iteraciones consecutivas. Analogía al video del acople entre las aeronaves. Si una coordenada aún no es la correcta .... menor a lo tolerado, pues NO hay acople...

    itera = 0

    x_1=\frac{6-(0)-(0)}{7} = 0.8571 y_1=\frac{-26+3(0)+(0)}{7} = -3,7142 z_1=\frac{1+2(0)-5(0)}{9} =0.111 diferencia = [ 0.8571 -0 , -3,7142 -0 ,0.1111- 0 ] = [ 0.8571 , -3,7142,0.1111 ] errado = max |diferencia| = 3,7142

    error es más alto que lo tolerado

    X_1= [ 0.8571 , -3.7142 ,0.1111 ]

    Itera = 1

    x_2=\frac{6-(-3,7142) -0.1111}{7} = 1.3718 y_2=\frac{-26+3(0.8571)+0.1111}{7} = -3.3310 z_2=\frac{1+2(0.8571)-5(-3,7142)}{9} = 2.3650 diferencia = [ 1.3718-0.8571 , -3.3310 -(-3.7142) , 2.3650 - 0.1111 ] = [0.5146 , 0.3832, 2.2539 ] errado = max |diferencia| = 2.2539

    el error disminuye, pero es más alto que lo tolerado

    X_2= [ 1.3718 , -3.3310, 2.3650 ]

    Itera = 2

    x_3 = \frac{6-(-3.3310)-2.3650}{7} = 0.9951 y_3 = \frac{-26+3(1.3718)+2.3650}{7} = -2.7884 z_3 = \frac{1+2(1.3718)-5(-3.3310)}{9} = 2.2665 diferencia = [ 0.9951 - 1.3718 , -2.7884 -(-3.3310) , 2.2665- 2.3650] diferencia = [ -0.3766 , 0.5425 , -0.0985] errado = max |diferencia| = 0.5425

    el error disminuye, pero es más alto que lo tolerado

    X_3= [ 0.9951 , -2.7884 , 2.2665]

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


    2. Solución numérica usando Python

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

    Matriz aumentada
    [[ -2.   5.   9.   1.]
     [  7.   1.   1.   6.]
     [ -3.   7.  -1. -26.]]
    Pivoteo parcial:
      1 intercambiar filas:  0 y 1
      2 intercambiar filas:  1 y 2
    AB
    Iteraciones Jacobi
    itera,[X],errado
    0 [0. 0. 0.] 1.0
    1 [ 0.85714286 -3.71428571  0.11111111] 3.7142857142857144
    2 [ 1.37188209 -3.33106576  2.36507937] 2.2539682539682544
    3 [ 0.99514091 -2.78846777  2.26656589] 0.5425979915775843
    4 [ 0.93170027 -2.96400162  1.8814023 ] 0.3851635892452223
    5 [ 1.0117999  -3.04621384  1.96482318] 0.0834208883015708
    6 [ 1.01162724 -2.99996816  2.02829656] 0.06347337312830126
    7 [ 0.99595309 -2.99097453  2.00256614] 0.025730417621558477
    8 [ 0.99834406 -3.0013678   1.99408654] 0.010393267289710018
    9 [ 1.00104018 -3.00155447  2.0003919 ] 0.006305364149447268
    10 [ 1.00016608 -2.99949822  2.00109475] 0.002056248145824391
    11 [ 0.99977193 -2.99977243  1.99975814] 0.0013366043333828959
    12 [ 1.00000204 -3.0001323   1.99982289] 0.0003598675094789172
    13 [ 1.0000442  -3.00002443  2.00007395] 0.0002510632800638568
    14 [ 0.99999292 -2.99997049  2.00002339] 5.3934767063168465e-05
    respuesta de A.X=B : 
    [ 0.99999292 -2.99997049  2.00002339]
    iterado, incluyendo X0:  15
    
    

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

    jacobi 1EIIT2011T2 diagonal dominante

    Se debe controlar el número de iteraciones para verificar la convergencia con iteramax.

    NO es necesario almacenar los valores de los puntos, solo el último valor y el contador itera permite determinar si el sistema converge.

    # 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
    # MATG1052 Métodos Numéricos
    # propuesta: edelros@espol.edu.ec,
    import numpy as np
    
    def jacobi(A,B,X0, tolera, 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)
        tabla = [np.copy(X0)]
    
        itera = 0
        if vertabla==True:
            print('Iteraciones Jacobi')
            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)
            
            tabla = np.concatenate((tabla,[X]),axis = 0)
    
            itera = itera + 1
            if vertabla==True:
                print(itera, xnuevo, errado)
    
        # No converge
        if (itera>iteramax):
            X=itera
            print('iteramax superado, No converge')
        return(X,tabla)
    
    def pivoteafila(A,B,vertabla=False):
        '''
        Pivotea parcial por filas, entrega matriz aumentada AB
        Si hay ceros en diagonal es matriz singular,
        Tarea: Revisar si diagonal tiene ceros
        '''
        A = np.array(A,dtype=float)
        B = np.array(B,dtype=float)
        # Matriz aumentada
        nB = len(np.shape(B))
        if nB == 1:
            B = np.transpose([B])
        AB  = np.concatenate((A,B),axis=1)
        
        if vertabla==True:
            print('Matriz aumentada')
            print(AB)
            print('Pivoteo parcial:')
        
        # Pivoteo por filas AB
        tamano = np.shape(AB)
        n = tamano[0]
        m = tamano[1]
        
        # Para cada fila en AB
        pivoteado = 0
        for i in range(0,n-1,1):
            # columna desde diagonal i en adelante
            columna = np.abs(AB[i:,i])
            dondemax = np.argmax(columna)
            
            # dondemax no es en diagonal
            if (dondemax != 0):
                # intercambia filas
                temporal = np.copy(AB[i,:])
                AB[i,:] = AB[dondemax+i,:]
                AB[dondemax+i,:] = temporal
    
                pivoteado = pivoteado + 1
                if vertabla==True:
                    print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
        if vertabla==True:
            if pivoteado==0:
                print('  Pivoteo por filas NO requerido')
            else:
                print('AB')
        return(AB)
    
    # PROGRAMA Búsqueda de solucion
    # INGRESO --------
    
    # INGRESO
    A = [[-2, 5, 9],
         [ 7, 1, 1],
         [-3, 7,-1]]
    
    B = [1, 6,-26]
    
    X0 = [0, 0,  0]
    tolera = 0.0001
    iteramax = 100
    
    # PROCEDIMIENTO --------
    AB = pivoteafila(A,B,vertabla=True)
    # separa matriz aumentada en A y B
    [n,m] = np.shape(AB)
    A = AB[:,:m-1]
    B = AB[:,m-1]
    
    [X, puntos] = jacobi(A,B,X0,tolera,vertabla=True)
    iterado = len(puntos)
    # SALIDA --------
    print('respuesta de A.X=B : ')
    print(X)
    print('iterado, incluyendo X0: ', iterado)
    
    # GRAFICA DE ITERACIONES
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    figura = plt.figure() # Una hoja para dibujar
    grafica = figura.add_subplot(111,projection = '3d')
    
    # Puntos de la iteración
    pxi = puntos[:,0]
    pyi = puntos[:,1]
    pzi = puntos[:,2]
    grafica.plot(pxi,pyi,pzi, color = 'magenta',
                 label = 'puntos[i]')
    grafica.scatter(pxi[0],pyi[0],pzi[0], color = 'green', marker='o')
    grafica.scatter(pxi[iterado-1],pyi[iterado-1],pzi[iterado-1],
                    color = 'red', marker='o')
    
    grafica.set_title('Puntos de iteración')
    grafica.set_xlabel('eje x')
    grafica.set_ylabel('eje y')
    grafica.set_zlabel('eje z')
    grafica.legend()
    plt.show()
    

    Gráfica de planos

    Dado que el sistema de ecuaciones es de tres incógnitas, la solución se puede interpretar como la intersección de los planos formados por cada una de las ecuaciones en el espacio X,Y,Z.

    Sistema Ecuaciones Diagonal Dominante 02

    se comenta la última instrucción del algoritmo anterior, #plt.show(), y se añaden las siguientes líneas al algoritmo anterior para obtener la gráfica:

    # Tema 2. Sistema de ecuaciones 3x3
    # Concepto como interseccion de Planos
    ##import matplotlib.pyplot as plt
    ##from mpl_toolkits.mplot3d import Axes3D
    
    ax = -5     # Intervalo X
    bx = 5
    ay = ax-2   # Intervalo Y
    by = bx+2
    
    muestras = 11
    
    # PROCEDIMIENTO --------
    # Ecuaciones de planos
    z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
    z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
    z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]
    
    xi = np.linspace(ax,bx, muestras)
    yi = np.linspace(ay,by, muestras)
    Xi, Yi = np.meshgrid(xi,yi)
    
    Z0 = z0(Xi,Yi)
    Z1 = z1(Xi,Yi)
    Z2 = z2(Xi,Yi)
    
    # solución al sistema
    punto = np.linalg.solve(A,B)
    
    # SALIDA
    print('respuesta de A.X=B : ')
    print(punto)
    
    # Interseccion entre ecuación 1 y 2
    # PlanoXZ, extremo inferior de y
    Aa  = np.copy(A[0:2,[0,2]])
    Ba  = np.copy(B[0:2])
    Ba  = Ba-ay*A[0:2,1]
    pta = np.linalg.solve(Aa,Ba)
    pa  = np.array([ay])
    pxza = np.array([pta[0],ay,pta[1]])
    
    # PlanoXZ, extremo superior de y
    Bb  = np.copy(B[0:2])
    Bb  = Bb-by*A[0:2,1]
    ptb = np.linalg.solve(Aa,Bb)
    pb  = np.array([by])
    pxzb = np.array([ptb[0],by,ptb[1]])
    
    # GRAFICA de planos
    fig3D = plt.figure()
    graf3D = fig3D.add_subplot(111, projection='3d')
    
    graf3D.plot_wireframe(Xi,Yi,Z0,
                           color ='blue',
                           label='Z1')
    graf3D.plot_wireframe(Xi,Yi,Z1,
                           color ='green',
                           label='Z2')
    graf3D.plot_wireframe(Xi,Yi,Z2,
                           color ='orange',
                           label='Z3')
    # recta intersección planos 1 y 2
    graf3D.plot([pxza[0],pxzb[0]],
                [pxza[1],pxzb[1]],
                [pxza[2],pxzb[2]],
                label='Sol 1y2',
                color = 'violet',
                linewidth = 4)
    # Punto solución del sistema 3x3
    graf3D.scatter(punto[0],punto[1],punto[2],
                    color = 'red',
                    marker='o',
                    label ='punto',
                    linewidth = 6)
    
    graf3D.set_title('Sistema de ecuaciones 3x3')
    graf3D.set_xlabel('x')
    graf3D.set_ylabel('y')
    graf3D.set_zlabel('z')
    graf3D.set_xlim(ax, bx)
    graf3D.set_xlim(ay, by)
    graf3D.legend()
    graf3D.view_init(45, 45)
    # rotacion de ejes
    for angulo in range(45, 360+45, 5 ):
        graf3D.view_init(45, angulo)
        plt.draw()
        plt.pause(.001)
    plt.show()
    
    

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


    Revisión de resultados

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

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

    Ejercicio: 1Eva_IT2011_T3_MN Precios unitarios en factura, k

    Las ecuaciones basadas en las sumas de cantidad.preciounitario representan el valor pagado en cada factura.

    Siendo Xi el precio unitario de cada material:

    2x_1 + 5x_2 + 4x_3 = 35 3x_1 + 9x_2 + 8x_3 = k 5x_1 + 3x_2 + x_3 = 17

    se escriben en la forma matricial Ax=B

    \begin{bmatrix} 2 && 5 && 4 \\ 3 && 9 && 8 \\ 5 && 3 && 1 \end{bmatrix}.\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}=\begin{bmatrix} 35 \\ k \\ 17 \end{bmatrix}

    luego se escribe la matriz aumentada:

    \begin{bmatrix} 2 && 5 && 4 && 35\\ 3 && 9 && 8 && k\\ 5 && 3 && 1 && 17\end{bmatrix}

    se pivotea por filas buscando tener una matriz diagonal dominante,

    \begin{bmatrix} 5 && 3 && 1 && 17 \\ 3 && 9 && 8 && k\\ 2 && 5 && 4 && 35\\\end{bmatrix}

    Luego se usa el procedimiento de eliminación hacia adelante,

    \begin{bmatrix} 5 && 3 && 1 && 17 \\ 3-5\frac{3}{5} && 9-3\frac{3}{5} && 8-1\frac{3}{5} && k-17\frac{3}{5} \\ 2-5\frac{2}{5} && 5-3\frac{2}{5} && 4-1\frac{2}{5} && 35-17\frac{2}{5} \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && \frac{36}{5} && \frac{37}{5} && k-\frac{51}{5} \\ 0 && \frac{19}{5} && \frac{18}{5} && \frac{141}{5} \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 19 && 18 && 141 \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 19-36\frac{19}{36} && 18-37\frac{19}{36} && 141-(5k-51)\frac{19}{36} \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 0 && \frac{-55}{36} && \frac{6045-5k}{36} \end{bmatrix}

    multiplicando la última fila por 36,

    \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 0 && -55 && 6045-5k \end{bmatrix}

    con lo que se pueden obtener cada precio unitario en función de k.
    Como variante, se continua siguiendo el procedimieno de Gauss, dejando como tarea el uso de Gauss-Jordan

    -55x_3 = 6045-5k x_3 = -\frac{6045-5k}{55} 36 x_2 + 37 x_3 = 5k-51 x_2 = \frac{1}{36}(5k-51 - 37 x_3) x_2 = \frac{1}{36} \Big( 5k-51 - 37 \big(-\frac{6045-5k}{55}\big) \Big) 5x_1 + 3 x_2 +x_3 = 17 x_1 = \frac{1}{5} \Big[ 17 - 3 x_2 - x_3 \Big] x_1 = \frac{1}{5} \Big[17-3\frac{1}{36} \Big( 5k-21 - 37 \big(-\frac{6045-5k}{55}\big) \Big) - \Big( -\frac{6045-5k}{55} \Big) \Big]

    para luego simplificar las expresiones (tarea).

    En el literal c se indica que el valor de k es 65, con lo que se requiere sustituir en la solución el valor de K para encontrar los precios unitarios.

    \begin{bmatrix} 5 && 3 && 1 && 17 \\ 3 && 9 && 8 && 65\\ 2 && 5 && 4 && 35\\\end{bmatrix}

    Se encuentra que:

    el vector solución X es:

    [[-0.18181818]
     [ 5.18181818]
     [ 2.36363636]]

    Lo que muestra que debe existir un error en el planteamiento del enunciado, considerando que los precios NO deberían ser negativos como sucede con el primer precio unitario de la respuesta.

    que es lo que suponemos ser trata de corregir en el literal d, al indicar que se cambie en la matriz el valor de 5 por 5.1. Los resultados en éste caso son más coherentes con el enunciado. Todas las soluciones son positivas.

    A = np.array([[ 5.1, 3  , 1],
                  [ 3. , 9  , 8],
                  [ 2. , 5.1, 4]])
    B1 = np.array([ 17, 65, 35])
    
    el vector solución X es:
    [[0.33596838]
     [3.88669302]
     [3.62648221]]
    

    El error relativo de los precios encontrados entre las ecuaciones planteadas es:

    diferencia = [0.335968-0.181818,
                  3.886693-5.181818, 
                  3.626482-2.363636]
               = [0.154150, -1.295125, 1.262845]
    error(dolares) = max|diferencia| = 1.295125
    Por las magnitudes de los precios, el error se aprecia
    usando el error relativo referenciado 
    al mayor valor de la nueva solución
    error relativo = 1.295125/3.886693 = 0.333220
    es decir de aproximadamente 33%
    

    Para revisar otra causa del error se analiza el número de condición de la matriz:

    >>> A
    array([[5.1, 3. , 1. ],
           [3. , 9. , 8. ],
           [2. , 5.1, 4. ]])
    >>> np.linalg.cond(A)
    60.28297696795716
    

    El número de condición resulta lejano a 1, por lo que para éste problema:
    pequeños cambios en la matriz de entrada producen grandes cambios en los resultados.

    por ejemplo: un 0.1/5= 0.02 que es un 2% de variación en la entrada produce un cambio del 33% en el resultado.

    Algoritmo en Python

    Resultados del algoritmo

    Matriz aumentada
    [[ 2.  5.  4. 35.]
     [ 3.  9.  8. 65.]
     [ 5.  3.  1. 17.]]
    Pivoteo parcial:
      1 intercambiar filas:  0 y 2
    [[ 5.  3.  1. 17.]
     [ 3.  9.  8. 65.]
     [ 2.  5.  4. 35.]]
    Elimina hacia adelante:
     fila 0 pivote:  5.0
       factor:  0.6  para fila:  1
       factor:  0.4  para fila:  2
     fila 1 pivote:  7.2
       factor:  0.5277777777777778  para fila:  2
     fila 2 pivote:  -0.3055555555555558
    [[ 5.00000000e+00  3.00000000e+00  1.00000000e+00  1.70000000e+01]
     [ 0.00000000e+00  7.20000000e+00  7.40000000e+00  5.48000000e+01]
     [ 0.00000000e+00 -4.44089210e-16 -3.05555556e-01 -7.22222222e-01]]
    Elimina hacia Atras:
     fila 2 pivote:  -0.3055555555555558
       factor:  -24.2181818181818  para fila:  1
       factor:  -3.2727272727272703  para fila:  0
     fila 1 pivote:  7.1999999999999895
       factor:  0.4166666666666671  para fila:  0
     fila 0 pivote:  5.0
    [[ 1.00000000e+00  0.00000000e+00  0.00000000e+00 -1.81818182e-01]
     [ 0.00000000e+00  1.00000000e+00  0.00000000e+00  5.18181818e+00]
     [-0.00000000e+00  1.45338287e-15  1.00000000e+00  2.36363636e+00]]
    solución X: 
    [-0.18181818  5.18181818  2.36363636]
    >>>

    Instrucciones en Python

    # 1Eva_IT2011_T3_MN Precios unitarios en factura, k
    # 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)
                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([[2,5,4],
                  [3,9,8],
                  [5,3,1]])
    
    B = np.array([[35],
                  [65],
                  [17]])
    
    # 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)
    

     

  • s1Eva_IT2011_T1_MN Fondo de inversión

    Ejercicio: 1Eva_IT2011_T1_MN Fondo de inversión

    Planteamiento

    Siendo C(t), considerando A como un millón, A=1

    C(t)=Ate^{-t/3}

    Se observa la gráfica de la función que muestra el valor máximo en el intervalo entre [2, 5] , alrededor de  3

    s1EIT2011T1 Fondo Inversion 01

    Para encontrar el valor del máximo, se requiere usar la derivada respecto al tiempo y encontrar el cruce por cero.

    \frac{dC}{dt} = -t e^{-t/3} +e^{-t/3}

    Adicionalmente, para el literal b en la gráfica se incluye como meta donde C(t) disminuye a 1/4 o 0.25. Lo que estaría en el intervalo entre [10, 14].

    Instrucciones en Python para observar la función:

    # 1ra Evaluación I Término 2011
    # Tema 1. Fondo de Inversion
    import numpy as np
    import matplotlib.pyplot as plt
    
    ft = lambda t: t*np.exp(-t/3)
    
    a=0
    b=20
    tolera = 0.0001
    muestras = 101
    meta = 0.25
    # PROCEDIMIENTO
    # Observar la función entre [a,b]
    ti = np.linspace(a,b,muestras)
    fti = ft(ti)
    
    # Salida
    # Gráfica
    plt.plot(ti,fti, label='f(t)')
    plt.axhline(meta, color = 'r')
    plt.axhline(0, color = 'k')
    plt.legend()
    plt.show()
    

    literal a

    Para encontrar el máximo se puede determinar la derivada con Sympy.

    Para la derivada se usa la forma simbólica de la función, que se convierte a forma numérica lambda para evaluarla de forma más fácil y obtener la gráfica.

    # Literal a) usando derivada simbólica
    import sympy as sp
    x = sp.Symbol('x')
    fxs = x*sp.exp(-x/3)
    dfxs = fxs.diff(x,1)
    
    # convierte la expresión a lambda
    dfxn = sp.utilities.lambdify(x,dfxs,'numpy')
    dfxni = dfxn(ti)
    print('derivada de la función: ')
    print(dfxs)
    # Gráfica de la derivada.
    plt.plot(ti,dfxni, label='df(t)/dt')
    plt.axhline(0, color = 'k')
    plt.legend()
    plt.show()
    

    s1EIT2011T1 Fondo Inversion 02

    derivada de la función: 
    -x*exp(-x/3)/3 + exp(-x/3)

    Se busca la raíz con algún método, por ejemplo bisección, siendo f(x) = 0

    f(t) = \frac{dC}{dt} = -t e^{-t/3} +e^{-t/3}

    Desarrollo analítico

    itera = 0, a = 2, b= 5

    c = \frac{a+b}{2} =\frac{2+5}{2}=3.5 f(2) = -2 e^{-2/3} +e^{-2/3} =0.1711 f(3.5) = -5 e^{-3.5/3} +e^{-3.5/3} = -0.0519 f(5) = -5 e^{-5/3} +e^{-5/3} = -0.1259

    cambio de signo del lado izquierdo

    a = 2, b = 3.5

    error = | 3.5-2 | = 1.5

    itera = 1, a = 2, b = 3.5

    c = \frac{2+3.5}{2}=2.75 f(2.75) = -2.75 e^{-2.75/3} +e^{-2.75/3} = 0.0333

    cambio de signo del lado derecho

    a = 2.75, b = 3.5

    error = | 3.5-2.75 | = 0.75

    itera = 2, a = 2.75, b = 3.5

    c = \frac{2.75+3.5}{2}=3.125 f(3.125) = -3.125 e^{-3.125/3} +e^{-3.125/3} = -0.0147

    cambio de signo del lado izquierdo

    a = 2.75, b = 3.1255

    error = | 3.12555-2.75 | = 0.375

    y se continúa hasta alcanzar la tolerancia dada para el ejercicio.

    Usando el algoritmo se encuentra luego de 14 iteraciones:

    método de Bisección
    i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
       tramo
    0 [2.  3.5 5. ] [ 0.1711 -0.0519 -0.1259]
       1.5
    1 [2.   2.75 3.5 ] [ 0.1711  0.0333 -0.0519]
       0.75
    2 [2.75  3.125 3.5  ] [ 0.0333 -0.0147 -0.0519]
       0.375
    3 [2.75   2.9375 3.125 ] [ 0.0333  0.0078 -0.0147]
       0.1875
    4 [2.9375 3.0312 3.125 ] [ 0.0078 -0.0038 -0.0147]
       0.09375
    5 [2.9375 2.9844 3.0312] [ 0.0078  0.0019 -0.0038]
       0.046875
    6 [2.9844 3.0078 3.0312] [ 0.0019 -0.001  -0.0038]
       0.0234375
    7 [2.9844 2.9961 3.0078] [ 0.0019  0.0005 -0.001 ]
       0.01171875
    8 [2.9961 3.002  3.0078] [ 0.0005 -0.0002 -0.001 ]
       0.005859375
    9 [2.9961 2.999  3.002 ] [ 0.0005  0.0001 -0.0002]
       0.0029296875
    10 [2.999  3.0005 3.002 ] [ 1.1979e-04 -5.9866e-05 -2.3935e-04]
       0.00146484375
    11 [2.999  2.9998 3.0005] [ 1.1979e-04  2.9941e-05 -5.9866e-05]
       0.000732421875
    12 [2.9998 3.0001 3.0005] [ 2.9941e-05 -1.4968e-05 -5.9866e-05]
       0.0003662109375
    13 [2.9998 2.9999 3.0001] [ 2.9941e-05  7.4847e-06 -1.4968e-05]
       0.00018310546875
    14 [2.9999 3.     3.0001] [ 7.4847e-06 -3.7422e-06 -1.4968e-05]
       9.1552734375e-05
    raíz en:  3.000030517578125

    Instrucciones en Python

    # 1Eva_IT2011_T1_MN Fondo de Inversión
    # Algoritmo de Bisección
    # [a,b] se escogen de la gráfica de la función
    # error = tolera
    import numpy as np
    
    def biseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
        '''
        Algoritmo de Bisección
        Los valores de [a,b] son seleccionados
        desde la gráfica de la función
        error = tolera
        '''
        fa = fx(a)
        fb = fx(b)
        tramo = np.abs(b-a)
        itera = 0
        cambia = np.sign(fa)*np.sign(fb)
        if cambia<0: # existe cambio de signo f(a) vs f(b)
            if vertabla==True:
                print('método de Bisección')
                print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
                print('  ','tramo')
                np.set_printoptions(precision)
                
            while (tramo>=tolera and itera<=iteramax):
                c = (a+b)/2
                fc = fx(c)
                cambia = np.sign(fa)*np.sign(fc)
                if vertabla==True:
                    print(itera,np.array([a,c,b]),np.array([fa,fc,fb]))
                if (cambia<0):
                    b = c
                    fb = fc
                else:
                    a = c
                    fa = fc
                tramo = np.abs(b-a)
                if vertabla==True:
                    print('  ',tramo)
                itera = itera + 1
            respuesta = c
            # Valida respuesta
            if (itera>=iteramax):
                respuesta = np.nan
    
        else: 
            print(' No existe cambio de signo entre f(a) y f(b)')
            print(' f(a) =',fa,',  f(b) =',fb) 
            respuesta=np.nan
        return(respuesta)
    
    # INGRESO
    fx  = lambda t: -t*np.exp(-t/3)/3 + np.exp(-t/3)
    a = 2
    b = 5
    tolera = 0.0001
    
    # PROCEDIMIENTO
    respuesta = biseccion(fx,a,b,tolera,vertabla=True)
    # SALIDA
    print('raíz en: ', respuesta)
    

    Determine el monto de la inversión inicial A necesaria para que el máximo sea igual a un millón de dólares. Como el máximo se encuentra en t=3, se tiene que:

    C(t)=Ate^{-t/3} 1=A(3)e^{-3/3} A=\frac{1}{(3)e^{-1}} = 0.906093

    considerando que las unidades se encuentran en millones.


    literal b

    Para el literal b, se busca la raíz usando el método de Newton-Raphson como se indica en el enunciado.
    En la función nueva se usa el valor de A encontrado, y la meta establecida.
    Se obtiene la misa derivada del problema anterior multiplicada por A, por ser solo un factor que multiplica a la función original. El valor de meta es una constante, que se convierte en cero al derivar.

    resultado con el algoritmo:

    ['xi', 'xnuevo', 'tramo']
    [[1.0000e+01 1.0880e+01 8.7987e-01]
     [1.0880e+01 1.1056e+01 1.7580e-01]
     [1.1056e+01 1.1076e+01 2.0098e-02]
     [1.1076e+01 1.1078e+01 1.9337e-03]
     [1.1078e+01 1.1078e+01 1.8202e-04]]
    raiz en:  11.077880437153812
    con error de:  0.00018201925869298918
    >>>

    Instrucciones en Python

    # 1Eva_IT2011_T1_MN Fondo de Inversión 
    # Método de Newton-Raphson
    # Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
    
    import numpy as np
    
    # INGRESO
    fx  = lambda t: 0.906093*t*np.exp(-t/3) - 0.25
    dfx = lambda t: -t*np.exp(-t/3)/3 + np.exp(-t/3)
    
    x0 = 10
    tolera = 0.001
    
    # PROCEDIMIENTO
    tabla = []
    tramo = abs(2*tolera)
    xi = x0
    while (tramo>=tolera):
        xnuevo = xi - fx(xi)/dfx(xi)
        tramo  = abs(xnuevo-xi)
        tabla.append([xi,xnuevo,tramo])
        xi = xnuevo
    
    # convierte la lista a un arreglo.
    tabla = np.array(tabla)
    n = len(tabla)
    
    # SALIDA
    print(['xi', 'xnuevo', 'tramo'])
    np.set_printoptions(precision = 4)
    print(tabla)
    print('raiz en: ', xi)
    print('con error de: ',tramo)
    

    s1EIT2011T1 Fondo Inversion 03

    el valor de t para la meta es: 11.0779035867