Autor: Edison Del Rosario

  • s1Eva_IIIT2007_T1_AN Container: Refrigeradoras y Cocinas

    Ejercicio: 1Eva_IIIT2007_T1 Container: Cocinas y Refrigeradoras

    literal a

    Considerando  como  x la cantidad de refrigeradoras, y cantidad de cocinas, las ecuaciones a plantear son:

    \begin{cases} 200 x + 100 y = 1000 \\ 2 x + 1.05 y = 10.4 \end{cases}

    La forma matricial del ejercicio se convierte a:

    \begin{pmatrix} 200 & 100 \\ 2 & 1.05\end{pmatrix} \begin{pmatrix}x \\ y \end{pmatrix} = \begin{pmatrix} 1000 \\10.4 \end{pmatrix}

    la matriz aumentada es

    \begin{pmatrix} 200 & 100 & 1000\\ 2 & 1.05 & 10.4\end{pmatrix}

    Para el pivoteo parcial por filas, dado que el mayor valor de la primera columna se encuentra en la diagonal, no se requiere y la matriz aumentada se mantiene igual.

    Para el proceso de eliminación hacia adelante, se incia con el pivote=200

    factor = \frac{2}{200} = 0.01

    que se aplica a la segunda fila

             [ 200.  100.     1000  ]
    -(2/200)*[   2.    1.05     10.4]
    _________________________________
           = [   0.     0.05     0.4]

    con lo que la matriz queda:

    \begin{pmatrix} 200 & 100 & 1000\\ 0 & 0.05 & 0.4\end{pmatrix}

    Se aplica sustitución hacia atrás, desde la última fila:

    0.05 y = 0.4 y = \frac{0.4}{0.05 } = 8

    para la primera fila:

    200 x+100(8)=1000 200 x=1000-100 (8) x=\frac{1000 - 100 (8)}{200} = 1

    siendo la respuesta [1,8]

    Con el algoritmo se obtiene:

    Matriz aumentada
    [[ 200.    100.   1000.  ]
     [   2.      1.05   10.4 ]]
    Pivoteo parcial:
      Pivoteo por filas NO requerido
    Elimina hacia adelante:
      factor:  0.01  para fila:  1
    [[2.e+02 1.e+02 1.e+03]
     [0.e+00 5.e-02 4.e-01]]
    solución: 
    [1. 8.]
    >>>

    literal b

    Matriz aumentada
    [[ 200.   100.  1000. ]
     [   2.     1.1   10.4]]
    Pivoteo parcial:
      Pivoteo por filas NO requerido
    Elimina hacia adelante:
     fila 0 pivote:  200
       factor:  0.01  para fila:  1
    [[2.e+02 1.e+02 1.e+03]
     [0.e+00 1.e-01 4.e-01]]
    solución: 
    [3. 4.]
    >>> 
    

    literal c

    Observación: el pequeño cambio de volumen de la cocina no es consistente con los resultados.

    El asunto es que la forma de la refrigeradora o cocina no se adapta al volumen disponible, pues son objetos rígidos. Por lo que el sistema de ecuaciones estaría mal planteado.
    Las ecuaciones tendrían sentido si esta diseñando el mejor "tamaño" para que entren la mayor cantidad dentro de un container, sin embargo los tamaños de las refrigeradoras y cocinas se encuentran estandarizados.

    Revisamos el número de condición, que resulta ser del orden de 104, lo que confirma que el sistema está mal condicionado.

    Usando la el valor de 1.05

    >>> C = np.concatenate((A,B),axis=1)
    >>> C
    array([[  200. ,   100. ,  1000.  ],
           [    2. ,     1.05,    10.4]])
    >>> np.linalg.cond(C)
    12926.000640466344
    

    Algoritmo en Python

    # 1Eva_IIIT2007_T1 Container: Cocinas y Refrigeradoras
    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)
        if vertabla==True:
            if pivoteado==0:
                print('  Pivoteo por filas NO requerido')
            else:
                print(AB)
        return(AB)
    
    def gauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
        ''' Gauss elimina hacia adelante
        tarea: verificar términos cero
        '''
        tamano = np.shape(AB)
        n = tamano[0]
        m = tamano[1]
        if vertabla==True:
            print('Elimina hacia adelante:')
        for i in range(0,n,1):
            pivote = AB[i,i]
            adelante = i+1
            if vertabla==True:
                print(' fila',i,'pivote: ', pivote)
            for k in range(adelante,n,1):
                if (np.abs(pivote)>=casicero):
                    factor = AB[k,i]/pivote
                    AB[k,:] = AB[k,:] - factor*AB[i,:]
                    if vertabla==True:
                        print('   factor: ',factor,' para fila: ',k)
                else:
                    print('  pivote:', pivote,'en fila:',i,
                          'genera division para cero')
        if vertabla==True:
            print(AB)
        return(AB)
    
    def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
        ''' Gauss sustituye hacia atras
        '''
        tamano = np.shape(AB)
        n = tamano[0]
        m = tamano[1]
        # Sustitución hacia atras
        X = np.zeros(n,dtype=float) 
        ultfila = n-1
        ultcolumna = m-1
        for i in range(ultfila,0-1,-1):
            suma = 0
            for j in range(i+1,ultcolumna,1):
                suma = suma + AB[i,j]*X[j]
            X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
        return(X)
    
    # INGRESO
    A = [[200, 100   ],
         [  2,   1.05]]
    
    B = [1000, 10.4]
    
    # PROCEDIMIENTO
    AB = pivoteafila(A,B,vertabla=True)
    
    AB = gauss_eliminaAdelante(AB,vertabla=True)
    
    X = gauss_sustituyeAtras(AB,vertabla=True)
    
    # SALIDA
    print('solución: ')
    print(X)
    
  • s1Eva_IIT2007_T3 Interpolación inversa, encontrar xi para fi dado

    Ejercicio: 1Eva_IIT2007_T3 Interpolación inversa

    Para determinar el valor de x, usando interpolación inversa.interpola inversa

    f(0.50) = 1.648
    f(0.65) = 1.915
    f( x  ) = 2.117
    f(0.80) = 2.225
    f(0.95) = 2.5857

    Para el algoritmo se intercambian las variables previo a usarlo.

    fi = [0.50 , 0.65 , 0.80,  0.95   ]
    xi = [1.648, 1.915, 2.225, 2.5857 ]
    

    Luego se evalúa en el punto buscado, en éste caso: fi=2.117, obteniendo que x es: 0.750321134121361

    Para obtener el polinomio se usa el método de Lagrange:

    término 1

    L_{0} (x) = \frac{(x-1.915)(x-2.225)(x-2.5857)}{(1.648-1.915)(1.648-2.225)(1.648-2.5857)}

    término 2

    L_{1} (x) = \frac{(x-1.648)(x-2.225)(x-2.5857)}{(1.915-1.648)(1.915-2.225)(1.915-2.5857)}

    término 3

    L_{2} (x) = \frac{(x-1.648)(x-1.915)(x-2.5857)}{(2.225-1.648)(2.225-1.915)(2.225-2.5857)}

    término 4

    L_{3} (x) = \frac{(x-1.648)(x-1.915)(x-2.225)}{(2.5857-1.648)(2.5857-1.915)(2.5857-2.225)}

    se construye el polinomio usando la fórmula para fn(x) para cada valor fi,

    p_3(x) = 0.5 L_{0} (x) + 0.65 L_{1} (x) + 0.8 L_{2} (x) + 0.95 L_{3} (x) p_3(x) = 0.5 \frac{(x-1.915)(x-2.225)(x-2.5857)}{-0.14446112} + +0.65 \frac{(x-1.648)(x-2.225)(x-2.5857)}{0.05551384 } + + 0.8 \frac{(x-1.648)(x-1.915)(x-2.5857)}{-0.06451841} + +0.95 \frac{(x-1.648)(x-1.915)(x-2.225)}{0.22684978} p_3(x) = 0.03588 x^3 - 0.34275 x^2 + 1.44073 x - 1.10404

    A partir del resultado del algoritmo se puede evaluar p(2.117)

        valores de fi:  [0.5  0.65 0.8  0.95]
    divisores en L(i):  [-0.14446112  0.05551384 -0.06451841  0.22684978]
    
    Polinomio de Lagrange, expresiones
    -3.46113878334255*(x - 2.5857)*(x - 2.225)*(x - 1.915) 
    + 11.7087921085767*(x - 2.5857)*(x - 2.225)*(x - 1.648) 
    - 12.3995618056856*(x - 2.5857)*(x - 1.915)*(x - 1.648) 
    + 4.18779332775953*(x - 2.225)*(x - 1.915)*(x - 1.648)
    
    Polinomio de Lagrange: 
    0.0358848473081546*x**3 - 0.342756582990933*x**2 
    + 1.44073214117569*x - 1.10404634485234
    >>> polisimple.subs(x,2.117)
    0.750321134121178
    >>> polisimple
    0.0358848473081546*x**3 - 0.342756582990933*x**2 
    + 1.44073214117569*x - 1.10404634485234
    >>> 
    

    Algoritmo en Python

    # 1Eva_IIT2007_T3 Interpolación inversa
    # Interpolacion de Lagrange
    # divisoresL solo para mostrar valores
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    # INGRESO , Datos de prueba
    fi = [0.50 , 0.65 , 0.80,  0.95   ]
    xi = [1.648, 1.915, 2.225, 2.5857 ]
    
    # PROCEDIMIENTO
    xi = np.array(xi,dtype=float)
    fi = np.array(fi,dtype=float)
    # Polinomio de Lagrange
    n = len(xi)
    x = sym.Symbol('x')
    polinomio = 0
    divisorL = np.zeros(n, dtype = float)
    for i in range(0,n,1):
        
        # Termino de Lagrange
        numerador = 1
        denominador = 1
        for j  in range(0,n,1):
            if (j!=i):
                numerador = numerador*(x-xi[j])
                denominador = denominador*(xi[i]-xi[j])
        terminoLi = numerador/denominador
    
        polinomio = polinomio + terminoLi*fi[i]
        divisorL[i] = denominador
    
    # simplifica el polinomio
    polisimple = polinomio.expand()
    
    # Salida
    print('    valores de fi: ',fi)
    print('divisores en L(i): ',divisorL)
    print()
    print('Polinomio de Lagrange, expresiones')
    print(polinomio)
    print()
    print('Polinomio de Lagrange: ')
    print(polisimple)

    para revisar con la gráfica, se añaden las líneas,

    # GRAFICA
    polinomio = sym.lambdify(x,polisimple)
    y = np.linspace(1,3,100)
    pyi= polinomio(y)
    plt.plot(pyi,y,label='p(x)')
    plt.plot(fi,xi,'o')
    plt.axhline(2.117,linestyle='dashed',
                label='f(x)=2.117', color='green')
    plt.legend()
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.show()
    
  • s1Eva_IIT2007_T2 Aplicar Gauss-Seidel 6x6

    Ejercicio: 1Eva_IIT2007_T2 Aplicar Gauss-Seidel 6x6

    1. Desarrollo Analítico

    - Verificar que la matriz es diagonal dominante

    No es necesario realizar el pivoteo por filas, ya la matriz tiene la diagonal dominante.

    A = [[7.63, 0.30, 0.15,  0.50, 0.34, 0.84],
         [0.38, 6.40, 0.70,  0.90, 0.29, 0.57],
         [0.83, 0.19, 8.33,  0.82, 0.34, 0.37],
         [0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
         [0.71, 0.30, 0.85,  0.82, 5.95, 0.55],
         [0.43, 0.54, 0.59,  0.66, 0.31, 9.25]]
    
    B = [ -9.44, 25.27, -48.01, 19.76, -23.63, 62.59]
    

    - revisar el número de condición

    cond(A) = || A||.||A-1||

    El número de condición no es "muy alto",  los valores de la diagonal son los mayores en toda la fila, por lo que el sistema converge.

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

    - realizar iteraciones con las expresiones:

    x_0 = \Big(-9.44 -0.3 x_1 -0.15x_2 -0.50 x_3 -0.34 x_4 -0.84 x_5\Big)\frac{1}{7.63} x_1 = \Big( 25.27 -0.38 x_0 -0.70 x_2 -0.90 x_3 -0.29 x_4 -0.57 x_5 \Big) \frac{1}{6.40} x_2 = \Big(-48.01 -0.83x_0 -0.19 x_1 -0.82x_3 -0.34x_4 -0.37x_5 \Big)\frac{1}{8.33} x_3 = \Big(19.76 -0.50 x_0 -0.68 x_1 -0.86 x_2 -0.53 x_4 -0.70x_5 \Big)\frac{1}{10.21} x_4 = \Big( -23.63 - 0.71 x_0 -0.30 x_1 -0.85 x_2 -0.82 x_3 -0.55x_5 \Big)\frac{1}{5.95} x_5 = \Big( 62.59 - 0.43x_0 -0.54x_1 - 0.59 x_2 -0.66 x_3 -0.31 x_4 \Big)\frac{1}{9.25}

    Dado que no se establece en el enunciado el vector inicial, se usará el vector cero. La tolerancia requerida es 10-5
    X0 = [ 0. 0. 0. 0. 0. 0.] = [x0,x1,x2,x3,x4,x5]

    iteración 0

    x_0 = \Big(-9.44 -0.3 (0) -0.15(0) -0.50 (0) -0.34(0) -0.84 (0)\Big)\frac{1}{7.63} = -1.23722 x_1 = \Big( 25.27 -0.38 (-1.23722) -0.70 (0) -0.90 (0) -0.29 (0) -0.57 (0) \Big) \frac{1}{6.40} = 4.0219 x_2 = \Big(-48.01 -0.83 (-1.23722) -0.19 (4.0219) -0.82(0) -0.34(0) -0.37(0) \Big)\frac{1}{8.33} = -5.73196 x_3 = \Big( 19.76 -0.50 (-1.23722) -0.68 (4.0219) -0.86 (-5.73196) -0.53 (0) -0.70(0) \Big)\frac{1}{10.21} = 2.21089 x_4 = \Big( -23.63 - 0.71 (-1.23722) -0.30 (4.0219) -0.85 (-5.73196) -0.82 (2.21089) -0.55(0) \Big)\frac{1}{5.95} = -3.51242 x_5 = \Big( 62.59 - 0.43(-1.23722) -0.54(4.0219) - 0.59 (-5.73196) -0.66 (2.21089) -0.31 (-3.51242) \Big)\frac{1}{9.25} = 6.91478 X_1 = [-1.23722, 4.0219, -5.73196, 2.21089, -3.51242, 6.91478 ] diferencia = X1 - [0,0,0,0,0,0] errado = max(|diferencia|) = 6.91478

    iteración 1

    x_0 = \Big(-9.44 -0.3 (4.0219) -0.15(4.0219) -0.50 (2.21089) -0.34 (-3.51242) -0.84(6.91478) \Big)\frac{1}{7.63} = -2.0323 x_1 = \Big( 25.27 -0.38 (-2.0323) -0.70 ( -5.73196) -0.90 (2.21089) -0.29 (-3.51242) -0.57 (6.91478) \Big) \frac{1}{6.40} = 3.92844 x_2 = \Big(-48.01 -0.83(-2.0323) -0.19 (3.92844) -0.82(2.21089) -0.34(-3.51242) -0.37(6.91478) \Big)\frac{1}{8.33} = -6.03203 x_3 = \Big(19.76 -0.50 (-2.0323) -0.68 (3.92844) -0.86 (-6.03203) -0.53 (-3.51242) -0.70(6.91478) \Big)\frac{1}{10.21} = 1.98958 x_4 = \Big( -23.63 - 0.71 (-2.0323) -0.30 (3.92844) -0.85 (-6.03203) -0.82 (1.98958) -0.55(6.91478) \Big)\frac{1}{5.95} = -3.97865 x_5 = \Big( 62.59 - 0.43(-2.0323) -0.54(3.92844) - 0.59 (-6.03203) -0.66 (1.98958) -0.31 (-3.97865) \Big)\frac{1}{9.25} = 7.00775 X_2 = [-2.0323, 3.92844, -6.03203, 1.98958, -3.97865, 7.00775] diferencia = X2 - [-1.23722, 4.0219, -5.73196, 2.21089, -3.51242, 6.91478 ] errado = max(|diferencia|) = 0.79507

    iteración 2

    Desarrollar como tarea.


    2. Algoritmo en Python

    La tabla de aproximaciones sucesivas para el vector X es:

    Pivoteo por filas NO requerido
    Iteraciones Gauss-Seidel
    itera,[X],[errado]
    0 [0. 0. 0. 0. 0. 0.] [1.]
    1 [-1.23722  4.0219  -5.73196  2.21089 -3.51242  6.91478] [6.91478]
    2 [-2.0323   3.92844 -6.03203  1.98958 -3.97865  7.00775] [0.79507]
    3 [-1.99768  4.00317 -6.00049  1.99808 -4.00082  6.9999 ] [0.07473]
    4 [-1.99994  4.00037 -5.99979  2.      -4.00005  6.99996] [0.00281]
    5 [-2.00001  3.99998 -6.       2.00001 -4.       7.     ] [0.00038]
    6 [-2.  4. -6.  2. -4.  7.] [1.60155e-05]
    7 [-2.  4. -6.  2. -4.  7.] [1.74554e-06]
    Matriz aumentada y pivoteada:
    [[  7.63   0.3    0.15   0.5    0.34   0.84  -9.44]
     [  0.38   6.4    0.7    0.9    0.29   0.57  25.27]
     [  0.83   0.19   8.33   0.82   0.34   0.37 -48.01]
     [  0.5    0.68   0.86  10.21   0.53   0.7   19.76]
     [  0.71   0.3    0.85   0.82   5.95   0.55 -23.63]
     [  0.43   0.54   0.59   0.66   0.31   9.25  62.59]]
    Vector Xi: 
    [-2.  4. -6.  2. -4.  7.]
    >>> 
    

    el gráfico de los iteraciones vs errores es:

    Gauss-Seidel errado_1EIIT2007T2

    La tabla obtiene aplicando la función de Gauss-Seidel, tomando como vector inicial el vector de ceros.

    Tarea: X=TX+C

    Instrucciones en Python

    # 1Eva_IIT2007_T2 Aplicar Gauss-Seidel
    # Algoritmo Gauss-Seidel
    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
        '''
        # Matriz aumentada
        nB = len(np.shape(B))
        if nB == 1:
            B = np.transpose([B])
        AB  = np.concatenate((A,B),axis=1)
        M = np.copy(AB)
        
        # Pivoteo por filas AB
        tamano = np.shape(M)
        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(M[i:,i])
            dondemax = np.argmax(columna)
            
            # dondemax no es en diagonal
            if (dondemax != 0):
                # intercambia filas
                temporal = np.copy(M[i,:])
                M[i,:] = M[dondemax+i,:]
                M[dondemax+i,:] = temporal
    
                pivoteado = pivoteado + 1
                if vertabla==True:
                    print(pivoteado, 'intercambiar: ',i,dondemax)
        if vertabla==True and pivoteado==0:
            print('Pivoteo por filas NO requerido')
        return(M)
    
    def gauss_seidel(A,B,tolera,X0, iteramax=100, vertabla=False, precision=5):
        ''' Método de Gauss Seidel, tolerancia, vector inicial X0
            para mostrar iteraciones: vertabla=True
        '''
        tamano = np.shape(A)
        n = tamano[0]
        m = tamano[1]
        diferencia = 2*tolera*np.ones(n, dtype=float)
        errado = np.max(diferencia)
        X = np.copy(X0)
    
        itera = 0
        if vertabla==True:
            print('Iteraciones Gauss-Seidel')
            print('itera,[X],[errado]')
            np.set_printoptions(precision)
            print(itera, X, np.array([errado]))
        while (errado>tolera and itera<iteramax):
            for i in range(0,n,1):
                xi = B[i]
                for j in range(0,m,1):
                    if (i!=j):
                        xi = xi-A[i,j]*X[j]
                xi = xi/A[i,i]
                diferencia[i] = np.abs(xi-X[i])
                X[i] = xi
            errado = np.max(diferencia)
            itera = itera + 1
            if vertabla==True:
                print(itera, X, np.array([errado]))        
        
        if (itera>iteramax): # No converge
            X = itera
            print('iteramax superado, No converge')
        return(X)
    
    # Programa de prueba #######
    # INGRESO
    A = np.array([[7.63, 0.30, 0.15,  0.50, 0.34, 0.84],
                  [0.38, 6.40, 0.70,  0.90, 0.29, 0.57],
                  [0.83, 0.19, 8.33,  0.82, 0.34, 0.37],
                  [0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
                  [0.71, 0.30, 0.85,  0.82, 5.95, 0.55],
                  [0.43, 0.54, 0.59,  0.66, 0.31, 9.25]])
    
    B = np.array([-9.44,25.27,-48.01,19.76,-23.63,62.59])
    
    tolera = 0.00001
    X = np.zeros(len(A), dtype=float)
    
    # PROCEDIMIENTO
    n = len(A)
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    
    
    AB = pivoteafila(A,B, vertabla=True)
    A = AB[:,:n]
    B = AB[:,n]
    respuesta = gauss_seidel(A,B, tolera, X, vertabla=True)
    
    # SALIDA
    print('Matriz aumentada y pivoteada:')
    print(AB)
    print('Vector Xi: ')
    print(respuesta)
    
    

    En el caso de la norma infinito, para la matriz A, se puede usar el algoritmo desarrollado en clase.
    Como valor para verificar su algoritmo, se obtuvo:

    >>> np.linalg.norm(A, np.inf)
    13.479999999999999
    

    Tarea: incluir la norma infinito para T

  • s1Eva_IIT2007_T1 Distribución binomial acumulada

    Ejercicio: 1Eva_IIT2007_T1 Distribución binomial acumulada

    Dado F=0.4, dado que n=5 y k=1

    F = \sum_{t=0}^{k} \binom{n}{t} p^t (1-p)^{n-t}

    La fórmula para el ejercicio se convierte en:

    F = \Bigg( \begin{array}{c} 5 \\ 0 \end{array} \Bigg) p ^0 (1-p)^{(5-0)} + \Bigg( \begin{array}{c} 5 \\ 1 \end{array} \Bigg) p ^1 (1-p)^{(5-1)} = 0.4

    Los valores de las combinatorias se calculan como:

    >>> import scipy.special as sts
    >>> sts.comb(5,0,repetition=False)
    1.0
    >>> sts.comb(5,1,repetition=False)
    5.0
    >>> 
    
    (1-p)^{5} + 5p (1-p)^{4} = 0.4

    La expresión para el ejercicio se convierte en:

    f(p) = (1-p)^{5} + 5p (1-p)^{4} - 0.4

    como referencia se revisa la gráfica para f(p)

    f(p) = (1-p)^5 + 5p(1-p)^4 - 0.4 = (1-p)^4 (1 - p + 5p) - 0.4 = (1-p)^4 (1 + 4p) - 0.4 = (1-p)^2 (1-p)^2 (1 + 4p) - 0.4 = (1-2p+p^2) (1-2p+p^2) (1 + 4p) - 0.4 = (1 - 4p + 6p^2 - 4p^3 +p^4 ) (1 + 4p) - 0.4 = 1 - 10p^2 + 20p^3 + 15p^4 + 4p^5 - 0.4 f(p) = 0.6 - 10p^2 + 20p^3 + 15p^4 + 4p^5

    y su derivada:

    f'(p) = - 20p + 60p^2 + 60p^3 +20p^4

    con lo que se puede desarrollar el método Newton-Raphson.

    Verificando el polinomio  obtenido a partir de la expresión inicial usando Sympy:

    >>> import sympy as sp
    >>> p = sp.Symbol('p')
    >>> poli = (1-p)**5 + 5*p*((1-p)**4) - 0.4
    >>> pol = poli.expand()
    >>> pol
    4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
    >>> pol.diff(p,1)
    20*p**4 - 60*p**3 + 60*p**2 - 20*p
    

    A partir de la gráfica, un punto inicial cercano a la raíz es X0 = 0.2

    itera = 0

    f(0.2) = 0.6 - 10(0.2)^2 + 20(0.2)^3 + 15(0.2)^4 + 4(0.2)^5 = 0.3373 f'(0.2) = - 20(0.2) + 60(0.2)^2 + 60(0.2)^3 +20(0.2)^4 = -2.048 x_1 = 0.2 -\frac{0.3373}{-2.048} = 0.3647 errado = |0.2 - 0.36469| = 0.1647

    itera = 1

    f(0.36469) = 0.6 - 10(0.36469)^2 + 20(0.36469)^3 + 15(0.36469)^4 + 4(0.36469)^5 = 0.0005566 f'(0.36469) = - 20(0.36469) + 60(0.36469)^2 + 60(0.36469)^3 +20(0.36469)^4 = -1.8703 x_1 = 0.36469 -\frac{0.0005566}{-1.8703} = 0.36499 errado = |0.36469 - 0.36499| = 0.000297

    itera = 3

    f(0.36499) = 0.6 - 10(0.36499)^2 + 20(0.36499)^3 + 15(0.36499)^4 + 4(0.36499)^5 = 1.6412291237166698e-07 f'(0.36499) = - 20(0.36499) + 60(0.36499)^2 + 60(0.36499)^3 +20(0.36499)^4 = -1.869204616112814 x_1 = 0.36469 -\frac{1.6412291237166698e-07}{ -1.8692} = 0.36499 errado = |0.36499 - 0.36499| = 8.7804e-08

    verificar con raíz: 0.3649852264049102

    Instrucciones para la gráfica

    # 1ra Eval II Término 2007
    # Tema 1. Distribución binomial acumulada
    import numpy as np
    import matplotlib.pyplot as plt
    import scipy.special as sts
    
    fp = lambda p: (1-p)**5 + 5*p*((1-p)**4) - 0.4
    
    a = 0
    b = 1
    pasos = 100
    
    # PROCEDIMIENTO
    xi = np.linspace(a,b,pasos+1)
    p_i = fp(xi)
    
    # SALIDA
    plt.plot(xi,p_i)
    plt.axhline(0)
    plt.show()
    

    Algoritmo en Python

    el resultado usando el algoritmo es:

    i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
    0 [ 0.2     0.3373 -2.048   0.3647  0.1647]
    1 [ 3.6469e-01  5.5668e-04 -1.8703e+00  3.6499e-01  2.9764e-04]
    2 [ 3.6499e-01  1.6412e-07 -1.8692e+00  3.6499e-01  8.7804e-08]
    raiz en:  0.3649852264049102
    

    con error de: 8.780360960525257e-08

    Instrucciones en Python para Newton-Raphson

    # 1Eva_IIT2007_T1 Distribución binomial acumulada
    # Método de Newton-Raphson
    
    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 p: 4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
    dfx = lambda p: 20*p**4 - 60*p**3 + 60*p**2 - 20*p
    
    x0 = 0.2
    tolera = 0.0000001
    
    # PROCEDIMIENTO
    respuesta = newton_raphson(fx,dfx,x0, tolera, vertabla=True)
    # SALIDA
    print('raiz en: ', respuesta)
    
    # GRAFICA
    a = 0
    b = 1
    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()
    
  • Medir el tiempo de ejecución en Python

    Con varios métodos para resolver un problema, una media de eficiencia es el tiempo de ejecución de un algoritmo.

    dibujo cronómetro y banderaPara determinar el tiempo de ejecución de un algoritmo, se toman las lecturas de tiempo antes y después del bloque, calculando el tiempo de ocupación de las diferencias entre tiempos.

    La librería "time" permite obtener lectura del reloj del computador.

    Como un algoritmo de prueba se usa la suma de los m primeros números enteros.

    # tiempos de ejecuci n de algoritmos
    import time
    
    # INGRESO
    m = 100
    
    # PROCEDIMIENTO
    t1 = time.clock()
    
    # Ejecuta Operaciones del algoritmo
    suma = 0
    for i in range(0,m,1):
        suma = suma + i
    
    t2 = time.clock()
    # Tiempo para ejecutar operaciones
    ocupado = t2-t1
    
    # SALIDA
    print('tiempos:')
    print(t1, t2)
    print('tiempo de ejecuci n:')
    print(ocupado)
    

    con lo que se obtienen los siguientes resultados:

    tiempos:
    8.210955878428587e-07 5.46028565915501e-05
    tiempo de ejecución:
    5.378176100370724e-05
    >>> 
    
  • Sympy - Fórmulas y funciones simbólicas

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]

    SymPy es una librería usada para manejar de forma algebraica las expresiones matemáticas. Se requiere definir el símbolo que representa la variable en la expresión, por ejemplo la letra 'x'. Si la librería Sympy no está disponible o muestra un error la puede revisar las instrucciones del enlace instalar con pip.

    Una formula o función matemática f(x) descrita como fx se puede derivar, integrar, simplificar sus términos. También se puede construir una expresión matemática, por ejemplo desarrollar los términos para un polinomio como Taylor.

    Referencia: https://www.sympy.org/es/

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    1. Expresiones con Sympy

    Una función f(x) como un polinomio

    f(x) = (1-x)^5 + 5 x ((1-x)^4) - 0.4

    Se pueden manejar de forma algebraica con Sympy al declarar la variable independiente 'x' como un símbolo. La expresión se asigna a fx

    import sympy as sym
    x = sym.Symbol('x')
    fx = (1-x)**5 + 5*x*((1-x)**4) - 0.4
    

    Se simplifican o se expande los términos del polinomio con solo una instrucción,

    fx = fx.expand()
    print(fx)
    

    mostrando el siguiente resultado:

    4*x**5 - 15*x**4 + 20*x**3 - 10*x**2 + 0.6
    

    o una presentación diferente con sym.pprint():

    >>> sym.pprint(fx)
               4          5      
    5*x*(1 - x)  + (1 - x)  - 0.4
    >>> 
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/printing.html#setting-up-pretty-printing

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    2. Operaciones, combinar otras expresiones

    A una función simbólica fx se pueden añadir más términos con el mismo símbolo

    >>> import sympy as sym
    >>> sym.Symbol('x')
    >>> fx = sym.cos(x)
    >>> gx = fx + x
    >>> gx
    x + cos(x)
    >>> gx = gx - 2
    >>> gx
    x + cos(x) - 2
    

    Por lo que las funciones simbólicas son útiles cuando se construyen expresiones, como en el caso de series de Taylor descritas en la Unidad01

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    3. Evaluar una expresión f(x) con Sympy

    Las expresiones simbólicas se pueden evaluar en un punto, ejemplo x0=0.1 usando la instrucción fx.subs(x,x0), que sustituye el símbolo de la variable x con el valor definido para x0.

    >>> import sympy as sym
    >>> sym.Symbol('x')
    >>> fx = sym.cos(x)
    >>> fx.subs(x,0.1)
    0.995004165278026
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/basic_operations.html#substitution

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    4. Conversión de forma simbólica a forma numérica Lambda

    Para evaluar varios puntos en la expresión  en forma numérica para usar 'x' con valores en un vector o matriz como un arreglo, se convierte a la forma numérica Lambda con la instrucción sym.lambdify(x,fx)

    >>> >>> import sympy as sym
    >>> sym.Symbol('x')
    >>> fx = sym.cos(x)
    &>>> f_x = sym.lambdify(x,fx)
    >>> f_x(0.1)
    0.9950041652780258
    >>> f_x([0.1,0.3,0.5])
    array([0.99500417, 0.95533649, 0.87758256])
    >>> 
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/basic_operations.html#lambdify

    Ejemplo: Polinomio de Taylor – Ejemplos con Sympy-Python

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    5. expresiones de Funciones matemáticas

    Las Funciones matemáticas usadas en otras librerías como Numpy tienen también su representación simbólica en Sympy. Ejemplo cos(x):

    f(x) = \cos (x)
    >>> import sympy as sym
    >>> x = sym.Symbol('x')
    >>> fx = sym.cos(x)
    

    Realice pruebas con otras funciones: sin(x), exp(x), log(x), log10(x), Heaviside(x)

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    6. Derivadas con Sympy

    Las expresiones de la derivada se obtienen con la expresión fx.diff(x,k), indicando la k-ésima  derivada de la expresión.

    >>> import sympy as sym
    >>> x = sym.Symbol('x')
    >>> fx = sym.cos(x)
    >>> x
    x
    >>> fx
    cos(x)
    >>> fx.diff(x,1)
    -sin(x)
    >>> fx.diff(x,2)
    -cos(x)
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#derivatives

    Ejemplos: Sistemas LTI CT – Respuesta a entrada cero con Sympy-Python

    Derivadas como expresión de Sympy sin evaluar

    Cuando se requiere expresar tan solo la operación de derivadas, que será luego usada o reemplazada con otra expresión, se usa la derivada sin evaluar. Ejemplo:

    y = \frac{d}{dx}f(x)

    Para más adelante definir f(x).

    En Sympy, la expresión de y se realiza indicando que f será una variable

    x = sym.Symbol('x', real=True)
    f = sym.Symbol('f', real=True)
    y = sym.diff(f,x, evaluate=False) # derivada sin evaluar
    g = sym.cos(x) + x**2
    yg = y.subs(f,g).doit() # sustituye f con g y evalua .doit()
    
    >>> y
    Derivative(f, x)
    >>> g
    x**2 + cos(x)
    >>> yg
    2*x - sin(x)

    Ejemplos:

    Polinomio de Taylor – Ejemplos con Sympy-Python

    EDO lineal – solución complementaria y particular con Sympy-Python

    EDO lineal – ecuaciones auxiliar, general y complementaria con Sympy-Python

    EDP Parabólica – analítico con Sympy-Python

    EDP Elípticas – analítico con Sympy-Python

    Evaluación de las expresiones Sympy

    Se define la expresión 'derivada' y se la usa con la instrucción fx.subs(x,valor)

    >>> fx.subs(x,0)
    1
    >>> fx.subs(x,1/3)
    0.944956946314738
    >>> derivada = fx.diff(x,1)
    >>> derivada
    -sin(x)
    >>> derivada.subs(x,1/3)
    -0.327194696796152
    >>> 
    

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]
    ..


    7. Integrales con Sympy

    Sympy incorpora la operación del integral en sus expresiones, que pueden ser integrales definidas en un intervalo o expresiones sin evaluar.

    >>> import sympy as sym
    >>> t = sym.Symbol('t',real=True)
    >>> x = 10*sym.exp(-3*t)
    >>> x
    10*exp(-3*t)
    >>> y = sym.integrate(x,(t,0,10))
    >>> y
    10/3 - 10*exp(-30)/3
    >>> y = sym.integrate(x,(t,0,sym.oo))
    >>> y
    10/3
    

    Referencia: https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#integrals

    Ejemplos:
    Transformada de Laplace – Integral con Sympy-Python

    Series de Fourier de señales periódicas con Sympy-Python

    Integral de Convolución entre x(t) y h(t) causal con Sympy-Python

    [ expresiones ] [ operaciones ] [ evaluar f(x) ] [ f(x) a lambda ] [ f(x) matemáticas ] [ derivadas ] [ integrales ]

  • Numpy - Vectores y matrices

    Resumen de algunas funciones para vectores y matrices como arreglos en la librería Numpy de Python, usadas en el curso para facilitan el cálculo numérico. El orden de las instrucciones es el que aparece en las entradas del blog.

    Lo primero es hacer el llamado a las librerías con el alias 'np'

    import numpy as np
    

    Vectores

    puntos espaciados entre [a,b] - np.linspace()

    para obtener puntos en el rango [a,b] para una cantidad de tramos. Se obtienen tramos+1 puntos .

    a = 1
    b = 2
    tramos = 4
    x = np.linspace(a,b,tramos+1)
    
    >>> x
    array([ 1.  ,  1.25,  1.5 ,  1.75,  2.  ])
    

    rango de valores en vector - np.arange(a,b,dt)

    crea un vector con valores en el rango [a,b) y espaciados dt.

    t=np.arange(0,10,2)
    >>> t
    array([0, 2, 4, 6, 8])

    transponer vector - np.transpose()

    >>> fila = np.array([1,2,3])
    >>> columna = np.transpose([fila])
    >>> columna
    array([[1],
           [2],
           [3]])
    

    indices para ordenar - np.argsort()

    para ordenar por una columna específica:
    - se obtiene un vector con la columna que se requiere ordenar: tabla[:,0]
    - con el vector se determinan los índices para ordenar la tabla por filas: np.argsort()
    - se aplica los índices a una copia de la matriz: tabla[orden]

    tabla = np.array([[5,3],
                      [4,2],
                      [3,1],
                      [2,4],
                      [1,5]])
    
    # ordenar por primera columna
    referencia = tabla[:,0]
    orden = np.argsort(referencia)
    ordenada = tabla[orden]
    
    print(orden)
    print(ordenada)
    

    resultado:

    [4 3 2 1 0]
    [[1 5]
     [2 4]
     [3 1]
     [4 2]
     [5 3]]
    
    # ordenar por segunda columna
    referencia = tabla[:,1]
    orden = np.argsort(referencia)
    ordenada = tabla[orden]
    print(orden)
    print(ordenada)
    

    resultado:

    [2 1 0 3 4]
    [[3 1]
     [4 2]
     [5 3]
     [2 4]
     [1 5]]
    

    Matrices en Numpy

    concatenar- np.concatenate((A,b), axis=1)

    concatenar usa el parámetro axis: 0 para filas, 1 para columnas.

    import numpy as np
    cantidad = np.array([[4,2,5],
                         [2,5,8],
                         [5,4,3]])
    
    pagado = np.array([[60.70],
                       [92.90],
                       [56.30]])
    
    matriz = np.concatenate((cantidad,pagado),axis=1)
    
    >>> matriz
    array([[  4. ,   2. ,   5. ,  60.7],
           [  2. ,   5. ,   8. ,  92.9],
           [  5. ,   4. ,   3. ,  56.3]])
    

    resolver sistema de ecuaciones - np.linalg.solve(A,B)

    Resuelve el sistema de ecuaciones dado por una matriz A y un vector B. siendo, por ejemplo:

     0 =  c1 +  c2
    -5 = -c1 - 2c2
    
    c1 = -5 
    c2 = 5
    >>> A = [[ 1, 1],
    	 [-1,-2]]
    >>> B = [0,-5]
    >>> np.linalg.solve(A,B)
    array([-5.,  5.])
    >>>

    Otras instrucciones

    np.pi constante con valor π

    >>> np.pi
    3.141592653589793
    np.sin(t)
    np.cos(t)
    función trigonométrica en radianes. La variable t puede ser un escalar o un arreglo.

    >>> t=0.65
    >>> np.sin(0.65)
    0.60518640573603955
    >>> t=[0, 0.3, 0.6]
    >>> np.sin(t)
    array([ 0. , 0.29552021, 0.56464247])
    np.abs() obtiene el valor absoluto de un número. En el caso de un número complejo obtiene la parte real.
    np.real(complejo)
    np.imag(complejo)
    obtiene la parte real de los números complejos en un vector. Se aplica lo mismo para la parte imaginaria del número complejo.
    complex(a,b) crea el número complejo a partir de los valores de a y b.
    a=2
    b=3
    el resultado es: 2+3j
    np.piecewise(t, t>=donde, [1,0]) función que crea a partir de t, los valores de la condición t>=donde, ubicando los valores de 1, para otro caso es 0. Usada en la función escalón.
    np.roots([a,b,c]) obtiene las raíces del polinomio:
    ax2+bx+c
    x2 + 3 x + 2 = (x+1)(x+2)

    >>> np.roots([1,3,2])
    array([-2., -1.])
  • Funciones lambda vs def-return

    Las dos formas de escritura de funciones matemáticas básicamente hacen lo mismo. Sin embargo, cuando la función se define por tramos, la forma def-return se convierte en la más usada.

    Se adjunta un ejemplo:

    Funciones Lambda

    Permite describir funciones sencillas de una sola línea, que no está conformada por intervalos.

    f(x) = x \cos(x) - 2 x^2 + 3 x -1
    import numpy as np
    fx = lambda x: x*np.cos(x) - 2*(x**2) + 3*x -1
    
    >>> fx(0.2)
    -0.28398668443175157
    

    Funciones def-return

    Permiten describir en detalle lo que sucede con una función matemática por intervalos. Tiene la ventaja de permitir definir valores por intervalos, puntos discontínuos, etc.

    f(x)= \begin{cases} x \cos(x) - 2 x^2 + 3 x -1, & x>0 \\1, & x \leq 0 \end{cases}
    import numpy as np
    def fx(x):
        if x>0:
            fi = x*np.cos(x) - 2*(x**2) + 3*x -1
        if x<=0:
            fi = 1 
        return(fi)
    
    >>> fx(0.2)
     -0.28398668443175157
    
  • 7.3 EDP hiperbólicas con Python

    [ concepto ] [ analítico ] [ algoritmo ]
    ..


    1. Concepto y ejercicio

    Referencia:  Chapra PT8.1 p860,  Rodríguez 10.4 p435

    Las Ecuaciones Diferenciales Parciales tipo hiperbólicas semejantes a la mostrada, para un ejemplo en particular, representa la ecuación de onda de una dimensión u[x,t], que representa el desplazamiento vertical de una cuerda.

    \frac{\partial ^2 u}{\partial t^2}=c^2\frac{\partial ^2 u}{\partial x^2}

    EDP Cuerda 01

    Los extremos de la cuerda de longitud unitaria están sujetos y referenciados a una posición 0 a la izquierda y 1 a la derecha.

    u[x,t] , 0<x<1, t≥0
    u(0,t) = 0 , t≥0
    u(1,t) = 0 , t≥0

    Al inicio, la cuerda está estirada por su punto central:

    u(x,0) = \begin{cases} -0.5x &, 0\lt x\leq 0.5 \\ 0.5(x-1) &, 0.5\lt x \lt 1 \end{cases}

    Se suelta la cuerda, con velocidad cero desde la posición inicial:

    \frac{\partial u(x,0)}{\partial t} = 0

    [ concepto ] [ analítico ] [ algoritmo ]
    ..


    2. Desarrollo analítico

    La solución se realiza de forma semejante al procedimiento para EDP parabólicas y elípticas. Se cambia a la forma discreta  usando diferencias finitas divididas:

    \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta t)^2} =c^2 \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

    El error es del orden O(\Delta x)^2 + O(\Delta t)^2.
    se reagrupa de la forma:

    u_{i,j+1}-2u_{i,j}+u_{i,j-1} = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \big( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \big)

    El término constante, por facilidad se simplifica con el valor de 1

    \lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1

    si c = 2 y Δx = 0.2, se deduce que Δt = 0.1

    que al sustituir en la ecuación, se simplifica anulando el término u[i,j]:

    u_{i,j+1}+u_{i,j-1} = u_{i+1,j}+u_{i-1,j}

    EDP Cuerda 02

    en los que intervienen solo los puntos extremos. Despejando el punto superior del rombo:

    u_{i,j+1} = u_{i+1,j}-u_{i,j-1}+u_{i-1,j}

    Convergencia:

    \lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \leq 1

    para simplificar aún más el problema, se usa la condición velocidad inicial de la cuerda igual a cero

    \frac{\delta u_{i,0}}{\delta t}=\frac{u_{i,1}-u_{i,-1}}{2\Delta t} = 0 u_{i,-1}=u_{i,1}

    que se usa para cuando el tiempo es cero, permite calcular los puntos para t[1]:

    j=0

    u_{i,1} = u_{i+1,0}-u_{i,-1}+u_{i-1,0} u_{i,1} = u_{i+1,0}-u_{i,1}+u_{i-1,0} 2 u_{i,1} = u_{i+1,0}+u_{i-1,0} u_{i,1} = \frac{u_{i+1,0}+u_{i-1,0}}{2}

    Aplicando solo cuando j = 0

    EDP Cuerda 03

    que al ponerlos en la malla se logra un sistema explícito para cada u[i,j], con lo que se puede resolver con un algoritmo:

    EDP hiperbolica 01

    [ concepto ] [ analítico ] [ algoritmo ]
    ..


    3. Algoritmo con Python

    # Ecuaciones Diferenciales Parciales
    # Hiperbólica. Método explicito
    import numpy as np
    
    def cuerdainicio(xi):
        n = len(xi)
        y = np.zeros(n,dtype=float)
        for i in range(0,n,1):
            if (xi[i]<=0.5):
                y[i]=-0.5*xi[i]
            else:
                y[i]=0.5*(xi[i]-1)
        return(y)
    
    # INGRESO
    x0 = 0 # Longitud de cuerda
    xn = 1
    y0 = 0 # Puntos de amarre
    yn = 0
    t0 = 0 # tiempo inicial
    c = 2  # constante EDP
    # discretiza
    tramosx = 16
    tramost = 32
    dx = (xn-x0)/tramosx 
    dt = dx/c
    
    # PROCEDIMIENTO
    xi = np.arange(x0,xn+dx/2,dx)
    tj = np.arange(0,tramost*dt+dt/2,dt)
    n = len(xi)
    m = len(tj)
    
    u = np.zeros(shape=(n,m),dtype=float)
    u[:,0] = cuerdainicio(xi)
    u[0,:] = y0
    u[n-1,:] = yn
    # Aplicando condición inicial
    j = 0
    for i in range(1,n-1,1):
        u[i,j+1] = (u[i+1,j]+u[i-1,j])/2
    # Para los otros puntos
    for j in range(1,m-1,1):
        for i in range(1,n-1,1):
            u[i,j+1] = u[i+1,j]-u[i,j-1]+u[i-1,j]
    
    # SALIDA
    np.set_printoptions(precision=2)
    print('xi =')
    print(xi)
    print('tj =')
    print(tj)
    print('matriz u =')
    print(u)
    

    con lo que se obtienen los resultados numéricos, que para mejor interpretación se presentan en una gráfica estática y otra animada.

    # GRAFICA
    import matplotlib.pyplot as plt
    
    for j in range(0,m,1):
        y = u[:,j]
        plt.plot(xi,y)
    
    plt.title('EDP hiperbólica')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.show()
    
    # **** GRÁFICO CON ANIMACION ***********
    import matplotlib.animation as animation
    
    # Inicializa parametros de trama/foto
    retardo = 70   # milisegundos entre tramas
    tramas  = m
    maximoy = np.max(np.abs(u))
    figura, ejes = plt.subplots()
    plt.xlim([x0,xn])
    plt.ylim([-maximoy,maximoy])
    
    # lineas de función y polinomio en gráfico
    linea_poli, = ejes.plot(xi,u[:,0], '-')
    plt.axhline(0, color='k')  # Eje en x=0
    
    plt.title('EDP hiperbólica')
    # plt.legend()
    # txt_x = (x0+xn)/2
    txt_y = maximoy*(1-0.1)
    texto = ejes.text(x0,txt_y,'tiempo:',
                      horizontalalignment='left')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    
    # Nueva Trama
    def unatrama(i,xi,u):
        # actualiza cada linea
        linea_poli.set_ydata(u[:,i])
        linea_poli.set_xdata(xi)
        linea_poli.set_label('tiempo linea: '+str(i))
        texto.set_text('tiempo['+str(i)+']')
        # color de la línea
        if (i<=9):
            lineacolor = 'C'+str(i)
        else:
            numcolor = i%10
            lineacolor = 'C'+str(numcolor)
        linea_poli.set_color(lineacolor)
        return linea_poli, texto
    
    # Limpia Trama anterior
    def limpiatrama():
        linea_poli.set_ydata(np.ma.array(xi, mask=True))
        linea_poli.set_label('')
        texto.set_text('')
        return linea_poli, texto
    
    # Trama contador
    i = np.arange(0,tramas,1)
    ani = animation.FuncAnimation(figura,
                                  unatrama,
                                  i ,
                                  fargs=(xi,u),
                                  init_func=limpiatrama,
                                  interval=retardo,
                                  blit=True)
    # Graba Archivo video y GIFAnimado :
    # ani.save('EDP_hiperbólica.mp4')
    ani.save('EDP_hiperbolica.gif', writer='imagemagick')
    plt.draw()
    plt.show()
    

    Una vez realizado el algoritmo, se pueden cambiar las condiciones iniciales de la cuerda y se observan los resultados.

    Se recomienda realizar otros ejercicios de la sección de evaluaciones para EDP Hiperbólicas y observar los resultados con el algoritmo modificado.

    [ concepto ] [ analítico ] [ algoritmo ]