Categoría: Soluciones

Ejercicios resueltos de examen

  • s2Eva_2023PAOII_T2 EDO Cable cuelga entre apoyos A y B

    Ejercicio: 2Eva_2023PAOII_T2 EDO Cable cuelga entre apoyos A y B

    Literal a

    La ecuación diferencial a resolver es:

    \frac{d^2y}{dx^2} = \frac{w_0}{T_0} \Big[ 1+ \sin \Big(\frac{\pi x}{2l_B} \Big) \Big]

    donde w0 = 1 000 lbs/ft, T0. = 0.588345×106.
    dy(0)/dx = 0 y lB=200 de la gráfica presentada.

    \frac{d^2y}{dx^2} = \frac{1000}{0.588345×10^6} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big]

    Para usar Runge Kutta para segunda derivada:

    z= y' = f(x,y,z) z' = (y')' = 0z + \frac{1000}{0.588345×10^6} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big] g(x,y,z) = \frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big]

    los valores iniciales para el ejercicio acorde al enunciado son: x0 = 0, y0=0, z0 = 0, con h=0.5

    literal b

    para itera 0

    K1y = h*z = 0.5*0 = 0 K1z = (0.5)\frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi (0)}{2(200)} \Big) \Big] = 0.0008498 K2y = h*(z+K1z) = (0.5) (0+0.00084984) = 0.0004249 K2z = (0.5)\frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi (0+0.5)}{2(200)} \Big) \Big] =0.0008531 y = 0+\frac{0+0.0004249}{2} = 0.0002124 z = 0+\frac{0.0008498+0.0008531}{2} = 0.0008515 x = 0 + 0.5 = 0.5

    ...

    Desarrollar dos iteraciones adicionales como tarea.

    Para las primeras iteraciones de un total de 400+1, los valores con Python y en resultados.txt :

    estimado[xi,yi,zi,K1y,K2y,K1z,K2z]
    [0.0000 0.0000e+00 0.0000e+00 
    0.0000e+00 0.0000e+00 
    0.0000e+00 0.0000e+00]
    
    [0.5000 2.124603761398499073e-04 8.515101601627471980e-04 
    0.000000000000000000e+00 4.249207522796998146e-04 
    8.498415045593996292e-04 8.531788157660947667e-04]
    
    [1.0000 8.515101601627470896e-04 1.706357605799455881e-03 
    4.257550800813735990e-04 8.523444879644209281e-04 
    8.531788157660947667e-04 8.565160755073224913e-04]
    
    [1.5000 1.918817981939305653e-03 2.564542259712321911e-03 
    8.531788028997279406e-04 1.281436840653389295e-03 
    8.565160755073224913e-04 8.598532323184093513e-04]
    
    [2.0000 3.416052419875068892e-03 3.426063993239661116e-03 
    1.282271129856160955e-03 1.712197746015365740e-03 
    8.598532323184093513e-04 8.631902347362692754e-04]
    ...
    

    literal c

    resultado en archivo.txt al ejecutar el algoritmo.

    literal d

    cable entre apoyos A y B

    Algoritmo con Python

    # 2Eva_2023PAOII_T2 Cable cuelga entre apoyos A y B
    import numpy as np
    
    def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,3+4),dtype=float)
    
        # incluye el punto [x0,y0,z0]
        estimado[0] = [x0,y0,z0,0,0,0,0]
        xi = x0
        yi = y0
        zi = z0
        for i in range(1,tamano,1):
            K1y = h * f(xi,yi,zi)
            K1z = h * g(xi,yi,zi)
            
            K2y = h * f(xi+h, yi + K1y, zi + K1z)
            K2z = h * g(xi+h, yi + K1y, zi + K1z)
    
            yi = yi + (K1y+K2y)/2
            zi = zi + (K1z+K2z)/2
            xi = xi + h
            
            estimado[i] = [xi,yi,zi,K1y,K2y,K1z,K2z]
        return(estimado)
    
    # PROGRAMA PRUEBA
    # Ref Rodriguez 9.1.1 p335 ejemplo.
    # prueba y'-y-x+(x**2)-1 =0, y(0)=1
    
    # INGRESO
    T0 = 0.588345e6
    LB = 200
    f = lambda x,y,z: z
    g = lambda x,y,z: (1000/T0)*(1+np.sin(np.pi*x/(2*LB)))
    x0 = 0
    y0 = 0
    z0 = 0
    h  = 0.5
    muestras = 401
    
    # PROCEDIMIENTO
    puntosRK2 = rungekutta2_fg(f,g,x0,y0,z0,h,muestras)
    xi = puntosRK2[:,0]
    yiRK2 = puntosRK2[:,1]
    
    # SALIDA
    np.set_printoptions(precision=4)
    print('estimado[xi,yi,zi,K1y,K2y,K1z,K2z]')
    print(puntosRK2)
    np.savetxt("tablaRk2.txt",puntosRK2)
    
    
    # Gráfica
    import matplotlib.pyplot as plt
    
    plt.plot(xi[0],yiRK2[0],
             'o',color='r', label ='[x0,y0]')
    plt.plot(xi[1:],yiRK2[1:],
             color='m',
             label ='y Runge-Kutta 2 Orden')
    
    plt.title('EDO: Solución con Runge-Kutta 2do Orden')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()
    

     

  • s2Eva_2023PAOII_T1 Volumen por solido de revolución

    Ejercicio: 2Eva_2023PAOII_T1 Volumen por solido de revolución

    Integra Revolucion 02El volumen se calcula a partir de la expresión:

    V = \int_{a}^{b} \pi (f(x))^2 dx

    literal a y c

    Para el volumen con f(x) con al menos 3 tramos y un método de Simpson, directamente se puede usar 3/8. Por lo que se Se reemplaza en la fórmula de volumen del sólido de revolución f(x) con:

    f(x) = \sqrt{\sin (x/2)}

    obteniendo:

    V_{fx} = \int_{a}^{b} \pi \Big(\sqrt{\sin (x/2)} \Big)^2 dx = \int_{a}^{b} \pi \sin (x/2) dx

    La expresión dentro del integral se denomina como fv:

    f_v (x)= \pi \sin (x/2)

    en el intervalo [0.1, 1.8],  con al menos 3 tramos, se requieren 4 muestras con tamaño de paso hf: y truncando a 4 decimales los resultados calculados con Python.

    h_f =\frac{b-a}{tramos} = \frac{1.8-0.1}{3} = 0.5666

    los puntos de muestra quedan np.linspace(0.1,1.8,3+1):

    xis= [0.1, 0.6666, 1.2333, 1.8 ]

    El integral se calcula con los puntos de muestra,

    V_{fx} = \frac{3}{8} (0.5666) \Big( f_v(0.1) +3 f_v(0.6666) + + 3 f_v(1.2333)+ f_v(1.8)\Big)

    recordando que se usa en radianes,

    V_{fx} = \frac{3}{8} (0.5666) \Bigg( \pi \sin \Big(\frac{0.1}{2}\Big) +3 \pi \sin \Big(\frac{0.6666}{2}\Big) + + 3 \pi \sin\Big(\frac{1.2333}{2}\Big)+ \pi \sin \Big(\frac{1.8}{2}\Big)\Bigg) = \frac{3}{8} (0.5666) \Big( 0.1570+3 (1.0279) + + 3 (1.8168)+ 2.4608\Big)

    literal d. el volumen generado por f(x) tiene como resultado:

    V_{fx} = 2.3698

    la cota de error para fx es el orden de O(h5) = O(0.56665) = O(0.05843), queda como tarea completar la cota de error total.

    literal b y c

    Para el volumen con g(x) con al menos 2 tramos y Cuadratura de Gauss de dos puntos, se reemplaza en la fórmula de volumen de sólido de revolución:

    g(x) = e^{x/3} - 1 V_{gx} = \int_{a}^{b} \pi (e^{x/3} - 1)^2 dx

    La expresión dentro del integral se denomina como gv:

    g_v = \pi (e^{x/3} - 1)^2

    en el intervalo [0.1, 1.8],  con al menos 2 tramos, se requieren 3 muestras con tamaño de paso hg:

    h_g =\frac{b-a}{tramos} = \frac{1.8-0.1}{2} = 0.85

    xic = [0.1, 0.95, 1.8 ]

    tramo 1: [0.1, 0.95] , a = 0.1 , b= 0.95, truncando a 4 decimales

    x_a = \frac{0.95+0.1}{2} - \frac{0.95-0.1}{2}\frac{1}{\sqrt{3}} = 0.2796 x_b = \frac{0.95+0.1}{2} + \frac{0.95-0.1}{2}\frac{1}{\sqrt{3}} = 0.7703 g_v(0.2796) = \pi (e^{0.2796/3} - 1)^2 = 0.02998 g_v(0.7703) = \pi (e^{0.7703/3} - 1)^2 = 0.2692 V_{c1} = \frac{0.95-0.1}{2}(g_v(0.2796) + g_v(0.7703)) V_{c1} = \frac{0.95-0.1}{2}(0.02998 + 0.2692) V_{c1} = 0.1271

    tramo 2: [0.95, 1.8] , a = 0.95 , b= 1.8

    x_a = \frac{1.8+0.95}{2} - \frac{1.8-0.95}{2}\frac{1}{\sqrt{3}} = 1.1296 x_b = \frac{1.8+0.95}{2} - \frac{1.8-0.95}{2}\frac{1}{\sqrt{3}} = 1.6203 g_v(1.1296) = \pi (e^{1.1296/3} - 1)^2 = 0.6567 g_v(1.6203) = \pi (e^{1.6203/3} - 1)^2 = 1.6115 V_{c2} = \frac{1.8-0.95}{2}(g_v(1.1296) + g_v(1.6203)) V_{c2} = \frac{1.8-0.95}{2}(0.6567 + 1.6115) V_{c2} = 0.9640

    literal d. volumen generado por g(x)

    V_{gx} = V_{c1} + V_{c2} = 0.1271 + 0.9640 = 1.0912

    completar la cota de error para cuadratura de Gauss de dos puntos.

    literal e. El volumen de revolución se genera como la resta del volumen de f(x) y volumen g(x)

    V = V_{fx} - V_{gx} = 2.3698 - 1.0912 = 1.2785

    Algoritmo con Python

    Los resultados usando el algoritmo con las operaciones usadas en el planteamiento son:

    para f(x):
    xis= [0.1        0.66666667 1.23333333 1.8       ]
    fiv= [0.15701419 1.02791246 1.81684275 2.46089406]
    Volumenfx:  2.369836948864926
    
    para g(x):
    Por tramos: [0.1  0.95 1.8 ]
    xab= [0.2796261355944091, 0.770373864405591,
          1.129626135594409, 1.620373864405591]
    gab= [0.02998177327492598, 0.26928904479784566,
          0.6567986343358181, 1.6115494735531555]
    Vc1= 0.12719009768092793  ; Vc2= 0.964047945852814
    Volumengx:  1.0912380435337419
    
    Volumen solido revolucion: 1.2785989053311841

    Considerando realizar los cálculos para cada sección:

    # 2Eva_2023PAOII_T1 Volumen por solido de revolución
    import numpy as np
    
    # INGRESO
    fx = lambda x: np.sqrt(np.sin(x/2))
    gx = lambda x: np.exp(x/3)-1
    a = 0.1
    b = 1.8
    tramosSimpson = 3
    tramosCGauss = 2
    
    # PROCEDIMIENTO
    # Volumen para f(x) con Simpson
    fv = lambda x: np.pi*np.sin(x/2)
    hs = (b-a)/tramosSimpson
    xis = np.linspace(a,b,tramosSimpson +1)
    fiv = fv(xis)
    Vs = (3/8)*hs*(fiv[0]+3*fiv[1]+3*fiv[2]+ fiv[3])
    
    # Volumen para g(x) con Cuadratura de Gauss
    gv = lambda x: np.pi*(np.exp(x/3)-1)**2
    hc = (b-a)/tramosSimpson
    xic = np.linspace(a,b,tramosCGauss +1)
    # tramo 1
    ac = xic[0]
    bc = xic[1]
    xa = (bc+ac)/2 + (bc-ac)/2*(-1/np.sqrt(3)) 
    xb = (bc+ac)/2 + (bc-ac)/2*(1/np.sqrt(3))
    Vc1 = (bc-ac)/2*(gv(xa)+gv(xb))
    xab = [xa,xb]
    gab = [gv(xa),gv(xb)]
    # tramo 2
    ac = xic[1]
    bc = xic[2]
    xa = (bc+ac)/2 + (bc-ac)/2*(-1/np.sqrt(3)) 
    xb = (bc+ac)/2 + (bc-ac)/2*(1/np.sqrt(3))
    Vc2 = (bc-ac)/2*(gv(xa)+gv(xb))
    Vc = Vc1+Vc2
    xab.append(xa)
    xab.append(xb)
    gab.append(gv(xa))
    gab.append(gv(xb))
    
    # Volumen solido revolucion
    Volumen = Vs - Vc
    
    # SALIDA
    print("para f(x):")
    print("xis=", xis)
    print("fiv=", fiv)
    print("Volumenfx: ",Vs)
    print()
    print("para g(x):")
    print("Por tramos:",xic)
    print("xab=", xab)
    print("gab=", gab)
    print("Vc1=",Vc1," ; Vc2=",Vc2) 
    print("Volumengx: ",Vc)
    print()
    print("Volumen solido revolucion:",Volumen)
    

    para la gráfica presentada en el enunciado (no requerida) , se complementa con las instrucciones:

    Integra Revolucion 02

    # para grafica -------------------
    import matplotlib.pyplot as plt
    muestras = 21 # grafica
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    gi = gx(xi)
    xig = np.linspace(a,b,tramosCGauss+1)
    fis = fx(xis)
    gig = gx(xig)
    
    # grafica
    plt.plot(xi,fi, label="f(x)")
    plt.plot(xi,gi, label="g(x)")
    plt.plot([0.0,2.0],[0,0], marker=".", color="blue")
    plt.fill_between(xi,fi,gi,color="lightgreen")
    plt.axhline(0)
    plt.axvline(a, linestyle="dashed")
    plt.axvline(b, linestyle="dashed")
    plt.xlabel('x')
    plt.ylabel('f(x), g(x)')
    plt.legend()
    plt.plot(xis,fis,'.b')
    plt.plot(xig,gig,'.r')
    plt.grid()
    plt.show()

    Gráfica de sólido de revolución en 3D

    sólido de revolución 3D

    Instrucciones en Python

    # 2Eva_2023PAOII_T1 Volumen por solido de revolución
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import axes3d
    
    # INGRESO
    f = lambda x: np.sqrt(np.sin(x/2))
    g = lambda x: np.exp(x/3)-1
    
    # eje x
    xa = 0.1
    xb = 1.8
    xmuestras = 31
    # angulo w de rotación
    w_a = 0
    w_b = 2*np.pi
    w_muestras = 31
    
    # PROCEDIMIENTO
    # muestreo en x y angulo w
    xi = np.linspace(xa, xb, xmuestras)
    wi = np.linspace(w_a, w_b, w_muestras)
    X, W = np.meshgrid(xi, wi)
    
    # evalua f(x) en 3D
    Yf = f(xi)*np.cos(W)
    Zf = f(xi)*np.sin(W)
    
    # evalua g(x) en 3D
    Yg = g(xi)*np.cos(W)
    Zg = g(xi)*np.sin(W)
    
    # SALIDA
    
    # grafica 3D
    figura = plt.figure()
    grafica = figura.add_subplot(111, projection='3d')
    
    grafica.plot_surface(X, Yf, Zf,
                         color='blue', label='f(x)',
                         alpha=0.3, rstride=6, cstride=12)
    grafica.plot_surface(X, Yg, Zg,
                         color='orange', label='g(x)',
                         alpha=0.3, rstride=6, cstride=12)
    
    grafica.set_title('Solidos de revolución')
    grafica.set_xlabel('x')
    grafica.set_ylabel('y')
    grafica.set_zlabel('z')
    # grafica.legend()
    eleva = 45
    rota = 45
    deltaw = 5
    grafica.view_init(eleva, rota)
    
    # rotacion de ejes
    for angulo in range(rota, 360+rota, deltaw ):
        grafica.view_init(eleva, angulo)
        plt.draw()
        plt.pause(.001)
    plt.show()
    
  • s2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

    Ejercicio: 2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

    Dada la EDP elíptica,

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = \Big( x^2 + y^2 \Big) e^{xy} 0 \lt x \lt 1 0 \lt y \lt 0.5

    Se convierte a la versión discreta usando diferencias divididas centradas, según se puede mostrar con la gráfica de malla.

    literal b

    EDP eliptica rectangular frontera variable

    literal a


    \frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = \Big( x^2 + y^2 \Big) e^{xy}

    literal c

    Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2

    \frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = \Big( x^2 + y^2 \Big) e^{xy}\frac{\Delta y^2}{\Delta x^2}

    los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

    \lambda= \frac{\Delta y^2}{\Delta x^2} = 1

    se simplifica el coeficiente en λ =1

    u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

    Se agrupan los términosiguales

    u[i-1,j]-4u[i,j]+ u[i+1,j] + u[i,j-1] +u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

    Se desarrollan las iteraciones para tres rombos y se genera el sistema de ecuacioens a resolver.
    para i=1,j=1

    u[0,1]-4u[1,1]+ u[2,1] + u[1,0] +u[1,2] = \Big( x[1]^2 + y[1]^2 \Big) e^{x[1]y[1]} 1-4u[1,1]+ u[2,1] + 1 + \frac{1}{8} = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} -4u[1,1]+ u[2,1] = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} - \frac{1}{8}

    para i=2, j=1

    u[1,1]-4u[2,1]+ u[3,1] + u[2,0] +u[2,2] = \Big( x[2]^2 + y[1]^2 \Big) e^{x[2]y[1]} u[1,1]-4u[2,1]+ u[3,1] + 1 + 0.25 = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} u[1,1]-4u[2,1]+ u[3,1] = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} -1.25

    para i=3, j=1

    u[2,1]-4u[3,1]+ u[4,1] + u[3,0] +u[3,2] = \Big( x[3]^2 + y[1]^2 \Big) e^{x[3]y[1]} u[2,1]-4u[3,1]+ 0.25 + 0 + \frac{3}{8} = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} u[2,1]-4u[3,1] = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} - 0.25 - \frac{3}{8}

    con lo que se puede crear un sistema de ecuaciones y resolver el sistema para cada punto desconocido

    \begin{pmatrix} -4 & 1 & 0 & \Big| & 0.008061807364732415 \\ 1 & -4 & 1 & \Big| & -0.8958911084166168 \\0 & 1 & -4 &\Big| & 0.1288939058881129 \end{pmatrix}

    se obtiene los resultados para:

    u[1,1] = 0.05953113
    u[2,1] = 0.24618634
    u[3,1] = 0.02932311

    >>> import numpy as np
    >>> (0.25**2+0.25**2)*np.exp(0.25*0.25) - 1/8
    0.008061807364732415
    >>> (0.5**2+0.25**2)*np.exp(0.5*0.25) - 1.25
    -0.8958911084166168
    >>> (0.75**2+0.25**2)*np.exp(0.75*0.25) - 0.25 -3/8
    0.1288939058881129
    >>> A=[[-4,1,0],[1,-4,1],[0.0,1.0,-4.0]]
    >>> B = [0.008061807364732415, -0.8958911084166168, 0.1288939058881129]
    >>> np.linalg.solve(A,B)
    array([0.05953113, 0.24618634, 0.02932311])
    
  • s2Eva_2023PAOI_T2 EDO Péndulo vertical amortiguado

    Ejercicio: 2Eva_2023PAOI_T2 EDO Péndulo vertical amortiguado

    literal a

    \frac{d^2 \theta}{dt^2} = -\mu \frac{d\theta}{ dt}-\frac{g}{L}\sin (\theta)

    Se simplifica su forma a:

    \frac{d\theta}{dt}= z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' = -\mu z -\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)

    se usan los valores dados: g = 9.81 m/s2, L = 2 m

    f_t(t,\theta,z) = z g_t(t,\theta,z) = - 0.5 z -\frac{9.81}{2}\sin (\theta)

    y los valores iniciales para la tabla: θ(0) = π/4 rad, θ' (0) = 0 rad/s, se complementan los valores en la medida que se aplica el desarrollo.

    ti θ(ti) θ'(ti)=z
    0 π/4 0
    0.2 0.7161 -0.6583
    0.4 0.5267 -0.1156
    0.6 0.2579 -0.1410
    ...

    literal b

    Iteración 1:  ti = 0 ; yi = π/4 ; zi = 0

    K1y = h * ft(ti,yi,zi) 
        = 0.2*(0) = 0
    K1z = h * gt(ti,yi,zi) 
        = 0.2*(-0.5(0) -(9.81/2)sin (π/4) = -0.6930
            
    K2y = h * ft(ti+h, yi + K1y, zi + K1z)
        = 0.2*(0-0.6930)= -0.1386
    K2z = h * gt(ti+h, yi + K1y, zi + K1z)
        = 0.2*(-0.5(0-0.6930) -(9.81/2)sin(π/4-0) 
        = -0.6237
    
    yi = yi + (K1y+K2y)/2 
       = π/4+ (0+(-0.1386))/2 = 0.7161
    zi = zi + (K1z+K2z)/2 
       = 0+(-0.6930-0.6237)/2 = -0.6583
    ti = ti + h = 0 + 0.2 = 0.2
    
    estimado[i] = [0.2,0.7161,-0.6583]
    

    Iteración 2:  ti = 0.2 ; yi = 0.7161 ; zi = -0.6583

    K1y = h * ft(ti,yi,zi) 
        = 0.2*(-0.6583) = -0.1317
    K1z = h * gt(ti,yi,zi) 
        = 0.2*(- 0.5 ( -0.6583) -(9.81/2)sin (0.7161) 
        = -0.5775
            
    K2y = h * ft(ti+h, yi + K1y, zi + K1z)
        = 0.2*(-0.6583 -0.5775)= -0.2472
    K2z = h * gt(ti+h, yi + K1y, zi + K1z)
        = 0.2*(- 0.5 (-0.6583 -0.5775) -(9.81/2)sin(0.7161-0.1317) 
        = -0.4171
    
    yi = yi + (K1y+K2y)/2 
       = 0.7161 + (-0.1317-0.2472)/2 = 0.5267
    zi = zi + (K1z+K2z)/2 
       = -0.6583+(-0.5775-0.4171)/2 = -0.1156
    ti = ti + h = 0.2 + 0.2 = 0.4
    
    estimado[i] = [0.4,0.5267,-0.1156]
    

    Iteración 3:  ti = 0.4 ; yi = 0.5267 ; zi = -1.156

    K1y = h * ft(ti,yi,zi) 
        = 0.2*(-1.156) = -0.2311
    K1z = h * gt(ti,yi,zi) 
        = 0.2*(- 0.5(-1.156) -(9.81/2)sin (0.5267) 
        = -0.3771
            
    K2y = h * ft(ti+h, yi + K1y, zi + K1z)
        = 0.2*(-1.156 -0.3771)= -0.3065
    K2z = h * gt(ti+h, yi + K1y, zi + K1z)
        = 0.2*(- 0.5 ( -1.156 -0.3771) -(9.81/2)sin(0.5267-0.2311) 
        = -0.1322
    
    yi = yi + (K1y+K2y)/2 
       = 0.5267 + (-0.2311-0.3065)/2 = 0.2579
    zi = zi + (K1z+K2z)/2 
       = -1.156+(-0.3771-0.1322)/2 = -1.410
    ti = ti + h = 0.4 + 0.2 = 0.6
    
    estimado[i] = [0.6,0.2579,-1.410]
    

    literal c

    resultados del algoritmo:

    [ t, 		 y, 	 dyi/dti=z,  K1y,	 K1z,	    K2y,      K2z]
    [[ 0.000e+00  7.854e-01  0.000e+00  0.000e+00  0.000e+00  0.000e+00   0.000e+00]
     [ 2.000e-01  7.161e-01 -6.583e-01  0.000e+00 -6.930e-01 -1.386e-01  -6.237e-01]
     [ 4.000e-01  5.267e-01 -1.156e+00 -1.317e-01 -5.775e-01 -2.472e-01  -4.171e-01]
     [ 6.000e-01  2.579e-01 -1.410e+00 -2.311e-01 -3.771e-01 -3.065e-01  -1.322e-01]
     [ 8.000e-01 -3.508e-02 -1.377e+00 -2.820e-01 -1.089e-01 -3.038e-01    1.756e-01]
    ...
    

    péndulo amortiguado 03

    con h=0.2 se tienen 1/0.2 = 5 tramos por segundo, por lo que para 10 segundo serán 50 tramos. La cantidad de muestras = tramos + 1(valor inicial) = 51

    con lo que se puede usar el algoritmo en EDO Runge-Kutta d2y/dx2

    literal d

    Se observa que la respuesta es oscilante y amortiguada en magnitud como se esperaba según el planteamiento. Con el tiempo se estabilizará en cero.

    Instrucciones en Python

    # 2Eva_2023PAOI_T2 Péndulo vertical amortiguado
    import numpy as np
    import matplotlib.pyplot as plt
    
    def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,7),dtype=float)
        # incluye el punto [x0,y0,z0]
        estimado[0] = [x0,y0,z0,0,0,0,0]
        xi = x0
        yi = y0
        zi = z0
        for i in range(1,tamano,1):
            K1y = h * f(xi,yi,zi)
            K1z = h * g(xi,yi,zi)
            
            K2y = h * f(xi+h, yi + K1y, zi + K1z)
            K2z = h * g(xi+h, yi + K1y, zi + K1z)
    
            yi = yi + (K1y+K2y)/2
            zi = zi + (K1z+K2z)/2
            xi = xi + h
            
            estimado[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
        return(estimado)
    
    # INGRESO theta = y
    g = 9.8
    L  = 2
    ft = lambda t,y,z: z
    gt = lambda t,y,z: -0.5*z +(-g/L)*np.sin(y)
    
    t0 = 0
    y0 = np.pi/4
    z0 = 0
    h = 0.2
    muestras = 51
    
    # PROCEDIMIENTO
    tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)
    
    # SALIDA
    np.set_printoptions(precision=3)
    print(' [ t, \t\t y, \t dyi/dti=z, K1y,\t K1z,\t K2y,\t K2z]')
    print(tabla)
    
    # Grafica
    ti = np.copy(tabla[:,0])
    yi = np.copy(tabla[:,1])
    zi = np.copy(tabla[:,2])
    plt.subplot(121)
    plt.plot(ti,yi)
    plt.grid()
    plt.xlabel('ti')
    plt.title('yi')
    plt.subplot(122)
    plt.plot(ti,zi, color='green')
    plt.xlabel('ti')
    plt.title('dyi/dti')
    plt.grid()
    plt.show()
    
  • s2Eva_2023PAOI_T1 Material para medalla de academia

    Ejercicio: 2Eva_2023PAOI_T1 Material para medalla de academia

    medalla área con integral numérico f(x) = 2-8\Big( \frac{1}{2} - x \Big)^2

    0 \le x \lt 1 g(x) = -\Big( 1-x\Big)\ln \Big( 1- x \Big)

    Para f(x) se usará Simpson de 1/3 que requiere al menos dos tramos para aplicar:

    a. Realice el planteamiento de las ecuaciones para el ejercicio.

    I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]

    b. Describa el criterio usado para determinar el número de tramos usado en cada caso.

    h = \frac{b-a}{2} = \frac{1-0}{2} = 0.5

    c. Desarrolle las expresiones completas del ejercicio para cada función.

    I_{fx}\cong \frac{0.5}{3}[f(0)+4f(0.5) + f(1)] f(0) = 2-8\Big( \frac{1}{2} - (0) \Big)^2 = 0 f(0.5) = 2-8\Big( \frac{1}{2} - (0.5) \Big)^2 = 2 f(1) = 2-8\Big( \frac{1}{2} - (1) \Big)^2 = 0 I_{fx} = \frac{1}{6}[0+4(2) + 0] = \frac{8}{6} = \frac{4}{3} = 1.3333

    cota de error O(h5) = O(0.55)= O(0.03125)

    Para g(x) se usará Simpson de 3/8 que requiere al menos tres tramos para aplicar:

    I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] h = \frac{b-a}{3} = \frac{1-0}{3} = 0.3333 I_{gx}\cong \frac{3(0.3333)}{8}[f(0)+3f(0.3333) +3 f(0.6666)+f(1)] g(0) = -\Big( 1-0\Big)\ln \Big( 1- 0 \Big) = 0 g(0.3333) = -\Big( 1-0.3333\Big)\ln \Big( 1- 0.3333 \Big) = 0.2703 g(0.6666) = -\Big( 1-0.6666\Big)\ln \Big( 1- 0.6666 \Big) = 0.3662 g(0.9999) = -\Big( 1-0.9999\Big)\ln \Big( 1- 0.9999 \Big) = 0

    para la evaluación numérica de 1 se usa un valor muy cercano desplazado con la tolerancia aplicada.

    I_{gx}\cong \frac{3(0.3333)}{8}[0+3(0.2703) + 3(0.3662)+0] = 0.2387

    d. Indique el resultado obtenido para el área requerida y la cota de error
    Area = I_{fx} - I_{gx} = 1.3333 - 0.2387 = 1.0945

    cota de error = O(0.03125) + O(0.00411) = 0.03536

    e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032

    Si el factor de mayor error es de Simpson 1/3, se considera como primera aproximación que:

    cota de error O(h5) = O(0.00032), h = (0.00032)(1/5) = 0.2
    es decir el número de tramos es de al menos (b-a)/tramos = 0.2 , tramos = 5.
    El número de tramos debe ser par en Simpson de 1/3, por lo que se toma el entero mayor tramos=6 y el tamaño de paso recomendado es al menos 1/6. EL error al aplicar 3 veces la formula es 3(O((1/6)5)) = 0.0003858.

    Lo que podría indicar que es necesario al menos dos tramos adicionales con h=1/8 y error O(0,00012) que cumple con el requerimiento.

    Se puede aplicar el mismo criterio para Simpson 3/8 y se combinan los errores para verificar que cumplen con el requerimiento.

    Algoritmo con Python

    Resultados

    Ifx:  1.3333332933359998
    Igx:  0.238779092876627
    Area:  1.094554200459373

    medalla area con integral numerico

    Instrucciones en Python usando las funciones

    # 2Eva_2023PAOI_T1 Material para medalla de academia
    import numpy as np
    import matplotlib.pyplot as plt
    
    def integrasimpson13_fi(xi,fi,tolera = 1e-10):
        ''' sobre muestras de fi para cada xi
            integral con método de Simpson 1/3
            respuesta es np.nan para tramos desiguales,
            no hay suficientes puntos.
        '''
        n = len(xi)
        i = 0
        suma = 0
        while not(i>=(n-2)):
            h = xi[i+1]-xi[i]
            dh = abs(h - (xi[i+2]-xi[i+1]))
            if dh<tolera:# tramos iguales
                unS13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
                suma = suma + unS13
            else:  # tramos desiguales
                suma = 'tramos desiguales'
            i = i + 2
        if i<(n-1): # incompleto, faltan tramos por calcular
            suma = 'tramos incompletos, faltan '
            suma = suma ++str((n-1)-i)+' tramos'
        return(suma)
    
    def integrasimpson38_fi(xi,fi,tolera = 1e-10):
        ''' sobre muestras de fi para cada xi
            integral con método de Simpson 3/8
            respuesta es np.nan para tramos desiguales,
            no hay suficientes puntos.
        '''
        n = len(xi)
        i = 0
        suma = 0
        while not(i>=(n-3)):
            h  = xi[i+1]-xi[i]
            h1 = (xi[i+2]-xi[i+1])
            h2 = (xi[i+3]-xi[i+2])
            dh = abs(h-h1)+abs(h-h2)
            if dh<tolera:# tramos iguales
                unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
                unS38 = (3/8)*h*unS38
                suma = suma + unS38
            else:  # tramos desiguales
                suma = 'tramos desiguales'
            i = i + 3
        if (i+1)<n: # incompleto, tramos por calcular
            suma = 'tramos incompletos, faltan '
            suma = suma +str(n-(i+1))+' tramos'
        return(suma)
    
    # INGRESO
    fx = lambda x: 2-8*(0.5-x)**2
    gx = lambda x: -(1-x)*np.log(1-x)
    a = 0
    b = 1-1e-4
    muestras1 = 2+1
    muestras2 = 3+1
    
    # PROCEDIMIENTO
    xi1 = np.linspace(a,b,muestras1)
    xi2 = np.linspace(a,b,muestras2)
    fi = fx(xi1)
    gi = gx(xi2)
    
    Ifx = integrasimpson13_fi(xi1,fi)
    Igx = integrasimpson38_fi(xi2,gi)
    Area = Ifx - Igx
    
    # SALIDA
    print('Ifx: ', Ifx)
    print('Igx: ', Igx)
    print('Area: ', Area)
    
    plt.plot(xi1,fi,'ob',label='f(x)')
    plt.plot(xi2,gi,'or', label='g(x)')
    plt.grid()
    plt.legend()
    plt.xlabel('xi')
    
    # curvas suave con mas muestras (no en evaluación)
    xi = np.linspace(a,b,51)
    fxi = fx(xi)
    gxi = gx(xi)
    plt.fill_between(xi,fxi,gxi,color='navajowhite')
    plt.plot(xi,fxi,color='blue',linestyle='dotted')
    plt.plot(xi,gxi,color='red',linestyle='dotted')
    
    plt.show()
    
  • s1Eva_2023PAOI_T3 Recoger los restos de sumergible

    ejercicio: 1Eva_2023PAOI_T3 Recoger los restos de sumergible

    literal a

    Para realizar el planteamiento del ejercicio, se realiza la gráfica para observar mejor los puntos.

    recorrido de polinomio

    Se podría observar los puntos y plantear un primer recorrido como el siguiente

    Rov Sumergible Polinomio01a

    los puntos corresponden a los índices j = [0 ,2, 5, 6, 11]:

    # INGRESO
    xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
                   3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
    fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
                   2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])
    
    # selecionando puntos
    j = [0, 2, 5, 6, 11]
    xi = xj[j]
    fi = fj[j]
    

    literal b

    Partiendo del algoritmo de Diferencias divididas de Newton, y añadiendo vectores xj y fj para todos los puntos dados en el ejercicio.

    Tabla Diferencia Dividida
    [['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]', 'F[5]']]
    [[ 0.      0.2478  1.8108  2.569  -1.3163  0.3389 -0.0561  0.    ]
     [ 1.      1.3802  4.7199 -0.7915 -0.1861  0.09    0.      0.    ]
     [ 2.      2.8008  3.5955 -1.2014  0.111   0.      0.      0.    ]
     [ 3.      3.583   2.6558 -0.9928  0.      0.      0.      0.    ]
     [ 4.      4.6794  1.5673  0.      0.      0.      0.      0.    ]]
    dDividida: 
    [ 2.569  -1.3163  0.3389 -0.0561  0.    ]
    polinomio: 
    2.56896856234546*x - 0.0561488275125463*(x - 3.583)*(x - 2.8008)*(x - 1.3802)*(x - 0.2478) 
    + 0.338875729488637*(x - 2.8008)*(x - 1.3802)*(x - 0.2478) 
    - 1.31628089036375*(x - 1.3802)*(x - 0.2478) + 1.17420959025079
    polinomio simplificado: 
    -0.0561488275125463*x**4 + 0.788728905753656*x**3 
    - 3.98331084054791*x**2 + 7.41286537452727*x 
    + 0.206696844067185
    >>>

    usando la tabla de diferencias divididas se plantea la expresión para el polinomio:

    p(x) = 1.8108 + 2.569 (x-0.2478) + + (-1.3163)(x-0.2478)(x-1.3802) + + 0.3389 (x-0.2478)(x-1.3802)(x-2.8008) + + (-0.0561)(x-0.2478)(x-1.3802)(x-2.8008)(x-3.583)

    el algoritmo simplifica la expresión al siguiente polinomio,

    p(x) = -0.0561 x^4 + 0.7887 x^3 - 3.9833 x^2 + 7.4128 x + 0.2066

    literal c

    Se usa el algoritmo Interpolación polinómica de Lagrange con Python para desarrollar el polinomio a partir de los puntos no usados en el literal anterior.

    j = [3,4,7,8,9,10]

    Para evitar oscilaciones grandes, se excluye el punto con índice 1, y se obtiene una opción mostrada:

    Rov Sumergible Polinomio 02a

    con el siguiente resultado:

        valores de fi:  [2.4529 2.3644 4.3657 4.1932 4.6998 4.7536]
    divisores en L(i):  [-1.32578996  0.60430667 -0.02841115  0.02632723 -0.09228243  0.13986517]
    
    Polinomio de Lagrange, expresiones
    -1.85014223337671*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351) 
    + 3.9125829809921*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.4744) 
    - 153.661523787291*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 2.7351)*(x - 2.4744) 
    + 159.272361555669*(x - 4.3757)*(x - 4.2796)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744) 
    - 50.928437685792*(x - 4.3757)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744) 
    + 33.9870187307222*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744)
    
    Polinomio de Lagrange: 
    -9.26814043907703*x**5 + 163.708008152209*x**4 
    - 1143.16428460497*x**3 + 3941.13583467614*x**2 
    - 6701.74628786718*x + 4496.64364271972
    >>> 
    

    Para presentar la parte analítica puede seleccionar menos puntos, se busca mostrar la aplicación del algoritmo, no necesariamente cuan largo puede ser.

    p(x) = + 2.4529\frac{(x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.7351) }{(2.4744 - 4.3757)(2.4744 - 4.2796)(2.4744 - 3.7213)(2.4744 - 3.6627)(2.4744 - 2.7351)} + 2.3644\frac{ (x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.4744) }{(2.7351 - 4.3757)(2.7351 - 4.2796)(2.7351 - 3.7213)(2.7351 - 3.6627)(2.7351 - 2.4744)} +4.3657\frac{(x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 2.7351)(x - 2.4744) }{(3.6627 - 4.3757)(3.6627 - 4.2796)(3.6627 - 3.7213)(3.6627 - 2.7351)(3.6627 - 2.4744)} + 4.1932 \frac{(x - 4.3757)(x - 4.2796)(x - 3.6627)(x - 2.7351)(x - 2.4744) }{(3.7213 - 4.3757)(3.7213 - 4.2796)(3.7213 - 3.6627)(3.7213 - 2.7351)(3.7213 - 2.4744)} +4.6998\frac{(x - 4.3757)(x - 3.7213)(x - 3.6627)(x - 2.7351)(x - 2.4744) }{(4.2796 - 4.3757)(4.2796 - 3.7213)(4.2796 - 3.6627)(4.2796 - 2.7351)(4.2796 - 2.4744)} + 4.7536 \frac{(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.7351)(x - 2.4744)}{(4.3757 - 4.2796)(4.3757 - 3.7213)(4.3757 - 3.6627)(4.3757 - 2.7351)(4.3757 - 2.4744)}

    literal d

    las gráficas se han presentado en el planteamiento para justificar los criterios usados.

    literal e

    Los errores de encuentran como la diferencia de los puntos no usados en el polinomio, y evaluando el polinomio en cada coordenada x.

    Como las propuestas de polinomio pueden ser variadas se obtendrán diferentes respuestas para cada literal.

    Se revisa la aplicación de conceptos en las propuestas.

     

    Algoritmos usados

    literal b. Diferencias divididas de Newton

    # 1Eva_2023PAOI_T3 Recoger los restos de sumergible
    # Polinomio interpolación
    # Diferencias Divididas de Newton
    # Tarea: Verificar tamaño de vectores,
    #        verificar puntos equidistantes en x
    import numpy as np
    import sympy as sym
    import matplotlib.pyplot as plt
    
    # INGRESO
    xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
                   3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
    fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
                   2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])
    
    # selecionando puntos
    j = [0, 2, 5, 6, 11]
    xi = xj[j]
    fi = fj[j]
    
    
    # PROCEDIMIENTO
    
    # Tabla de Diferencias Divididas Avanzadas
    titulo = ['i   ','xi  ','fi  ']
    n = len(xi)
    ki = np.arange(0,n,1)
    tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
    tabla = np.transpose(tabla)
    
    # diferencias divididas vacia
    dfinita = np.zeros(shape=(n,n),dtype=float)
    tabla = np.concatenate((tabla,dfinita), axis=1)
    
    # Calcula tabla, inicia en columna 3
    [n,m] = np.shape(tabla)
    diagonal = n-1
    j = 3
    while (j < m):
        # Añade título para cada columna
        titulo.append('F['+str(j-2)+']')
    
        # cada fila de columna
        i = 0
        paso = j-2 # inicia en 1
        while (i < diagonal):
            denominador = (xi[i+paso]-xi[i])
            numerador = tabla[i+1,j-1]-tabla[i,j-1]
            tabla[i,j] = numerador/denominador
            i = i+1
        diagonal = diagonal - 1
        j = j+1
    
    # POLINOMIO con diferencias Divididas
    # caso: puntos equidistantes en eje x
    dDividida = tabla[0,3:]
    n = len(dfinita)
    
    # expresión del polinomio con Sympy
    x = sym.Symbol('x')
    polinomio = fi[0]
    for j in range(1,n,1):
        factor = dDividida[j-1]
        termino = 1
        for k in range(0,j,1):
            termino = termino*(x-xi[k])
        polinomio = polinomio + termino*factor
    
    # simplifica multiplicando entre (x-xi)
    polisimple = polinomio.expand()
    
    # polinomio para evaluacion numérica
    px = sym.lambdify(x,polisimple)
    
    # Puntos para la gráfica
    muestras = 101
    a = np.min(xi)
    b = np.max(xi)
    pxi = np.linspace(a,b,muestras)
    pfi = px(pxi)
    
    # SALIDA
    np.set_printoptions(precision = 4)
    print('Tabla Diferencia Dividida')
    print([titulo])
    print(tabla)
    print('dDividida: ')
    print(dDividida)
    print('polinomio: ')
    print(polinomio)
    print('polinomio simplificado: ' )
    print(polisimple)
    
    # Gráfica
    plt.plot(xj,fj,'o',color='red')
    plt.plot(xi,fi,'o', label = 'Puntos')
    ##for i in range(0,n,1):
    ##    plt.axvline(xi[i],ls='--', color='yellow')
    
    plt.plot(pxi,pfi, label = 'Polinomio')
    plt.legend()
    plt.grid()
    plt.xlabel('xi')
    plt.ylabel('fi')
    plt.title('Diferencias Divididas - Newton')
    plt.show()
    
    

    literal c. Polinomio de Lagrange

    # 1Eva_2023PAOI_T3 Recoger los restos de sumergible
    # Interpolacion de Lagrange
    # divisoresL solo para mostrar valores
    import numpy as np
    import sympy as sym
    import matplotlib.pyplot as plt
    
    # INGRESO , Datos de prueba
    xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
                   3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
    fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
                   2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])
    
    # selecionando puntos
    #j = [0, 2, 5, 6, 11]
    j = [3,4,7,8,9,10]
    xi = xj[j]
    fi = fj[j]
    
    # PROCEDIMIENTO
    # 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()
    
    # para evaluación numérica
    px = sym.lambdify(x,polisimple)
    
    # Puntos para la gráfica
    muestras = 101
    a = np.min(xi)
    b = np.max(xi)
    pxi = np.linspace(a,b,muestras)
    pfi = px(pxi)
    
    # 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)
    
    # Gráfica
    plt.plot(xj,fj,'o',color='red')
    plt.plot(xi,fi,'o', label = 'Puntos')
    plt.plot(pxi,pfi, label = 'Polinomio')
    plt.legend()
    plt.grid()
    plt.xlabel('xi')
    plt.ylabel('fi')
    plt.title('Interpolación Lagrange')
    plt.show()

     

  • s1Eva_2023PAOI_T2 Productos en combo por subida de precio

    Ejercicio: 1Eva_2023PAOI_T2 Productos en combo por subida de precio

    literal a

    Para el ejercicio solo es posible plantear tres ecuaciones, se muestran datos solo para tres semanas. Se tienen 4 incógnitas, por lo que tendría infinitas soluciones y no se podría resolver.

    500 a + 600 b + 400 c + 90 d = 1660 800 a + 450 b + 300 c + 100 d = 1825 400 a + 300 b + 600 c + 80 d = 1430

    Sin embargo desde el punto de vista práctico se puede usar una variable libre, considerando algunos criterios como:

    - que en la nota se da el precio para la docena de huevos a 2.5 USD
    - que la variable para huevos es la de menor cantidad se puede incurrir en un menor error si el precio ha variado en los últimos días respecto a lo que se podría haber tenido como referencia.

    500 a + 600 b + 400 c = 1660 - 90 d 800 a + 450 b + 300 c = 1825 - 100 d 400 a + 300 b + 600 c = 1430 - 80 d

    que al sustituir con d = 2.5

    500 a + 600 b + 400 c = 1660 - 90 (2.5) =1435 800 a + 450 b + 300 c = 1825 - 100 (2.5) = 1575 400 a + 300 b + 600 c = 1430 - 80 (2.5) = 1230

    Nota: el estudiante puede realizar una justificación diferente para el uso de otra variable libre.

    literal b

    La forma matricial del sistema de ecuaciones se convierte a:

    \begin{pmatrix} 500 & 600 & 400 & \Big| & 1435 \\ 800 & 450 & 300 & \Big| & 1575 \\ 400 & 300 & 600 &\Big| & 1230 \end{pmatrix}

    literal c

    pivoteando la matriz aumentada, se encuentra que para la primera columna se pivotean la fila 1 y 2, por ser 800 el valor de mayor magnitud:

    \begin{pmatrix} 800 & 450 & 300 & \Big| & 1575 \\ 500 & 600 & 400 & \Big| & 1435 \\ 400 & 300 & 600 &\Big| & 1230 \end{pmatrix}

    luego, observando la segunda columna desde el elemento de la diagonal hacia abajo, se observa que el mayor valor en magnitud ya se encuentra en la diagonal. No se requieren intercambios de fila para la tercera columna.

    literal d

    Para el método iterativo de Gauss-Seidel se requiere el vector inicial, tomando las observaciones de la nota. Se da un precio de 50 USD por un quintal métrico de 100Kg, por lo que se estima que el precio por kilogramo ronda 50/100= 0.5 USD. Se indica además que los demás valores se encuentran en el mismo orden de magnitud, por lo que se tiene como vector inicial

    X0 = [0.5, 0.5, 0.5]

    Para un método iterativo se despejan las ecuaciones como:

    a = \frac {1575 - 450 b - 300 c}{800} b = \frac {1435 -500 a - 400 c}{600} c = \frac {1230 - 400 a - 300 b}{600}

    Para el cálculo del error se considera que se busca un precio, lo regular es considerar el centavo, es decir una tolerancia 10-2. Para el desarrollo se considera usar cuatro decimales por si se considera truncar o redondear.

    itera = 0

    a = \frac {1575 - 450 (0.5) - 300 (0.5)}{800} = 1.5 b = \frac {1435 -500(1.5) - 400 (0.5)}{600} = 0.8083 c = \frac {1230 - 400(1.5) - 300 (0.8083)}{600} =0.6458

    errado = max|[1.5-0.5, 0.8083-0.5, 0.6458-0.5]|
    errado = max|[1, 0.3083, 0.1458]| = 1

    itera = 1

    X1 = [1.5, 0.8083, 0.6458]

    a = \frac {1575 - 450 (0.8083) - 300 (0.6458)}{800} = 1.2719 b = \frac {1435 -500(1.2719) - 400 (0.6458)}{600} = 0.9012 c = \frac {1230 - 400 (1.2719) - 300 (0.9012)}{600} = 0.7514

    errado = max|[1.2719-1.5, 0.9012-0.8083, 0.7514-0.6458]|
    errado = max|[-0.2281, 0.0929, 0.1056]| = 0.2281

    el error disminuye entre iteraciones.

    itera = 2

    X2 = [1.2719 , 0.9012, 0.7514]

    a = \frac {1575 - 450 (0.9012) - 300 (0.7514)}{800} = 1.1805 b = \frac {1435 -500 (1.1805) - 400 (0.7514)}{600} = 0.9069 c = \frac {1230 - 400 (1.1805) - 300 (0.9069)}{600} = 0.8095

    errado = max|[ 1.1805-1.2719 , 0.9069-0.9012, 0.8095-0.7514]|
    errado = max|[-0.0914, 0.0057, 0.1056]| = 0.1056

    el error disminuye entre iteraciones.

    literal e

    Si el error entre iteraciones disminuye, se considera que el método converge.

    usando el algoritmo Método de Gauss-Seidel con Python, se tiene como resultado en 13 iteraciones:

    [1.5        0.80833333 0.64583333] 1.0
    [1.271875   0.90121528 0.75147569] 0.2281249999999999
    [1.18001302 0.90733869 0.80965531] 0.09186197916666661
    [1.15475125 0.88960375 0.83536396] 0.025708648304880177
    [1.1550864  0.87218536 0.84384972] 0.01741839593330019
    [1.16170209 0.86101511 0.84502438] 0.011170246538965367
    [1.16754486 0.85536303 0.84395525] 0.005842764350612484
    [1.17112508 0.85309227 0.84270381] 0.0035802211655899807
    [1.17287167 0.85247107 0.84185002] 0.0017465903731879173
    [1.17354127 0.85248226 0.84139802] 0.0006695985845832642
    [1.17370447 0.85264759 0.84120656] 0.00019146597344232852
    [1.17368327 0.8527929  0.84114804] 0.00014530950399282982
    [1.17362348 0.85288174 0.84114348] 8.884049018109685e-05
    respuesta X: 
    [[1.17362348]
     [0.85288174]
     [0.84114348]]
    verificar A.X=B: 
    [[1575.03861029]
     [1434.99817609]
     [1230.        ]]
    iteraciones 13
    >>> 
    

    instrucciones en Python

    Se debe recordar que las matrices y vectores deben ser tipo real (float) .

    # Método de Gauss-Seidel
    # solución de sistemas de ecuaciones
    # por métodos iterativos
    
    import numpy as np
    
    # INGRESO
    A = np.array([[800, 450, 300],
                  [500, 600, 400],
                  [400, 300, 600]], dtype=float)
    
    B = np.array([1575, 1435,1230], dtype=float)
    
    X0  = np.array([0.5,0.5,0.5])
    
    tolera = 0.0001
    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)
        print(X,errado)
        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)
    print('iteraciones',itera)
    
    
  • s1Eva_2023PAOI_T1 Desacople de cohete de dos etapas

    Ejercicio: 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas

    literal a. Planteamiento

    La ecuación a usar según el enunciado es y usando los valores dados es:

    v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt 800 = 1870 \ln\Big(\frac{195000}{195000-2500 t}\Big) - 9.8 t

    con lo que la función para buscar la raíz es:

    f(t) = 1870 \ln\Big(\frac{195000}{195000-2500 t}\Big) - 9.8 t -800

    Literal b. Intervalo de búsqueda

    Para el intervalo de búsqueda se puede usar una gráfica e interpretar el punto a buscar alrededor de 800 m/s. Que de la gráfica se observa que un intervalo alrededor de 35 sería válido para el método de la Bisección. Para otros métodos abiertos, también es posible deducir un punto t0.

    La validación se muestra con la primera iteración al evaluar f(30) que es negativo y f(40) que es de signo positivo.

    Desacople de cohete de dos etapasInstrucciones en Python

    # 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    v = lambda t: 1870*np.log(195000/(195000-2500*t))-9.8*t
    a = 0
    b = 40
    tramos = 51
    
    # PROCEDIMIENTO
    ti = np.linspace(a,b,tramos)
    vi = v(ti)
    
    # SALIDA
    plt.plot(ti,vi)
    plt.xlabel('ti')
    plt.ylabel('vi')
    plt.title('Velocidad vertical vs tiempo')
    plt.grid()
    plt.show()
    

    literal c. Desarrollo con algoritmo de Bisección

    itera =0, intervalo [30,40]

    f(30) = 1870 \ln\Big(\frac{195000}{195000-2500 (30)}\Big) - 9.8 (30) -800 = -186.100 f(40) = 1870 \ln\Big(\frac{195000}{195000-2500 (40)}\Big) - 9.8 (40) -800 = 152.759 c= \frac{a+b}{2}= \frac{30+40}{2} = 35 f(35) = 1870 \ln\Big(\frac{195000}{195000-2500 (35)}\Big) - 9.8 (35) -800 = -29.399

    error = tramo = |40-30| = 10

    como f(a) y f(c) son del mismo signo, el intervalo nuevo será [35,40]

    itera =1 , intervalo [35,40]

    c= \frac{35+40}{2} = 37.5 f(37.5) = 1870 \ln\Big(\frac{195000}{195000-2500 (37.5)}\Big) - 9.8 (37.5) -800 = 58.111

    error = tramo = |40-35| = 5

    como f(c) y f(b) son del mismo signo, el intervalo nuevo será [35,37.5]

    itera =2 , intervalo [35,37.5]

    c= \frac{37.5+35}{2} = 36.25 f(36.25) = 1870 \ln\Big(\frac{195000}{195000-2500 (36.25)}\Big) - 9.8 (36.25) -800 = 13.518

    error = tramo = |37.5-35| = 2.5

    como f(c) y f(b) son del mismo signo, el intervalo nuevo será [35,36.25]

    [ i, a,    c,     b,     f(a),    f(c),   f(b),   tramo]
      1 30.000 35.000 40.000 -186.100 -29.399 152.759 10.000 
      2 35.000 37.500 40.000 -29.399   58.111 152.759  5.000 
      3 35.000 36.250 37.500 -29.399   13.518 58.111   2.500 
      4 35.000 35.625 36.250 -29.399   -8.144 13.518   1.250 
      5 35.625 35.938 36.250 -8.144     2.635 13.518   0.625 
      6 35.625 35.781 35.938 -8.144    -2.767 2.635    0.312 
      7 35.781 35.859 35.938 -2.767    -0.069 2.635    0.156 
      8 35.859 35.898 35.938 -0.069     1.282 2.635    0.078 
    raiz:  35.8984375
    >>>
    

    literal d. tolerancia y errores

    La tolerancia depende de la escala a la que se mide y el instrumento de medición. Si consideramos décimas de segundo la tolerancia será de 10-1. ó 0.1

    Los errores entre iteraciones se muestran en el literal anterior.

    literal e. convergencia

    Los errores en cada iteración disminuye, lo que muestra que el método converge. Luego de 8 iteraciones se encuentra el tiempo ti a usar como 35.8 s.

    Se adjunta el algoritmo de la bisección ajustado para el ejercicio. La gráfica es complementaria a la presentada en el literal b, que puede ser incorporada al mismo algoritmo para la presentación.

    Algoritmo en Python

    # Algoritmo de Bisección
    # 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas
    import numpy as np
    
    # INGRESO
    fx = lambda t: 1870*np.log(195000/(195000-2500*t))-9.8*t -800
    a = 30
    b = 40
    tolera = 0.1
    
    # PROCEDIMIENTO
    tabla = []
    tramo = b-a
    
    fa = fx(a)
    fb = fx(b)
    i = 1
    while (tramo>tolera):
        c = (a+b)/2
        fc = fx(c)
        tabla.append([i,a,c,b,fa,fc,fb,tramo])
        i = i + 1
                     
        cambia = np.sign(fa)*np.sign(fc)
        if (cambia<0):
            b = c
            fb = fc
        else:
            a=c
            fa = fc
        tramo = b-a
    c = (a+b)/2
    fc = fx(c)
    tabla.append([i,a,c,b,fa,fc,fb,tramo])
    tabla = np.array(tabla)
    
    raiz = c
    
    # SALIDA
    np.set_printoptions(precision = 4)
    print('[ i, a, c, b, f(a), f(c), f(b), tramo]')
    # print(tabla)
    
    # Tabla con formato
    n=len(tabla)
    for i in range(0,n,1):
        unafila = tabla[i]
        formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
        unafila = formato.format(*unafila)
        print(unafila)
        
    print('raiz: ',raiz)
    
  • s2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

    Ejercicio: 2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

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

    2Eva_2022PAOII_T3 EDP Parabolica Malla

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

    agrupando variables

    \frac{\Delta t}{b} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = \frac{\Delta t}{b}b\frac{u_{i,j+1}-u_{i,j}}{\Delta t} λ = \frac{\Delta t}{b(\Delta x)^2} λ = \frac{0.002}{2(0.2)^2} =0.025

    como λ<0.5 el método converge.

    \lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] = u[i,j+1]-u[i,j]

    por el método explícito:

    u[i,j+1] =\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] + u[i,j] u[i,j+1] =\lambda u[i+1,j]+(1-2\lambda)u[i,j]+\lambda u[i-1,j]

    iteración i=1, j=0

    u[1,1] =\lambda u[0,0]+(1-2\lambda)u[1,0]+\lambda u[2,0] u[1,1] =0.025(1) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.2\Big)+0.025 \cos \Big( \frac{3π}{2}0.4\Big)

    iteración i=2, j=0

    u[2,1] =\lambda u[1,0]+(1-2\lambda)u[2,0]+\lambda u[3,0] u[2,1] =0.025\cos \Big( \frac{3π}{2}0.2\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.4\Big) +0.025 \cos \Big( \frac{3π}{2}0.6\Big)

    iteración i=3, j=0

    u[3,1] =\lambda u[2,0]+(1-2\lambda)u[3,0]+\lambda u[4,0] u[3,1] =0.025\cos \Big( \frac{3π}{2}0.4\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.6\Big) +0.025 \cos \Big( \frac{3π}{2}0.8\Big)

    iteración i=4, j=0

    u[4,1] =\lambda u[3,0]+(1-2\lambda)u[4,0]+\lambda u[5,0] u[4,1] =0.025\cos \Big( \frac{3π}{2}0.6\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.8\Big) +0.025 (0)

    continuar con las iteraciones en el algoritmo

    Resultados con el algoritmo

    Tabla de resultados
    [[ 1.    1.    1.    1.    1.    1.    1.    1.    1.    1.  ]
     [ 0.59  0.58  0.56  0.55  0.54  0.53  0.53  0.52  0.51  0.5 ]
     [-0.31 -0.3  -0.3  -0.29 -0.28 -0.28 -0.27 -0.27 -0.26 -0.26]
     [-0.95 -0.93 -0.91 -0.89 -0.88 -0.86 -0.84 -0.82 -0.81 -0.79]
     [-0.81 -0.79 -0.78 -0.76 -0.74 -0.73 -0.71 -0.7  -0.68 -0.67]
     [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]]
    

    2Eva2022PAOII_T3 EDP Parabolica 02

    Instrucciones en Python

    # EDP parabólicas d2u/dx2  = K du/dt
    # método explícito,usando diferencias divididas
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    # Valores de frontera
    Ta = 1
    Tb = 0
    #T0 = 25
    fx = lambda x: np.cos(3*np.pi/2*x)
    # longitud en x
    a = 0
    b = 1
    # Constante K
    K = 2
    # Tamaño de paso
    dx = 0.2
    dt = dx/100
    # iteraciones en tiempo
    n = 10
    
    # PROCEDIMIENTO
    # iteraciones en longitud
    xi = np.arange(a,b+dx,dx)
    fi = fx(xi)
    m = len(xi)
    ultimox = m-1
    
    # Resultados en tabla u[x,t]
    u = np.zeros(shape=(m,n), dtype=float)
    
    # valores iniciales de u[:,j]
    j=0
    ultimot = n-1
    u[0,:]= Ta
    u[1:ultimox,j] = fi[1:ultimox]
    u[ultimox,:] = Tb
    
    # factores P,Q,R
    lamb = dt/(K*dx**2)
    P = lamb
    Q = 1 - 2*lamb
    R = lamb
    
    # Calcula U para cada tiempo + dt
    j = 0
    while not(j>=ultimot): # igual con lazo for
        for i in range(1,ultimox,1):
            u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
        j=j+1
    
    # SALIDA
    print('Tabla de resultados')
    np.set_printoptions(precision=2)
    print(u)
    
    # Gráfica
    salto = int(n/10)
    if (salto == 0):
        salto = 1
    for j in range(0,n,salto):
        vector = u[:,j]
        plt.plot(xi,vector)
        plt.plot(xi,vector, '.r')
        
    plt.xlabel('x[i]')
    plt.ylabel('t[j]')
    plt.title('Solución EDP parabólica')
    plt.show()
    
  • s2Eva_2022PAOII_T2 EDO - población de protestantes en una sociedad

    Ejercicio: 2Eva_2022PAOII_T2 EDO - población de protestantes en una sociedad

    \frac{\delta}{\delta t}x(t) = b x(t) - d (x(t))^2 \frac{\delta}{\delta t}y(t) = b y(t) - d (y(t))^2 +r b (x(t)-y(t))

    literal a

    simplificando la nomenclatura

    x' = b x - d x^2 y' = b y - d y^2 +r b (x-y)

    sustituyendo constantes, y considerando x(0)=1 ; y(0)=0.01 ; h=0.5

    x' = 0.02 x - 0.015 x^2 y' = 0.02 y - 0.015 y^2 +0.1(0.02) (x-y)

    el planteamiento de Runge Kutta se hace junto a la primera iteración, además de encontrarse en las instrucciones con Python.

    literal b

    Se describen 3 iteraciones usando los resultados de la tabla con el algoritmo, para mostrar la comprensión del algoritmo.

    t = 0
    K1x = 0.5 (0.02 (1) - 0.015 (1)2 = 0.0025
    K1y = 0.5(0.02 (0.01) - 0.015 (0.01)2 +0.1(0.02) (1-0.01)= 0.001089

    K2x = 0.5 (0.02 (1+0.0025) - 0.015 (1+0.0025)2= 0.00248
    K2y = 0.5(0.02 (0.01+0.00108) - 0.015 (0.01+0.00108)2+0.1(0.02) ((1+0.0025)-(0.01+0.00108)) = 0.001101

    x1 = 1 + (1/2)(0.0025+0.00248) = 1.0025
    y1 = 0.01 + (1/2)(0.001089+0.001101) = 0.01109
    t1 = 0 + 0.5 =0.5

    t=0.5
    K1x = 0.5 (0.02 (1.0025) - 0.015 (1.0025)2 = 0.002487
    K1y = 0.5(0.02 (0.01109) - 0.015 (0.01109)2 +0.1(0.02) (1.0025-0.01109)= 0.001101

    K2x = 0.5 (0.02 (1.0025+ 0.002487) - 0.015 (1.0025+ 0.002487)2= 0.002474
    K2y = 0.5(0.02 (0.01109+0.001101) - 0.015 (0.01109+0.001101)2+0.1(0.02) ((1.0025+ 0.002487)-(0.01109+0.001101)) = 0.001113

    x2 = 1.0025 + (1/2)(0.002487+0.002474) = 1.0050
    y2 = 0.01109 + (1/2)(0.001101+0.001113) = 0.01220
    t2 = 0.5 + 0.5 = 1

    t=1
    K1x = 0.5 (0.02 (1.0050) - 0.015 (1.0050)2 = 0.002474
    K1y = 0.5(0.02 (0.01220) - 0.015 (0.01220)2 +0.1(0.02) (1.0050-0.01220)= 0.001113

    K2x = 0.5 (0.02 (1.0050+ 0.002474) - 0.015 (1.0050+ 0.002474)2= 0.002462
    K2y = 0.5(0.02 (0.01220+0.001113) - 0.015 (0.01220+0.001113)2+0.1(0.02) ((1.0050+ 0.002474)-(0.01220+0.001113)) = 0.001126

    x3 = 1.0050 + (1/2)(0.002474+0.002462) = 1.0074
    y3 = 0.01220 + (1/2)(0.001113+0.001126) = 0.01332
    t3 = 1 + 0.5 = 1.5

    Resultado con el algoritmo

    Para obtener los datos de las iteraciones, primero se ejecuta el algoritmo para pocas iteraciones.
    Para la pregunta sobre 200 años, se incrementa las iteraciones a 2 por año y las condiciones iniciales, es decir 401 muestras.

     [ ti, xi, yi]
     [ ti, xi, yi]
    [[0.0000e+00 1.0000e+00 1.0000e-02 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
     [5.0000e-01 1.0025e+00 1.1095e-02 2.5000e-03 1.0892e-03 2.4875e-03 1.1014e-03]
     [1.0000e+00 1.0050e+00 1.2203e-02 2.4875e-03 1.1014e-03 2.4749e-03 1.1136e-03]
     [1.5000e+00 1.0074e+00 1.3323e-02 2.4749e-03 1.1137e-03 2.4623e-03 1.1260e-03]
     [2.0000e+00 1.0099e+00 1.4455e-02 2.4624e-03 1.1260e-03 2.4497e-03 1.1384e-03]
     [2.5000e+00 1.0123e+00 1.5600e-02 2.4498e-03 1.1384e-03 2.4371e-03 1.1509e-03]
     [3.0000e+00 1.0148e+00 1.6757e-02 2.4371e-03 1.1509e-03 2.4245e-03 1.1634e-03]
     [3.5000e+00 1.0172e+00 1.7926e-02 2.4245e-03 1.1635e-03 2.4118e-03 1.1761e-03]
     [4.0000e+00 1.0196e+00 1.9109e-02 2.4118e-03 1.1761e-03 2.3991e-03 1.1888e-03]
    ...
     [1.9950e+02 1.3252e+00 1.1561e+00 ... 1.7202e-03 8.1217e-05 1.7059e-03]
     [2.0000e+02 1.3252e+00 1.1578e+00 ... 1.7060e-03 8.0418e-05 1.6918e-03]
     [2.0050e+02 1.3253e+00 1.1595e+00 ... 1.6919e-03 7.9628e-05 1.6778e-03]]
    >>> 
    

    Observación: La población identificada como protestante, continua creciendo, mientras que la proporción de "conformistas" se reduce según los parámetros indicados en el ejercicio. Los valores de natalidad y defunción cambian con el tiempo mucho más en años por otras variables, por lo que se deben realizar ajustes si se pretende extender el modelo.

    2Eva2022PAOII_T2 poblacion protestantes
    Instrucciones en Python

    # Modelo predador-presa de Lotka-Volterra
    # Sistemas EDO con Runge Kutta de 2do Orden
    import numpy as np
    
    def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
        tamano = muestras +1
        tabla = np.zeros(shape=(tamano,7),dtype=float)
        tabla[0] = [t0,x0,y0,0,0,0,0]
        ti = t0
        xi = x0
        yi = y0
        for i in range(1,tamano,1):
            K1x = h * f(ti,xi,yi)
            K1y = h * g(ti,xi,yi)
            
            K2x = h * f(ti+h, xi + K1x, yi+K1y)
            K2y = h * g(ti+h, xi + K1x, yi+K1y)
    
            xi = xi + (1/2)*(K1x+K2x)
            yi = yi + (1/2)*(K1y+K2y)
            ti = ti + h
            
            tabla[i] = [ti,xi,yi,K1x,K1y,K2x,K2y]
        tabla = np.array(tabla)
        return(tabla)
    
    # PROGRAMA ------------------
    
    # INGRESO
    # Parámetros de las ecuaciones
    b = 0.02
    d = 0.015
    r = 0.1
    
    # Ecuaciones
    f = lambda t,x,y : (b-d*x)*x
    g = lambda t,x,y : (b-d*y)*y + r*b*(x-y)
    
    # Condiciones iniciales
    t0 = 0
    x0 = 1
    y0 = 0.01
    
    # parámetros del algoritmo
    h = 0.5
    muestras = 401
    
    # PROCEDIMIENTO
    tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
    ti = tabla[:,0]
    xi = tabla[:,1]
    yi = tabla[:,2]
    
    # SALIDA
    np.set_printoptions(precision=6)
    print(' [ ti, xi, yi, K1x, K1y, K2x, K2y]')
    print(tabla)
    
    # Grafica tiempos vs población
    import matplotlib.pyplot as plt
    
    plt.plot(ti,xi, label='xi poblacion')
    plt.plot(ti,yi, label='yi protestante')
    
    plt.title('población y protestantes')
    plt.xlabel('t años')
    plt.ylabel('población')
    plt.legend()
    plt.grid()
    plt.show()
    
    # gráfica xi vs yi
    plt.plot(xi,yi)
    
    plt.title('población y protestantes [xi,yi]')
    plt.xlabel('x población')
    plt.ylabel('y protestantes')
    plt.grid()
    plt.show()