Categoría: Soluciones

Ejercicios resueltos de examen

  • s3Eva_IT2018_T1 Intersección de dos círculos

    Ejercicio: 3Eva_IT2018_T1 Intersección de dos círculos

    Para la solución se presentan dos secciones:

    1. Solución particular de intersección de círculos

    2. Solución General de intersección de círculos

    _


    1. Solución Particular de intersección de círculos

    La solución particular se enfoca en el enunciado del ejercicio presentado

    Literal a

    Se grafica las funciones usando Python, para encontrar el rango de búsqueda de raíces.

    intersecta circulos 01

    De la gráfica se usa el 'zoom' y se puede aproximar los valores para la intersección de las curvas estimando raíces en x=1.80 y x=3.56

    Desarrollo numérico

    Se usan las ecuaciones para encontrar la diferencia entre las funciones.

    (x-4)^2 + (y-4)^2 = 5 x^2 + y^2 = 16

    Se despeja la variable y para la primera ecuación:

    (y-4)^2 = 5 - (x-4)^2 y-4 = \sqrt{5 - (x-4)^2} f(x) = y = \sqrt{5 - (x-4)^2} + 4

    la segunda ecuación se transforma en

    x^2 + y^2 = 16 y^2 = 16 - x^2 g(x) = y = \sqrt{16 - x^2}

    La intersección se obtiene restando las ecuaciones, para f(x) se usa la parte inferior del circulo y para g(x) la parte superior de circulo.

    Para buscar las raíces se analiza en el rango de existencia entre las dos funciones:

    [-4,4]\text{ y } [4 -\sqrt{5} ,4 + \sqrt{5}] [-4,4] \text{ y } [1.7639 , 6.2360]

    por lo que la diferencia existe en el rango:

    [1.7639 ,4] \text{diferencia}(x) = f(x)-g(x)

    que es el que se usa para el literal b

    intersecta circulos 02 diferencia


    Literal b

    Las ecuaciones para la diferencia entre las funciones son :

    f_{2} (x) = -\sqrt{5-(x-4)^2}+4 g_{1} (x) = \sqrt{16-x^2}

    Para el método de Newton-Raphson se requieren las derivadas:

    \frac{d f_2}{dx} = \frac{x-4}{ \sqrt{5-(x-4)^2} } \frac{d g_{1}}{dx} = \frac{-x}{ \sqrt{16-x^2} }

    por lo que:

    \frac{d \text{diferencia}}{dx} = \frac{d f_{2}}{dx} - \frac{d g_{1}}{dx}

    Usando el algoritmo con Python se obtienen las raices:

     usando Newton-Raphson
    raices en:  1.80582463574 3.56917099898

    Desarrollo en Python:

    El desarrollo se realiza por partes, en el mismo orden del planteamiento de  los literales

    # 3ra Evaluación I Término 2018
    # Tema 1. Intersección de círculos
    import numpy as np
    import matplotlib.pyplot as plt
    
    # literal a
    
    fx1 = lambda x: np.sqrt(5-(x-4)**2)+4
    fx2 = lambda x: -np.sqrt(5-(x-4)**2)+4
    gx1 = lambda x: np.sqrt(16-x**2)
    gx2 = lambda x: -np.sqrt(16-x**2)
    
    # Rango inicial de análisis (visual)
    a = -5; b = 7
    muestras = 501
    
    # PROCEDIMIENTO
    # Evalua los puntos en el rango
    xi = np.linspace(a,b,muestras)
    fx1i = fx1(xi)
    fx2i = fx2(xi)
    gx1i = gx1(xi)
    gx2i = gx2(xi)
    
    # SALIDA - Gráfica
    plt.plot(xi,fx1i)
    plt.plot(xi,fx2i)
    plt.plot(xi,gx1i)
    plt.plot(xi,gx2i)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Intersección de círculos')
    plt.grid()
    plt.show()
    
    # GRAFICAR las diferencias
    a = 4 - np.sqrt(5)
    b = 4 + np.sqrt(5)
    # PROCEDIMIENTO
    xi = np.linspace(a,b,muestras)
    diferencia = fx2(xi) - gx1(xi)
    # GRAFICA
    plt.plot(xi,diferencia)
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('diferencia entre círculos')
    plt.grid()
    plt.show()
    
    # literal b -----------------------
    def newton_raphson(funcionx, fxderiva, xi, tolera):
        # funciónx y fxderiva en forma numérica
        # xi es el punto inicial de búsqueda
        tramo = abs(2*tolera)
        while (tramo>=tolera):
            xnuevo = xi - funcionx(xi)/fxderiva(xi)
            tramo = abs(xnuevo-xi)
            xi = xnuevo
        return(xi)
    
    funcionx = lambda x: fx2(x) - gx1(x)
    fxderiva = lambda x: (x-4)/np.sqrt(5-(x-4)**2)+x/np.sqrt(16-x**2)
    
    tolera = 0.001
    xi1 = a + tolera
    xi2 = 3.5
    
    raiz1 = newton_raphson(funcionx, fxderiva, xi1, tolera)
    raiz2 = newton_raphson(funcionx, fxderiva, xi2, tolera)
    
    # SALIDA
    print('\n usando Newton-Raphson')
    print('raices en: ', raiz1,raiz2)
    

    _


    2. Solución General de intersección de círculos

    Una solución más general de la intersección de círculos, considerada como para una actividad de mayor duración, revisa previamente si existe un cruce de áreas entre los dos círculos y estima el intervalo donde se encuentran las raíces [xa,xb].

    De existir esa posibilidad, con el intervalo anterior  [xa,xb] busca por un método de búsqueda de raíces las coordenadas de la intersección de las circunferencias.

    2.1 Buscar cruce de áreas entre dos círculos

    El cruce de áreas entre dos círculos se determina comparando si la distancia entre la suma de los radios es mayor o igual a la distancia entre los centros de los círculos.

    De cumplirse la condición anterior, es posible encontrar las intersecciones de los círculos. El valor xa se obtiene como el mayor entre los límites x hacia la izquierda de cada círculo, mientras que xb se obtiene como el límite x hacia la derecha entre los círculos.

    intersecta Circulos 00

    Lo siguiente que hay que reconocer es cuál de las partes (superior e inferior) de cada círculo es necesario usar para encontrar las intersecciones. Esta sección es necesaria puesto que la fórmula que describe el círculo contiene una raiz cuadrada que puede se positiva o negativa, generando dos segmentos en cada círculo.

    Por ejemplo, partiendo de la fórmula general de un círculo con centro en [x1,y1] y radio r1:

    (x-x_1)^2 + (y-y_1)^2 = r_1^2 (y-y_1)^2 = r_1^2 - (x-x_1)^2 \sqrt{(y-y_1)^2} = \sqrt{r_1^2 - (x-x_1)^2} y = \sqrt{r_1^2 - (x-x_1)^2} + y_1

    Con lo que se muestra la necesidad de identificar para cada círculo el sector arriba y abajo que interviene para encontrar las intersecciones. El orden del sector se establece con las posibles combinaciones de:

    tabla de signos en raíz cuadrada para círculo
    círculo 2 abajo círculo2 arriba
    círculo 1 abajo [-1,-1] [-1,1]
    círculo 1 arriba [ 1,-1] [ 1,1]

    El uso de cada combinación se estrablece en el vector de 1 y 0 con el siguiente orden:

    sector = [ abajo1*abajo2,  abajo1*arriba2,
              arriba1*abajo2, arriba1*arriba2]

    las instrucciones en Python para lo descrito se muestran como una función:

    import numpy as np
    import scipy.optimize as sp
    def cruce2circulos(x1,y1,r1,x2,y2,r2):
        ''' Revisa intervalo de area de cruce
            entre dos círculos de centro y radio
            x1,y1,r1 // x2,y2,r2
        '''
        intersecta = []
        dx = x2 - x1
        dy = y2 - y1
        d_centros = np.sqrt(dx**2 + dy**2)
        d_cruce   = r2 + r1
        
        # los circulos se cruzan o tocan
        if d_cruce >= d_centros:
    
            # intervalos de cruce
            xa = np.max([x1-r1,x2-r2])
            xb = np.min([x1+r1,x2+r2])
            ya = np.max([y1-r1,y2-r2])
            yb = np.min([y1+r1,y2+r1])
            
            # cada circulo arriba, abajo
            abajo1 = 0 ; arriba1 = 0
            abajo2 = 0 ; arriba2 = 0
            if ya<=y1:
                abajo1  = 1
            if yb>=y1:
                arriba1 = 1
            if ya<=y2:
                abajo2  = 1
            if yb>=y2:
                arriba2 = 1
            sector  = [ abajo1*abajo2, abajo1*arriba2,
                       arriba1*abajo2, arriba1*arriba2]
            uncruce = [xa,xb,ya,yb,sector]
        return(uncruce)
    

    El resultado para los círculos del ejercicio son:

    >>> x1=4; y1=4; r1=np.sqrt(5)
    >>> x2=0; y2=0; r2=np.sqrt(16)
    >>> uncruce = cruce2circulos(x1,y1,r1,x2,y2,r2)
    >>> uncruce
    [1.7639320225002102, 4.0, 
     1.7639320225002102, 2.23606797749979, 
    [0, 1, 0, 0]]
    >>> 
    

    2.2 Raíces como coordenadas de intersección entre dos círculos

    Las coordenadas de intersección entre dos círculos se obtienen aplicando un método de búsqueda de raíces. Por ejemplo bisección, que para esta parte se usa el algoritmo de SciPy con la instrucción sp.bisect(fx,xa,xb,xtol=2e-12).

    Para el caso más general, donde existen dos raíces que buscar, se divide el intervalo de búsqueda [xa,xb] en dos medios segmentos [xa,xc] y [xc,xb]. Se aplica un método de búsqueda de raíces para cada subintervalo. Para minimizar errores de truncamiento, en cada búsqueda de desplaza dx/10 cada xc hacia el lado que amplia el subintervalo de búsqueda.

    intersecta Circulos 01

    Para el caso donde los círculos solo tienen un punto de contacto, se realiza una revisión considerando que el intervalo de búsqueda podría ser menor al valor de tolerancia del radio.

    Por ejemplo, cuando la linea que une los centros de los círculos resulta paralelos al eje de las x,  adicionalmente se topan en un solo punto, el algoritmo anterior indica que se usan todos los sectores de los círculos, dando como resultado cuatro raíces iguales. El caso se corrige realizando la parte de sectores solo cuando la distancia entre [xa,xb] es mayor a cero.

    El resultado se presenta como los vectores raizx y raizy.

    Las instrucciones en Python para esta sección se describen a continuación:

    def raices2circulos(x1,y1,r1,x2,y2,r2,tolera=2e-12):
        ''' busca las intersección entre 2 circulos
            de centro y radio: x1,y1,r1 || x2,y2,r2
            revisa con cruce2circulos()
        '''
        uncruce = cruce2circulos(x1,y1,r1,x2,y2,r2)
        raizx = []; raizy = []
    
        # si hay cruce de circulos
        if len(uncruce)>0:
            sectores = [[-1,-1],[-1,1], 
                        [ 1,-1],[ 1,1]]
            [xa,xb,ya,yb,sector] = uncruce
            xc = (xa+xb)/2
            dx = np.abs(xb-xa)
            dy = np.abs(yb-ya)
            k = 1    # se tocan en un punto
            if dx>0: # se tocan en mas de un punto
                k = len(sector)
            for j in range(0,k,1):
                if sector[j]==1:
                    s1 = sectores[j][0]
                    s2 = sectores[j][1]
                    fx1 = lambda x: s1*np.sqrt(r1**2-(x-x1)**2)+y1
                    fx2 = lambda x: s2*np.sqrt(r2**2-(x-x2)**2)+y2
                    fx  = lambda x: fx1(x)-fx2(x)
                    fa = fx(xa)
                    fb = fx(xb)
                    raiz1 = np.nan
                    raiz2 = np.nan
                    
                    # intervalo/2 izquierda
                    xc = xc + dx/10
                    fc = fx(xc)
                    cambio = np.sign(fa)*np.sign(fc)
                    if cambio<0:
                        raiz1 = sp.bisect(fx,xa,xc,xtol=tolera)
                        
                    # intervalo/2 derecha
                    xc = xc - 2*dx/10
                    fc = fx(xc)
                    cambio = np.sign(fc)*np.sign(fb)
                    if cambio<0:
                        raiz2 = sp.bisect(fx,xc,xb,xtol=tolera)
                        
                    # si hay contacto en un borde
                    if dx<tolera*r1 and dy>0:
                        raiz1 = xa
                    if dy<tolera*r1 and dx>0:
                        raiz1 = x1
                        
                    # Añade si existe raiz
                    if not(np.isnan(raiz1)):
                        raizx.append(raiz1)
                        raizy.append(fx1(raiz1))
                    if not(np.isnan(raiz2)):
                        raizx.append(raiz2)
                        raizy.append(fx1(raiz2))
            raices = [raizx,raizy]
        return(raices)
    

    El resultado del algoritmo para el ejercicio es:

    >>> raices = raices2circulos(x1,y1,r1,x2,y2,r2,tolera=2e-12)
    >>> raices
    [[1.805829001269906, 3.569170998730207],
     [3.569170998734088, 1.8058290012706681]]
    >>>
  • s2Eva_IT2018_T3 EDP Eliptica

    Ejercicio: 2Eva_IT2018_T3 EDP Eliptica

    Generar las ecuaciones a resolver usando diferencias finitas divididas centradas:

    \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = 2(x^2+y^2)

    por facilidad se sustituye también en la forma discreta de la ecuación:

    f[i,j] = f(x_i,y_j) = 2(x_{i} ^2+y_{j} ^2)
    \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} = f[i,j]
    \frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j]\Big) + + u[i,j-1]-2u[i,j]+u[i,j+1] = \Delta y^2 f[i,j]

    dado que hx = hy = 1/3

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

    La última ecuación puede ser usada de forma iterativa, para lo cual hay que definir los valores iniciales de la matriz u.

    Al conocer el rango de operación para los ejes x, y, hx, hy se realizan los cálculos para:

    1. Evaluar los valores para cada eje x[i], y[j]

    x[i]:
    [ 0.    0.33  0.67  1.  ]
    y[j]:
    [ 0.    0.33  0.67  1.  ]
    

    2. Evaluar en cada punto generando una matriz f(i,j):

    f[i,j]:
    [[ 0.    0.22  0.89  2.  ]
     [ 0.22  0.44  1.11  2.22]
     [ 0.89  1.11  1.78  2.89]
     [ 2.    2.22  2.89  4.  ]]
    

    3. Se evaluan las funciones indicadas para la frontera y se tiene la matriz inicial para u:

    matriz inicial u[i,j]:
    [[ 1.    1.33  1.67  2.  ]
     [ 1.33  0.    0.    2.44]
     [ 1.67  0.    0.    3.11]
     [ 2.    2.44  3.11  4.  ]]
    

    con lo que se puede trabajar cada punto i,j de forma iterativa, teniendo como resultado para la matriz u:

    resultado para u, iterando: 
    converge =  1
    [[ 1.    1.33  1.67  2.  ]
     [ 1.33  1.68  2.05  2.44]
     [ 1.67  2.05  2.53  3.11]
     [ 2.    2.44  3.11  4.  ]]
    

    La gráfica usando una mayor resolución para tener una idea de la solución:


    Los resultados se obtienen usando las siguientes instrucciones:

    # 2da Evaluación I Término 2018
    # Tema 3. EDP Eliptica
    import numpy as np
    
    # INGRESO
    # ejes x,y
    x0 = 0 ; xn = 1 ; hx = (1/3)# (1/3)/10
    y0 = 0 ; yn = 1 ; hy = (1/3) # (1/3)/10
    # Fronteras
    fux0 = lambda x: x+1
    fu0y = lambda y: y+1
    fux1 = lambda x: x**2 + x + 2
    fu1y = lambda y: y**2 + y + 2
    
    fxy = lambda x,y: 2*(x**2+y**2)
    
    # PROCEDIMIENTO
    xi = np.arange(x0,xn+hx,hx)
    yj = np.arange(y0,yn+hy,hy)
    n = len(xi)
    m = len(yj)
    # funcion f[xi,yi]
    fij = np.zeros(shape=(n,m), dtype = float)
    for i in range(0,n,1):
        for j in range(0,m,1):
            fij[i,j]=fxy(xi[i],yj[j])
    # matriz inicial u[i,j]
    u = np.zeros(shape=(n,m), dtype = float)
    u[:,0] = fux0(xi)
    u[0,:] = fu0y(yj)
    u[:,m-1] = fux1(xi)
    u[n-1,:] = fu1y(yj)
    
    uinicial = u.copy()
    
    # Calcular de forma iterativa
    maxitera = 100
    tolera = 0.0001
    # valor inicial de iteración
    promedio = (np.max(u)+np.min(u))/2
    u[1:n-1,1:m-1] = promedio
    # iterar puntos interiores
    itera = 0
    converge = 0
    erroru = 2*tolera # para calcular al menos una matriz
    while not(erroru=maxitera):
        itera = itera +1
        nueva = np.copy(u)
        for i in range(1,n-1):
            for j in range(1,m-1):
                u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-(hy**2)*fij[i,j])/4
        diferencia = nueva-u
        erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1
    
    # SALIDA
    np.set_printoptions(precision=2)
    print('x[i]:')
    print(xi)
    print('y[j]:')
    print(yj)
    print('f[i,j]:')
    print(fij)
    print('matriz inicial u[i,j]:')
    print(uinicial)
    print('resultado para u, iterando: ')
    print('converge = ', converge)
    print('iteraciones = ', itera)
    print(u)
    

    para obtener la gráfica se debe añadir:

    # Gráfica
    import matplotlib.pyplot as plt
    from matplotlib import cm
    from mpl_toolkits.mplot3d import Axes3D
    
    X, Y = np.meshgrid(xi, yj)
    figura = plt.figure()
    ax = Axes3D(figura)
    U = np.transpose(u) # ajuste de índices fila es x
    ax.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds)
    plt.title('EDP elíptica')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.show()
    
  • s2Eva_IT2018_T2 Deducir Simpson 1/3

    Ejercicio: 2Eva_IT2018_T2 Deducir Simpson 1/3

    Para el problema, se usan los puntos:  [a,f(a)], [b,f(b)] y [c,f(c)]
    por donde pasa la curva f(x) aproximada a un polinomio de grado 2, f(x) \approx p(x)

    \int_a^b f(x) dx \approx \int_a^b p(x) dx p(x) = L_a f(a) + L_c f(c) + L_b f(b) \int_a^b p(x) dx = \int_a^b \Big[ L_a f(a) +L_c f(c) + L_b f(b) \Big] dx = = \int_a^b L_a f(a) dx + \int_a^b L_c f(c) dx + \int_a^b L_b f(b) dx \int_a^b p(x) dx = I_1 + I_2 + I_3

    Como referencia se usa la gráfica para relacionar a, b, c y h:

    regla simpson 01


    Primer Integral

    Para el primer integral I_1= \int_a^b L_a f(a) dx se tiene que:

    L_a = \frac{(x-b)(x-c)}{(a-b)(a-c)} = \frac{(x-b)(x-c)}{(-2h)(-h)} L_a = \frac{(x-b)(x-c)}{2h^2}

    se convierte en:

    I_1 = \int_a^b \frac{(x-b)(x-c)}{2h^2} f(a) dx = \frac{f(a)}{2h^2} \int_a^b (x-b)(x-c)dx

    Para dejar la parte del integral en función de h, a y b, teniendo que c está en la mitad de [a,b], es decir c = (a+b)/2 , se usa:

    u = x-a

    por lo que \frac{du}{dx}=1 y du = dx

    x-c = (u+a) - \frac{a+b}{2} = u+ \frac{a-b}{2} = u - \frac{b-a}{2} x-c = u-h x-b = (u+a)-b = u-2\Big(\frac{b-a}{2}\Big) = u-2h

    Se actualiza el integral de x entre [a,b]  usando u = x-a, que se convierte el rango para u en [0, b-a] que es lo mismo que [0,2h]

    \int_a^b (x-b)(x-c)dx = \int_0^{2h} (u-2h)(u-h)du = = \int_0^{2h} \Big( u^2 - 2hu - uh + 2h^2 \Big) du = \int_0^{2h} \Big( u^2 - 3hu + 2h^2 \Big) du = \frac{u^3}{3}- 3h\frac{u^2}{2}+ 2h^2u \Big|_0^{2h} = \frac{(2h)^3}{3}- 3h\frac{(2h)^2}{2} + 2h^2(2h) -(0-0+0) = \frac{8h^3}{3}- 6h^3 + 4h^3 =\frac{8h^3}{3}- 2h^3 = \frac{2h^3}{3}

    resultado que se usa en I1

    I_1= \frac{f(a)}{2h^2}\frac{2h^3}{3} =\frac{h}{3} f(a)

    que es el primer término de la fórmula general de Simpson 1/3


    Segundo Integral

    Para el Segundo integral I_2= \int_a^b L_c f(c) dx se tiene que:

    L_c = \frac{(x-a)(x-b)}{(c-a)(c-b)} = \frac{(x-a)(x-b)}{(h)(-h)} L_c = \frac{(x-a)(x-b)}{-h^2}

    se convierte en:

    I_2 = \frac{f(c)}{-h^2} \int_a^b (x-a)(x-b) dx = \frac{f(c)}{-h^2} \int_0^{2h} (u)(u-2h) du

    siendo:

    \int_0^{2h}(u^2-2hu)du=\Big(\frac{u^3}{3}-2h\frac{u^2}{2}\Big)\Big|_0^{2h} =\frac{(2h)^3}{3}-h(2h)^2-(0-0) =\frac{8h^3}{3}-4h^3 = -\frac{4h^3}{3}

    usando en I2

    I_2 = \frac{f(c)}{-h^2}\Big(-\frac{4h^3}{3}) = \frac{h}{3}4f(c)

    Tarea: Continuar las operaciones para y tercer integral para llegar a la fórmula general de Simpson 1/3:

    I = \frac{h}{3} \Big( f(a)+4f(c) + f(b) \Big)
  • s2Eva_IT2018_T4 Dragado acceso marítimo

    Ejercicio: 2Eva_IT2018_T4 Dragado acceso marítimo

    a) La matriz para remover sedimentos se determina como la diferencia entre la profundidad y la matriz de batimetría.

    Considere el signo de la profundidad para obtener el resultado:

    matriz remover sedimentos: 
    [[ 4.21  0.    0.96  0.    3.76  3.09]
     [ 2.15  0.11  2.05  3.77  0.    3.07]
     [ 0.    1.14  1.65  0.    1.62  1.35]
     [ 3.7   0.    0.59  2.33  0.    4.23]
     [ 0.    1.38  3.53  4.49  1.98  1.4 ]
     [ 0.    0.77  0.32  1.06  4.24  3.54]]
    

    se obtiene con la instrucciones:

    # 2da Evaluación I Término 2018
    # Tema 4. canal acceso a Puertos de Guayaquil
    import numpy as np
    
    # INGRESO
    profcanal = 11
    
    xi = np.array([ 0.,  50., 100., 150., 200., 250.])
    yi = np.array([ 0., 100., 200., 300., 400., 500.])
    
    batimetria = [[ -6.79,-12.03,-10.04,-11.60, -7.24,-7.91],
                  [ -8.85,-10.89, -8.95, -7.23,-11.42,-7.93],
                  [-11.90, -9.86, -9.35,-12.05, -9.38,-9.65],
                  [ -7.30,-11.55,-10.41, -8.67,-11.84,-6.77],
                  [-12.17, -9.62, -7.47, -6.51, -9.02,-9.60],
                  [-11.90,-10.23,-10.68, -9.94, -6.76,-7.46]]
    
    batimetria = np.array(batimetria)
    # PROCEDIMIENTO
    [n,m] = np.shape(batimetria)
    
    # Matriz remover sedimentos
    remover = batimetria + profcanal
    for i in range(0,n,1):
        for j in range(0,m,1):
            if remover[i,j]<0:
                remover[i,j]=0
    # SALIDA
    print('matriz remover sedimentos: ')
    print(remover)
    

    b) el volumen se calcula usando el algoritmo de Simpson primero por un eje, y luego con el resultado se continúa con el otro eje,

    Considere que existen 6 puntos, o 5 tramos integrar en cada eje.

    • Al usar Simpson de 1/3 que usan tramos pares, faltaría integrar el último tramo.
    • En el caso de Simpson de 3/8 se requieren tramos múltiplos de 3, por lo que faltaría un tramo para volver a usar la fórmula.

    La solución por filas se implementa usando una combinación de Simpson 3/8 para los puntos entre remover[i, 0:(2+1)] y Simpson 1/3 para los puntos entre remover[i, 3:(4+1)].

    Luego se completa el integral del otro eje con el resultado anterior, aplicando el mismo método.

    Se obtuvieron los siguientes resultados:

    Integral en eje x: 
    [435.1  346.5  287.44 255.58 590.54 440.52]
    Volumen:  199161.11111111115
    

    que se obtiene usando las instrucciones a continuación de las anteriores:

    # literal b) ---------------------------
    def simpson13(fi,h):
        s13 = (h/3)*(fi[0] + 4*fi[1] + fi[2])
        return(s13)
    def simpson38(fi,h):
        s38 = (3*h/8)*(fi[0] + 3*fi[1] + 3*fi[2]+ fi[3])
        return(s38)
    
    Integralx = np.zeros(n,dtype = float)
    
    for i in range(0,n,1):
        hx = xi[1]-xi[0]
        fi = remover[i, 0:(0+4)]
        s38 = simpson38(fi,hx)
        fi = remover[i, 3:(3+3)]
        s13 = simpson13(fi,hx)
        Integralx[i] = s38 + s13
    
    hy = yi[1] - yi[0]
    fj = Integralx[0:(0+4)]
    s38 = simpson38(fj,hy)
    fj = Integralx[3:(3+3)]
    s13 = simpson13(fj,hy)
    volumen = s38 + s13
    
    # Salida
    np.set_printoptions(precision=2)
    print('Integral en eje x: ')
    print(Integralx)
    print('Volumen: ', volumen)
    

    Para el examen escrito, se requieren realizar al menos 3 iteraciones/ filas para encontrar el integral.

  • s2Eva_IT2018_T1 EDO Paracaidista wingsuit

    Ejercicio: 2Eva_IT2018_T1 EDO Paracaidista wingsuit

    Plantear con: [ Runge-Kutta para f''(x) ] [ Runge-Kutta para f'(x) ]

    ..


    a. Planteamiento con Runge-Kutta 2do Orden para Segunda derivada

    La expresión:

    \frac{dv}{dt} = g - \frac{cd}{m} v^2

    se puede plantear sustituir la variable con v = -\frac{dy}{dt} al considerar el sentido contrario entre la velocidad al caer y la referencia de altura hacia arriba. Ver figura.

    \frac{d^2y}{dt^2} = g - \frac{cd}{m} \Big( \frac{dy}{dt} \Big) ^2

    Que es una EDO de 2do orden o como 2da derivada.diagrama cuerpo libre

    La solución se propone resolver de forma simultanea para t,y,v con Runge Kutta para segunda derivada donde:

    f(t,y,v) = -v g(t,y,v) = g - \frac{cd}{m} v^2

    Al sustituir los valores de las constantes en la ecuación como gravedad, masa e índice de arrastre se tiene:

    f(t,y,v) = -v g(t,y,v) = 9.8 - \frac{0.225}{90} v^2

    con las condiciones iniciales del ejercicio  t0 = 0 , y0 = 1000, v0 = 0
    la velocidad se inicia con cero, si el paracaidista se deja caer desde el borde el risco, como en el video adjunto al enunciado.

    Para las iteraciones, recuerde que
    t se corrige con t+h (en el algoritmo era la posición para x)
    y se corrige con y+K1y
    v se corrige con v+K1v (en el algoritmo era la posición para z)

    itera = 0

    K1y = h f(t,y,v) = 2(-(0)) = 0 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (0)^2) = 19.6

    ..
    K2y = h f(t+h, y+K1y, v + K1v) = 2(-(0 + 19.6)) = -39.2

    K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (0+19.6)^2) =17.6792

    ..
    y_1 = y_0 + \frac{K1y+K2y}{2} = 1000 + \frac{0-39.2}{2}= 980.4

    v_1 = v_0 + \frac{K1v+K2v}{2} = 0 + \frac{19.6-17.67}{2} = 18.63 t_1 =t_0 + h = 0+2 = 2

     

    ti yi vi K1y K1v K2y K2v
    0 1000 0 0 0 0 0
    2 980.4 18.63 0 19.6 -39.2 17.6792

    itera = 1

    K1y = 2(-(18.63)) = -37.2792 K1v = 2(9.8 - \frac{0.225}{90} (18.63)^2) = 17.8628

    ..
    K2y =2(-(18.6396+17.8628)) =-73.00

    K2v = 2(9.8 - \frac{0.225}{90} (18.6396+17.8628)^2) =12.9378

    ..
    y_2 =980.4 + \frac{ -37.2792+(-73.00)}{2}= 925.25

    v_2 = 18.63 + \frac{17.8628+12.9378}{2} = 34.0399 t_2 =t_1 + h = 2+2 = 4
    ti yi vi K1y K1v K2y K2v
    0 1000 0 0 0 0 0
    2 980.4 18.63 0 19.6 -39.2 17.6792
    4 925.25 34.0399 -37.2792 17.8628 -73.00 12.9378

    itera = 2

    K1y = h f(t,y,v) = 2(-(34.0399)) = -68.0798 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (34.0399)^2) = 13.8064

    ..
    K2y = h f(t+h, y+K1y, v + K1v) = 2(-(34.0399+13.8064)) =-95.6927

    K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (34.0399+13.8064)^2) =8.1536

    ..
    y_2 = 925.25 + \frac{ -68.0798+(-95.6927)}{2}= 843.3716

    v_2 = 34.0399 + \frac{13.8064+8.1536}{2} = 45.0199 t_2 =t_1 + h = 4+2 = 6
    ti yi vi K1y K1v K2y K2v
    0 1000 0 0 0 0 0
    2 980.4 18.63 0 19.6 -39.2 17.6792
    4 925.25 34.0399 -37.2792 17.8628 -73.00 12.9378
    6 843.3716 45.0199 -68.0798 13.8064 -95.6927 8.1536

    Algoritmo con Python

    Resultados
    La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
    Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de  62 m/s, que son algo mas de 223 Km/h, el objeto se estrella...!

    Resultados con el algoritmo:

    Runge-Kutta Segunda derivada
    i  [ xi,  yi,  zi ]
       [ K1y,  K1z,  K2y,  K2z ]
    0 [   0. 1000.    0.]
       [0. 0. 0. 0.]
    1 [  2.     980.4     18.6396]
      [  0.      19.6    -39.2     17.6792]
    2 [  4.       925.257973  34.039945]
      [-37.2792    17.862827 -73.004853  12.937864]
    3 [  6.       843.371672  45.019966]
      [-68.079891  13.806411 -95.692712   8.153631]
    4 [  8.       743.865726  52.131168]
      [ -90.039933    9.466013 -108.971959    4.75639 ]
    5 [ 10.       633.591684  56.485537]
      [-104.262336    6.011707 -116.285749    2.697031]
    6 [ 12.       516.97369   59.069216]
      [-112.971073    3.646921 -120.264915    1.520438]
    7 [ 14.       396.681119  60.575537]
      [-118.138432    2.154139 -122.446709    0.858504]
    8 [ 16.       274.277023  61.445121]
      [-121.151075    1.253021 -123.657117    0.486147]
    9 [ 18.       150.664295  61.944336]
      [-122.890243    0.722485 -124.335213    0.275943]
    10 [20.       26.361127 62.230024]
       [-123.888671    0.414496 -124.717664    0.15688 ]
    11 [ 22.       -98.336041  62.393224]
       [-1.244600e+02  2.371205e-01 -1.249343e+02  8.927924e-02]
    12 [  24.       -223.257917   62.486357]
       [-1.247864e+02  1.354280e-01 -1.250573e+02  5.083841e-02]
    13 [  26.       -348.307907   62.539475]
       [-1.249727e+02  7.727584e-02 -1.251273e+02  2.895913e-02]
    14 [  28.       -473.430927   62.56976 ]
       [-1.250789e+02  4.407055e-02 -1.251671e+02  1.649935e-02]
    15 [  30.       -598.595572   62.587023]
       [-1.251395e+02  2.512592e-02 -1.251898e+02  9.401535e-03]
    16 [  32.       -723.783942   62.596863]
       [-1.251740e+02  1.432256e-02 -1.252027e+02  5.357469e-03]
    >>> 
    

    Instrucciones en Python

    # 2Eva_IT2018_T1 Paracaidista wingsuit
    import numpy as np
    import matplotlib.pyplot as plt
    
    def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                       vertabla=False, precision = 6):
        ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
            x0,y0 son valores iniciales, h es tamano de paso,
            muestras es la cantidad de puntos a calcular.
        '''
        tamano = muestras + 1
        tabla = np.zeros(shape=(tamano,3+4),dtype=float)
    
        # incluye el punto [x0,y0,z0]
        tabla[0] = [x0,y0,z0,0,0,0,0]
        xi = x0
        yi = y0
        zi = z0
        i=0
        if vertabla==True:
            print('Runge-Kutta Segunda derivada')
            print('i ','[ xi,  yi,  zi',']')
            print('   [ K1y,  K1z,  K2y,  K2z ]')
            np.set_printoptions(precision)
            print(i,tabla[i,0:3])
            print('  ',tabla[i,3:])
        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
            
            tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
            if vertabla==True:
                txt = ' '
                if i>=10:
                    txt='  '
                print(str(i)+'',tabla[i,0:3])
                print(txt,tabla[i,3:])
        return(tabla)
    
    
    # PROGRAMA PRUEBA
    # INGRESO
    f = lambda t,y,v: -v # el signo, revisar diagrama cuerpo libre
    g = lambda t,y,v: 9.8 - (0.225/90)*(v**2)
    t0 = 0
    y0 = 1000
    v0 = 0
    h  = 2
    muestras = 15+1
    
    # PROCEDIMIENTO
    tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras, vertabla=True)
    ti = tabla[:,0]
    yi = tabla[:,1]
    vi = tabla[:,2]
    
    # SALIDA
    # print('tabla de resultados')
    # print(tabla)
    
    # GRAFICA
    plt.subplot(211)
    plt.plot(ti,vi,label='velocidad v(t)', color='green')
    plt.plot(ti,vi,'o')
    plt.ylabel('velocidad (m/s)')
    plt.title('paracaidista Wingsuit con Runge-Kutta')
    plt.legend()
    plt.grid()
    
    plt.subplot(212)
    plt.plot(ti,yi,label='Altura y(t)',)
    plt.plot(ti,yi,'o',)
    plt.axhline(0, color='red')
    plt.xlabel('tiempo (s)')
    plt.ylabel('Altura (m)')
    plt.legend()
    plt.grid()
    
    plt.show()
    

    Plantear con: [ Runge-Kutta para f''(x) ] [ Runge-Kutta para f'(x) ]

    ..


    b. Usando Runge-Kutta 2do Orden para Primera Derivada o velocidad,

    El problema para un tiempo de observación t>0, se puede dividir en dos partes: velocidad y altura.

    1. Determinar velocidad: Se aplica Runge-Kutta a la expresión con primera derivada o velocidad. Use  los valores iniciales dados, descarte calcular las alturas.
    2. Determinar las altura:  Con los valores de velocidades y la altura inicial de 1 km = 1000 m puede integrar con trapecio para obtener la tabla de alturas. Se integra tramo a tramo.

    Observe las unidades de medida y que la  velocidad es contraria  al eje de altura dy/dt = -v.

     

    La expresión:

    \frac{dv}{dt} = g - \frac{cd}{m} v^2 f(t,v) = g - \frac{cd}{m} v^2

    itera = 0

    K1v = h f(t,v) = 2(9.8)- \frac{0.225}{90}(0)^2) = 19.6 K2v = h f(t+h , v + K1v) = 2( 9.8 - \frac{0.225}{90}(0+19.6)^2) = 17.6792 v_1 = v_0 + \frac{K1y+K2y}{2} = 0 + \frac{19.6+17.6792}{2}= 18.6396 t_1 =t_0 + h = 0+2 = 2

    itera = 1

    K1v = 2(9.8 - \frac{0.225}{90} (18.6396)^2) = 17.8628 K2v = 2( 9.8 - \frac{0.225}{90} (18.6396+17.8628)^2) = 12.9379 v_1 = 18.6396 + \frac{17.8628+12.9379}{2}= 34.0399 t_1 =2+2 = 4

    itera = 2

    K1v = 2(9.8 - \frac{0.225}{90} (34.0399)^2) = 13.8064 K2v = 2( 9.8 - \frac{0.225}{90} (34.0399+13.8064)^2) = 8.1536 v_1 = 34.0399 + \frac{13.8064+8.1536}{2}= 45.02 t_1 = 4+2 = 6

    las siguientes iteraciones se completan con el algoritmo.

    Resultados

    La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
    Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de  62 m/s, que son algo mas de 223 Km/h, el objeto se estrella...!

    paracaidista wingsuit 02

    velocidad con Runge-Kutta primera derivada
     [tiempo, velocidad, K1,K2]
    [[ 0.      0.      0.      0.    ]
     [ 2.     18.6396 19.6    17.6792]
     [ 4.     34.0399 17.8628 12.9379]
     [ 6.     45.02   13.8064  8.1536]
     [ 8.     52.1312  9.466   4.7564]
     [10.     56.4855  6.0117  2.697 ]
     [12.     59.0692  3.6469  1.5204]
     [14.     60.5755  2.1541  0.8585]
     [16.     61.4451  1.253   0.4861]
     [18.     61.9443  0.7225  0.2759]
     [20.     62.23    0.4145  0.1569]
     [22.     62.3932  0.2371  0.0893]]
    velocidad con RK2 y altura con trapecio
     [tiempo, velocidad, altura]
    [[   0.      0.   1000.  ]
     [   2.     18.64  981.36]
     [   4.     34.04  928.68]
     [   6.     45.02  849.62]
     [   8.     52.13  752.47]
     [  10.     56.49  643.85]
     [  12.     59.07  528.3 ]
     [  14.     60.58  408.65]
     [  16.     61.45  286.63]
     [  18.     61.94  163.24]
     [  20.     62.23   39.07]
     [  22.     62.39  -85.55]]
    >>> 
    

    Los cálculos se realizaron usando las instrucciones en Python:

    # 2da Evaluación I Término 2018
    # Tema 1. Paracaidista wingsuit
    import numpy as np
    
    def rungekutta2(d1y,x0,y0,h,muestras):
        # Runge Kutta de 2do orden
        tamano = muestras + 1
        tabla = np.zeros(shape=(tamano,2+2),dtype=float)
        
        # incluye el punto [x0,y0]
        tabla[0] = [x0,y0,0,0]
        xi = x0
        yi = y0
        for i in range(1,tamano,1):
            K1 = h * d1y(xi,yi)
            K2 = h * d1y(xi+h, yi + K1)
    
            yi = yi + (1/2)*(K1+K2)
            xi = xi + h
            
            tabla[i] = [xi,yi,K1,K2]
        return(tabla)
    
    def integratrapecio_fi_tabla(xi,fi,y0):
        tamano = len(xi)
        yi = np.zeros(tamano,dtype = float)
        yi[0] = y0
        for i in range(1,tamano,1):
            h = xi[i]-xi[i-1]
            trapecio = h*(fi[i]+fi[i-1])/2
            yi[i]= yi[i-1] + trapecio
        return(yi)
    
    # PROGRAMA -------------------------
    
    # INGRESO
    g = 9.8
    cd = 0.225
    m = 90
    d1v = lambda t,v: g - (cd/m)*(v**2)
    
    t0 = 0
    v0 = 0
    h = 2
    y0 = 1000
    muestras = 11
    
    # PROCEDIMIENTO
    velocidad = rungekutta2(d1v,t0,v0,h,muestras)
    ti = velocidad[:,0]
    vi = velocidad[:,1]
    
    # Altura, velocidad es contraria altura,
    # integrar en cada tramo por trapecios o Cuadratura de Gauss
    altura = integratrapecio_fi_tabla(ti,-vi,y0)
    
    # Tabla de resultados de tiempo, velocidad, altura
    altura = np.transpose([altura])
    tabla = np.concatenate((velocidad[:,0:2],altura), axis = 1)
    
    # SALIDA
    np.set_printoptions(precision=4)
    print('velocidad con Runge-Kutta primera derivada')
    print(' [tiempo, velocidad, K1,K2]')
    print(velocidad)
    np.set_printoptions(precision=2)
    print('velocidad con RK2 y altura con trapecio')
    print(' [tiempo, velocidad, altura]')
    print(tabla)
    

    Plantear con: [ Runge-Kutta para f''(x) ] [ Runge-Kutta para f'(x) ]

  • s1Eva_IT2018_T1 Tanque esférico canchas deportivas

    Ejercicio: 1Eva_IT2018_T1 Tanque esférico canchas deportivas

    a) Para seleccionar el rango para h=[a,b], se observa que el tanque puede estar vacío, a=0 o lleno al máximo, b=2R = 2(3)=6, con lo que se obtiene:

    h =[0.0, 6.0]

    conociendo la proporción con el valor máximo, se tiene un valor inicial para h0 para las iteraciones.

    Vmax = \frac{\pi}{3} (2R)^2 (3R-2R) Vmax = \frac{4\pi }{3}R^{3}= 113.09 h_0 = (6)*30/113.09 = 1.59

    b) Usar Newton-Raphson con tolerancia 1e-6

    f(h) = \frac{\pi }{3}h^2 (3(3)-h)-30 f(h) = \frac{\pi }{3}(9h^2 -h^3-28.647) f'(h) = \frac{\pi }{3}(18h-3h^2)

    el punto siguiente de iteración es:

    h_{i+1} = h_{i} -\frac{f(h)}{f'(h)} = h_{i}-\frac{ \frac{\pi }{3}(9h^2 -h^3-28.647)}{ \frac{\pi }{3}(18h-3h^2)} h_{i+1} = h_{i} -\frac{(9h^2 -h^3-28.647)}{(18h-3h^2)}

    con lo que se realiza la tabla de iteraciones:

    hi hi+1 error orden
    1.590 2.061 0.47 10-1
    2.061 2.027 -0.034 10-2
    2.027 2.02686 -0.00014 10-4
    2.02686 2.0268689 -2.32E-09 10-9

    En valor práctico 2.028 m usando flexómetro, a menos que use medidor laser con precisión 10-6 usará más dígitos con un tanque de 6 metros de altura ganará una precisión de una gota de agua para usar en duchas o regar el césped .

    c) El orden de convergencia del error observando las tres primeras iteraciones es cuadrático

    Tarea: Realizar los cálculos con Python, luego aplique otro método. Añada sus observaciones y conclusiones.

  • s1Eva_IT2018_T3 Temperatura en nodos de placa

    Ejercicio: 1Eva_IT2018_T3 Temperatura en nodos de placa

    a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior: Placa Temp 03

    a=\frac{50+c+100+b}{4} b=\frac{a+30+50+d}{4} c=\frac{a+60+100+d}{4} d=\frac{b+60+c+30}{4}

    que reordenando se convierte en:

    4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90

    simplificando:

    4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90

    que a forma matricial se convierte en:

    A = [[ 4, -1, -1, 0.0],
         [ 1, -4,  0, 1.0],
         [ 1,  0, -4, 1.0],
         [ 0,  1,  1,-4.0]]
    B = [[ 150.0],
         [ -80.0],
         [-160.0],
         [ -90.0]]
    

    Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas.  Se aumentó el punto decimal a los valores de la matriz A y el vector B  para que sean considerados como números reales.

    El número de condición es: np.linalg.cond(A) = 3.0

    que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es "estable" y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.

    b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial

     
    X(0) = [[60.0],
            [40], 
            [70],
            [50]] 
    

    reemplazando los valores iniciales en cada ecuación sin cambios.

    iteración 1
    a=\frac{50+70+100+40}{4} = 65

    b=\frac{60+30+50+50}{4} = 47.5 c=\frac{60+60+100+50}{4} = 67.5 d=\frac{40+60+70+30}{4} = 50
    X(1) = [[65],
            [47.5],
            [67.5],
            [50]]
    
    vector de error = 
         [|65-60|,
          |47.5-40|,
          |67.5-70|,
          |50-50|]
      = [|5|,
         |7.5|,
         |-2.5|,
         |0|]
    errormax = 7.5
    

    iteración 2
    a=\frac{50+67.5+100+47.5}{4} = 66.25

    b=\frac{65+30+50+50}{4} = 48.75 c=\frac{65+60+100+50}{4} = 68.75 d=\frac{47.5+60+67.7+30}{4} = 51.3
    X(2) = [[66.25],
            [48.75],
            [68.75],
            [51.3]]
    
    vector de error = 
           [|66.25-65|,
            |48.75-47.5|,
            |68.75-67.5|,
               |51.3-50|] 
      = [|1.25|,
         |1.25|,
         |1.25|,
         |1.3|]
    errormax = 1.3
    

    iteración 3
    a=\frac{50+68.75+100+48.75}{4} = 66.875

    b=\frac{66.25+30+50+51.3}{4} = 49.38 c=\frac{66.25+60+100+51.3}{4} = 69.3875 d=\frac{48.75+60+68.75+30}{4} = 51.875
    X(2) = [[66.875],
            [49.38],
            [69.3875],
            [51.875]]
    
    vector de error = 
          [|66.875-66.25|,
           |49.38-48.75|,
           |69.3875-68.75|,
           |51.875-51.3|]
        = [|0.655|,
           |0,63|,
           |0.6375|,
            |0.575|]
    errormax = 0.655
    con error relativo de:
    100*0.655/66.875 = 0.97%

    siguiendo las iteraciones se debería llegar a:

    >>> np.linalg.solve(A,B)
    array([[ 67.5],
           [ 50. ],
           [ 70. ],
           [ 52.5]])
    
  • s1Eva_IT2018_T4 El gol imposible

    Ejercicio: 1Eva_IT2018_T4 El gol imposible

    Tabla de datos:

    ti = [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20]
    xi = [0.00, 0.50, 1.00, 1.50, 1.80, 2.00, 1.90, 1.10, 0.30]
    yi = [0.00, 4.44, 8.88,13.31,17.75,22.19,26.63,31.06,35.50]
    zi = [0.00, 0.81, 1.40, 1.77, 1.91, 1.84, 1.55, 1.03, 0.30]
    

    Observe que, un gol simplificado con física básica es un tiro parabólico cuya trayectoria se compone de movimientos en los ejes, Y y Z. Sin embargo, lo "imposible" del gol mostrado es añadir el movimiento en X. (Referencia de la imagen en el enunciado)

    El movimiento de "profundidad" o dirección hacia el arco y(t) es semejante a un polinomio de primer grado, y el movimiento de "altura" z(t) es un polinomio de segundo grado. El movimiento "imposible" en el eje X, podría ser descrito por un polinomio de segundo o mayor grado.

    a) Encontrar t para altura máxima, que se encuentra al igualar la derivada dz/dt a cero. Para interpolar el polinomio z(t), de segundo grado, se puede usar tres puntos de los sugeridos: 0, 0.3 y 0.6, que además son equidistantes en t (facilita usar cualquier método de interpolación).

    Por ejemplo, con diferencias finitas avanzadas:

    t z(t) d1 d2 d3
    0.00 0.00 1.40 -0.89
    0.30 1.40 0.51
    0.60 1.91
    z(t) = 0 + 1.40\frac{(t-0)}{1!(0.3)} - 0.89 \frac{(t-0)(t-0.3)}{2!(0.3)^2} = 0 + 4.66 t - 4.94(t^2-0.3t) = 4.66 t + 1.48 t - 4.94 t^2 z(t) = 6.42 t - 4.94 t^2

    para encontrar el valor máximo se encuentra \frac{dz}{dt} = 0

    \frac{dz}{dt} = 6.42 - 2(4.94) t 6.42 - 2(4.94) t = 0 t = \frac{6.42}{2(4.94)} t = 0.649

    Observe que el resultado tiene sentido, pues según la tabla, el máximo debería estar entre 0.60 y 0.75

    b) El cruce por la "barrera", corresponde al desplazamiento del balón en el eje Y a 9 metros: y(t) = 9.
    El polinomio modelo de primer grado usa dos puntos para la interpolación, de los sugeridos se escogen dos, por ejemplo: 0.0 y 0.3.

    Usando diferencias finitas avanzadas :

    d1 = (8.88-0) = 8.88 y(t) = 0 + 8.88\frac{(t-0)}{1!(0.3)} y(t) = 29.6 t

    usando y(t) = 9

    29.6 t = 9
    t = 0.30
    z(0.30) = 1.40
    (de la tabla de datos)

    cuya respuesta con solo dos dígitos decimales es coherente, al estar cercano el valor a una distancia y=8.88 o aproximado a 9.
    La respuesta puede ser mejorada usando más digitos significativos en las operaciones.

    c)  La desviación máxima en eje X.
    Determine un polinomio para la trayectoria en el eje X y obtenga el máximo igualando la derivada del polinomio x(t) a cero.

    Por simplicidad, se usa un polinomio de segundo grado, alrededor de los valores máximos en el eje X

    t x(t) d1 d2 d3
    0.60 1.80 0.20 -0.30
    0.75 2.00 -0.10
    0.90 1.90
    x(t) = 1.80 + 0.20 \frac{(t-0.60)}{1!(0.15)} -0.30 \frac{(t-0.60)(t-0.75)}{2!(0.15)^2} x(t) = 1.80 + 1.33 (t-0.60) - 6.67(t-0.60)(t-0.75)

    como se busca el máximo, se usa la derivada \frac{dx}{dt} = 0

    \frac{dx}{dt} = 1.33 - 6.67(2t +(-0.60-0.75)) 1.33 - 13.34t + 9.00 = 0 -13.34t + 10.33 = 0

    t = 0.77

    x(0.77) = 1.80 + 1.33(0.77-0.60) - 6.67(0.77-0.60)(0.77-0.75) x(0.77) = 2.003

    lo que es coherente con la tabla para el eje x, pues el máximo registrado es 2, y el próximo valor es menor, la curva será decreciente.


    Desarrollo extra para observar y verificar resultados:

    Usando los puntos y un graficador 3D se puede observar la trayectoria:

    Gol Imposible 02

    Tarea: Realice el ejercicio usando los algoritmos en Python, muestre los polinomios obtenidos y compare.

    Nota: La gráfica 3D presentada usa interpolación de Lagrange con todos los puntos. Realice las observaciones y recomendaciones necesarias y presente su propuesta como tarea.

  • s1Eva_IIT2017_T1 Aproximar a polinomio usando puntos

    Ejercicio: 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos

    Se dispone de tres puntos para la gráfica.

    x  f(x)
     0  1
     0.2  1.6
     0.4  2.0

    Si el polinomio de Taylor fuera de grado 0, sería una constante, que si se evalúa en x0 = 0 para eliminar los otros términos, se encuentra que sería igual a 1

    Como se pide el polinomio de grado 2, se tiene la forma:

    p(x) = a + bx + c x ^2 p(x) = 1 + bx + c x^2

    Se disponen de dos puntos adicionales que se pueden usar para determinar b y c:

    p(0.2) = 1 + 0.2 b + (0.2)^2 c = 1.6 p(0.4) = 1 + 0.4 b + (0.4)^2 c = 2.0

    simplificando:

    0.2 b + (0.2)^2 c = 1.6 - 1 = 0.6 0.4 b + (0.4)^2 c = 2.0 - 1 = 1

    multiplicando la primera ecuación por 2 y restando la segunda ecuación:

    0 - 0.08 c = 1.2-1 = 0.2 c = - 0.2/0.08 = -2.5

    sustituyendo el valor de c obtenido en la primera ecuación

    0.2 b + 0.04(-2.5) = 0.6 0.2 b = 0.6 - 0.04(-2.5) = 0.6 + 0.1 = 0.7 b = 0.7/0.2 = 3.5

    con lo que el polinomio queda:
    p(x) = 1 + 3.5 x - 2.5 x^2

    validando con python:
    tomando los puntos de prueba:

    xi = [ 0, 0.2, 0.4]
    fi = [ 1, 1.6, 2 ]
    >>>

    se obtiene la gráfica:

    se adjunta las instrucciones usadas para validar que el polinomio pase por los puntos requeridos.

    # 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    xi = [ 0, 0.2, 0.4]
    fi = [ 1, 1.6, 2 ]
    
    px = lambda x: 1 + 3.5*x - 2.5*(x**2)
    a = -0.5
    b = 1
    muestras = 21
    
    # PROCEDIMIENTO
    xj = np.linspace(a,b,muestras)
    pxj = px(xj)
    
    # SALIDA
    print(xj)
    print(pxj)
    
    # Gráfica
    plt.plot(xj,pxj,label='p(x)')
    plt.plot(xi,fi,'o', label='datos')
    
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.grid()
    plt.legend()
    plt.show()
    
    

    Nota: Se puede intentar realizar los polinomios aumentando el grado, sin embargo cada término agrega un componente adicional a los términos anteriores por la forma (x - x0)k

    literal b

    se requiere el integral aproximado de f(x) en el intervalo limitado por los 3 puntos de la tabla:

    \int_{0}^{0.4}f(x) dx

    Esta aproximación con un polinomio es el concepto de integración numérica con la regla de Simpson de 1/3, tema desarrollado en la unidad 5

    I_{S13} = \frac{0.2}{3} \Big(1+4(1.6)+ 2 \Big) = 0.62666
  • s2Eva_IIT2017_T4 Problema con valor de frontera

    Ejercicio: 2Eva_IIT2017_T4 Problema con valor de frontera

    \frac{d^2T}{dx^2} + \frac{1}{x}\frac{dT}{dx} +S =0 0 \leq x \leq 1

    Las diferencias finitas a usar son:

    \frac{dT}{dx} =\frac{T_{i+1} - T_{i-1}}{2h}+O(h^2) \frac{d^2T}{dx^2}=\frac{T_{i+1} - 2T_i + T_{i-1}}{h^2}+O(h^2)

    que al reemplazar el la ecuación:

    \frac{T_{i+1} - 2T_i + T_{i-1}}{h^2} + \frac{1}{x_i}\frac{T_{i+1} -T_{i-1}}{2h}+S=0 2x_i (T_{i+1} - 4h x_i T_i + T_{i-1}) + h (T_{i+1} - T_{i-1}) = -2h^2 S x_i T_{i+1}(2x_i + h) - 4x_i T_i + T_{i-1}(2x_i - h) = -2h^2 S x_i T_{i-1}(2x_i - h) - 4x_i T_i + T_{i+1}(2x_i + h)= -2h^2 S x_i

    con lo que se puede crear un sistema de ecuaciones para cada valor xi con T0=2 y T4 =1 que son parte del enunciado del problema.

    Siendo h = 0.25 = 1/4,  y se indica al final que S=1, se crea un sistema de ecuaciones a resolver,

    x = [0, 1/4, 1/4, 3/4, 1]
    T_{i-1}\Big[2x_i - \frac{1}{4} \Big] - 4x_i T_i + T_{i+1}\Big[2x_i + \frac{1}{4} \Big] = -2\Big(\frac{1}{4}\Big)^2 (1) x_i T_{i-1}\Big[2x_i -\frac{1}{4}\Big] - 4x_i T_i + T_{i+1}\Big[2x_i + \frac{1}{4}\Big] = -\frac{1}{8} x_i

    se sustituye con los valores conocidos para cada i:

    i=1: 
    T0[2(1/4) - 1/4] - 4(1/4)T1 + T2[2(1/4) + 1/4] = -(1/8)(1/4)
    
         - T1 + (3/4)T2 = -1/32 - (1/4)(2)
         - T1 + (3/4)T2 = -17/32
    
    i=2: 
    T1[2(1/2) - 1/4] - 4(1/2)T2 + T3[2(1/2) + 1/4] = -(1/8)(1/2)
    
         (3/4)T1 - 2T2 + (5/4)T3 = -1/16
    
    i=3: 
    T2[2(3/4) - 1/4] - 4(3/4)T3 + T4[2(3/4) + 1/4] = -(1/8)(3/4)
    
         (5/4)T2 - 3T3 = -3/32 - (7/4)(1)
         (5/4)T2 - 3T3 = -59/32
    

    se ponen las ecuaciones en matrices para resolver con algun metodo numérico:

    A = [[ -1, 3/4,   0],
         [3/4,  -2, 5/4],
         [  0, 5/4,  -3]]
    B = [-17/32, -1/16, -59/32]
    np.linalg.solve(A,B)
    array([ 1.54,  1.34,  1.17])