Etiqueta: análisis numérico

  • s2Eva_IIT2018_T2 EDO d2x/dt2 Kunge Kutta 2do Orden

    Ejercicio: 2Eva_IIT2018_T2 EDO d2x/dt2 Kunge Kutta 2do Orden

    \frac{\delta ^2 x}{\delta t^2} + 5t\frac{\delta x}{\delta t} +(t+7)\sin (\pi t) = 0 x'' + 5tx' +(t+7)\sin (\pi t) = 0 x'' = -5tx' -(t+7)\sin (\pi t) = 0

    si se usa z=x'

    z' = -5tz -(t+7)\sin (\pi t) = 0

    se convierte en:

    f(t,x,z) = x' = z g(t,x,z) = x'' = z' = -5tz -(t+7)sin (\pi t) = 0

    Tarea: Desarrollar 3 iteraciones en Papel.

    Donde se aplica el algoritmo de Runge Kutta
    https://blog.espol.edu.ec/analisisnumerico/8-2-2-runge-kutta-d2y-dx2/

       t,              x,              z
    [[ 0.          6.          1.5       ]
     [ 0.2         6.3         0.92679462]
     [ 0.4         6.38218195 -0.27187703]
     [ 0.6         6.19792527 -1.17287944]
     [ 0.8         5.88916155 -1.23638799]
     [ 1.          5.6491005  -0.61819399]
     [ 1.2         5.5872811   0.17288691]
     [ 1.4         5.69750883  0.69945284]
     [ 1.6         5.8992535   0.77223688]
     [ 1.8         6.09372469  0.43437943]
     [ 2.          6.20586248 -0.12630953]]
    

    Instrucciones en Python

    # 2Eva_IIT2018_T2 Kunge Kutta 2do Orden x''
    import numpy as np
    
    def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
        tamano = muestras + 1
        estimado = np.zeros(shape=(tamano,3),dtype=float)
        # incluye el punto [x0,y0]
        estimado[0] = [x0,y0,z0]
        xi = x0
        yi = y0
        zi = z0
        for i in range(1,tamano,1):
            K1y = h * 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]
        return(estimado)
    
    # PROGRAMA
    # INGRESO
    f = lambda t,x,z: z
    g = lambda t,x,z: -5*t*z-(t+7)*np.sin(np.pi*t)
    t0 = 0
    x0 = 6
    z0 = 1.5
    h = 0.2
    muestras = 10
    
    # PROCEDIMIENTO
    tabla = rungekutta2_fg(f,g,t0,x0,z0,h,muestras)
    
    # SALIDA
    print(tabla)
    # GRAFICA
    import matplotlib.pyplot as plt
    plt.plot(tabla[:,0],tabla[:,1])
    plt.xlabel('t')
    plt.ylabel('x(t)')
    plt.show()
    
  • s2Eva_IIT2018_T1 Masa entra o sale de un reactor

    Ejercicio: 2Eva_IIT2018_T1 Masa entra o sale de un reactor

    a) Se pueden combinar los métodos para realizar la integral. Se usa el método de Simpson 1/3 para los primeros dos tramos y Simpson 3/8 para los 3 tramos siguientes.  Siendo f(x) equivalente a Q(t)C(t). El tamaño de paso h es constante para todo el ejercicio con valor 5.

    a.1 Simpson 1/3, tramos 2, puntos 3:

    I_1 \cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)] I_1 \cong \frac{5}{3}[(10)(4)+4(18)(6) + (27)(7)] I_1 \cong 1101,66

    a.2 Simpson de 3/8, tramos 3, puntos 4:

    I_2 \cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] I_2 \cong \frac{3(5)}{8}[(27)(7)+3(35)(6) +3(40)(5)+(30)(5)] I_2 \cong 2941,88 I_1 + I_2 \cong 4043,54

    b) El error se calcula por tramo y se acumula.

    b.1 se puede estimar como la diferencia entre la parábola del primer tramo y simpson 1/3
    b.2 siguiendo el ejemplo anterior, como la diferencia entre la interpolación de los tramos restantes y simpson 3/8.

  • s2Eva_IT2010_T2 EDO Movimiento angular

    Ejercicio: 2Eva_IT2010_T2 EDO Movimiento angular

    Para resolver, se usa Runge-Kutta_fg de segundo orden como ejemplo

    y'' + 10 \sin (y) =0

    se hace

    y' = z = f(t,y,z)

    y se estandariza:

    y'' =z'= -10 \sin (y) = g(t,y,z)

    teniendo como punto de partida t0=0, y0=0 y z0=0.1

    y(0)=0, y'(0)=0.1

    Se desarrolla el algoritmo para obtener los valores:

     [ t, 		 y, 	 dyi/dti=z]
    [[ 0.          0.          0.1       ]
     [ 0.2         0.02        0.08000133]
     [ 0.4         0.03200053  0.02401018]
     [ 0.6         0.03040355 -0.04477916]
     [ 0.8         0.01536795 -0.09662411]
     [ 1.         -0.00703034 -0.10803459]]

    que permiten generar la gráfica de respuesta:

    Movimiento Angular 01


    Algoritmo en Python

    # 2Eva_IT2010_T2 Movimiento angular
    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,3),dtype=float)
        # incluye el punto [x0,y0,z0]
        estimado[0] = [x0,y0,z0]
        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]
        return(estimado)
    
    # INGRESO theta = y
    ft = lambda t,y,z: z
    gt = lambda t,y,z: -10*np.sin(y)
    
    t0 = 0
    y0 = 0
    z0 = 0.1
    h=0.2
    muestras = 5
    
    # PROCEDIMIENTO
    tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)
    
    # SALIDA
    print(' [ t, \t\t y, \t dyi/dti=z]')
    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.xlabel('ti')
    plt.title('yi')
    plt.subplot(122)
    plt.plot(ti,zi, color='green')
    plt.xlabel('ti')
    plt.title('dyi/dti')
    plt.show()
    
  • s1Eva_IIT2018_T4 Tasa de interés en hipoteca

    Ejercicio: 1Eva_IIT2018_T4 Tasa de interés en hipoteca

    literal a

    Siguiendo el desarrollo analítico tradicional, para adecuar la ecuación para los algoritmo de búsquda de raíces de ecuaciones,  se reemplazan los valores en la fórmula.

    P = A\Big(\frac{1-(1+i)^{-n}}{i} \Big) 70000 = 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)

    Como ambos lados de la ecuación deben ser iguales, si se restan ambos se obtiene una ecuación que tiene como resultado cero, que es la forma ideal para usar en el algoritmo que representa f(x) o en este caso f(i)

    70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) = 0

    Para evitar inconvenientes con la división para cero en caso que i tome el valor de cero, dado se multiplica toda la ecuación por i:

    i \Big[70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) \Big]= i (0) 70000 i - 1200 (1-(1+i)^{-300}) = 0

    La ecuación es la utilizada en el algoritmo de búsqueda de raíces pueden ser:

    fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) fx(i) = 70000i - 1200(1-(1+i)^{-300})

    literal b

    El intervalo de existencia correspondería a la tasa de interés mínimo y el interés máximo.

    [izquierda, derecha] = [a,b]

    Para el intervalo se deben tomar en cuenta algunas consideraciones descritas a continuación:

    izquierda:

    En el extremo izquierdo, las tasas no son negativas, lo que se interpreta en que un banco paga por que le presten dinero.

    Tampoco tiene mucho sentido el valor cero, que son prestamos sin intereses. A menos que sean sus padres quienes le dan el dinero.

    Un valor inicial para el interés puede ser por ejemplo 1% ó 0.01, siempre que se cumpla que existe cambio de signo en la función a usar.

    derecha:

    En el extremo derecho, si se propone por ejemplo i con 100%, o 1.00, no tendría mucho sentido un préstamo con intereses al 100% anual, que resulta en el doble del valor inicial en tan solo un periodo o año.

    La tasa de interés de consumo que son de las más alto valor, se encuentran reguladas. En Ecuador es un valor alrededor del 16% anuales o 0.16.

    Considerando las observaciones iniciales del problema, se propone empezar el análisis para la búsqueda de la raíz en el intervalo en un rango más amplio:

    [ 0.01, 0.50]

    Ser realiza la comprobación que existe cambio de signo en los extremos del intervalo.

    fx(0.001) =- 43935.86

    fx(0.50) = 67600.0

    Para el ejercicio se hace notar que la es tasa nominal anual, pero los pagos son mensuales. Por lo que se debe unificar las tasas de interes a mensuales. Una aproximación es usar las tasas anuales divididas para los 12 meses del año.

    Tolerancia al error

    La tolerancia se considera en éste ejercicio como el valor de diferencias  (tramo) entre iteraciones con precisión satisfactoria.

    Por ejemplo si no negociaremos más con el banco por variaciones de tasas del 0.1% , entonces la tolerancia será de 0.001.

    Las publicaciones de tasas en el mercado incluyen dos decimales, por lo que para el ejercicio aumentamos la precisión a : 0.0001

    tolera = 1x10-4


    Literal c


    Se presentan dos formas se solución para el litera c:

    - c.1 la requerida en el enunciado con Newton-Raphson

    - c.2 una alterna con el método de la Bisección.


    c.1. Desarrollo del ejercicio con el método del enunciado Newton-Raphson


    Para el método de Newton-Raphson se tiene que:

    x_{i+1} = x_{i} - \frac{f(x_0i)}{f'(x_i)}

    Se requiere la derivada de la función planteada en el literal a:

    fx(i) = 70000i - 1200(1-(1+i)^{-300}) f'x(i) = 70000 + 1200(300)(1+i)^{-301})

    tomando como valor inicial xi = 0.16/12 ≈ 0.013

    Se realizan las iteraciones suponiendo que tolera = 1x10-4

    iteración 1

    fx(0.013) = 70000(0.013) - 1200(1-(1+0.013)^{-300})

     = -265.09

    f'x(0.013) = 70000 + 1200(300)(1+0.013)^{-301})

    = 62623.3

    x_{2} = 0.013 - \frac{-265.09}{62623.34} = 0.017233

    error = |0.013 - 0.01723| = 0.004331

    iteración 2

    fx(0.01723) = 70000i - 1200(1-(1+0.0.01723)^{-300})

    = 13.446

    f'x(0.01723) = 70000 + 1200(300)(1+0.01723)^{-301}

    = 67897.5

    x_{3} = 0.017233 - \frac{13.446}{67897.5} = 0.017031

    error = |0.017233 - 0.017031| = 0.000198

    cuyo valor de error está casi dentro del valor de tolerancia,

    que permite tomar el último valor como respuesta de tasa mensual

    raiz = tasa mensual = 0.01703

    Convirtiendo a la tasa tasa anual que es la publicada por las instituciones financieras se tiene que:

    tasa anual nominal =  0.01703*12 = 0.2043

    Teniendo como resultado una tasa anual de 20.43%

    E2_IIT2018_T4 Tasa Interes Hipoteca 01


    Algoritmo en Python

    El resultado con el algoritmo es:

    método de Newton-Raphson
    i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
    0 [ 1.30000e-02 -2.65091e+02  6.26233e+04  1.72331e-02  4.23311e-03]
    1 [1.72331e-02 1.34468e+01 6.78975e+04 1.70351e-02 1.98045e-04]
    2 [1.70351e-02 1.24433e-02 6.77706e+04 1.70349e-02 1.83609e-07]
    raiz encontrada en:  0.017034880749732726
    tasa anual:  0.20441856899679273

    Instrucciones en Python

    # 1ra Evaluación II Término 2018
    # Tema 4. Tasa de interes para hipoteca
    import numpy as np
    
    def newton_raphson(fx,dfx,xi, tolera, iteramax=100,
                       vertabla=False, precision=4):
        '''fx y dfx en forma numérica lambda
        xi es el punto inicial de búsqueda
        '''
        itera=0
        tramo = abs(2*tolera)
        if vertabla==True:
            print('método de Newton-Raphson')
            print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
            np.set_printoptions(precision)
        while (tramo>=tolera):
            fi = fx(xi)
            dfi = dfx(xi)
            xnuevo = xi - fi/dfi
            tramo = abs(xnuevo-xi)
            if vertabla==True:
                print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
            xi = xnuevo
            itera = itera + 1
    
        if itera>=iteramax:
            xi = np.nan
            print('itera: ',itera,
                  'No converge,se alcanzó el máximo de iteraciones')
    
        return(xi)
    
    # INGRESO
    P = 70000.00
    A = 1200.00
    n = 25*12
    fx  = lambda i: P*i - A*(1-(1+i)**(-n))
    dfx = lambda i: P + A*(-n)*(i+1)**(-n-1)
    
    x0 = 0.013 # 0.16/12
    tolera = 0.0001
    
    # PROCEDIMIENTO
    raiz   = newton_raphson(fx, dfx, x0, tolera, vertabla=True, precision=5)
    tanual = 12*raiz
    
    # SALIDA
    print('raiz encontrada en: ', raiz)
    print('tasa anual: ',tanual)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    a = 0.01/12
    b = 0.25/12
    muestras = 21
    
    tasa = np.linspace(a,b,muestras)
    fi   = fx(tasa)
    
    plt.plot(tasa*12,fi, label="tasa anual")
    plt.axhline(0, color='green')
    plt.title('tasa anual de interes para Hipoteca')
    plt.xlabel('tasa')
    plt.ylabel('fx(tasa)')
    plt.grid()
    plt.legend()
    plt.show()
    

    c.2. Desarrollo con el método de la Bisección


    Desarrollo Analítico con Bisección

    Como parte del desarrollo del ejercicio se presenta las iteraciones para el algoritmo, tradicionalmente realizadas con una calculadora.

    fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)

    iteración 1

    a = 0.01, b = 0.5 c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 fx(0.01) = 70000 - 1200\Big(\frac{1-(1+(0.01))^{-300}}{0.01} \Big) = -43935.86 fx(0.255) = 70000 - 1200\Big(\frac{1-(1+(0.255))^{-300}}{0.255} \Big) = 65294.11 fx(0.5) = 70000 - 1200\Big(\frac{1-(1+(0.5))^{-300}}{0.5} \Big) = 67600.0 tramo = 0.5-0.01 =0.49

    cambio de signo a la izquierda

    a = 0.01, b=0.255

    iteración 2

    c = \frac{a+b}{2} = \frac{0.01+0.225}{2} = 0.1325 fx(0.1325) = 70000 - 1200\Big(\frac{1-(1+(0.1325))^{-300}}{0.1325} \Big) = 60943.39 tramo = 0.225-0.01 =0.215

    cambio de signo a la izquierda

    a = 0.01, b=0.1325

    iteración 3

    c = \frac{a+b}{2} = \frac{0.01+0.1325}{2} = 0.07125 fx(0.07125) = 70000 - 1200\Big(\frac{1-(1+(0.07125))^{-300}}{0.07125} \Big) = 53157.89 tramo = 0.1325-0.01 =0.1225

    cambio de signo a la izquierda

    a = 0.01, b=0.07125

    y se continuaría con las iteraciones, hasta cumplir que tramo<=tolera

    Tabla de datos obtenidos

    tabla para Bisección
    i a c b f(a) f(c) f(b) tramo
    1 0.01 0.255 0.5 -43935.86 65294.11 67600.0 0.49
    2 0.01 0.1325 0.255 -43935.86 60943.39 65294.11 0.215
    3 0.01 0.07125 0.1325 -43935.86 53157.89 60943.39 0.1225

    hasta lo calculado la raiz se encontraría en el intervalo [0.01,0.07125] con error estImado de 0.06125, aún por mejorar con más iteraciones.

    Algoritmo en Python para Bisección

    • El algoritmo bisección usa las variables a y b, por lo que los limites en el intervalo usados son [La,Lb]
    • para el problema la variable 'i' se usa en el eje x.
    • La selección de cambio de rango [a,b] se hace usando solo el signo del valor.
    • El algoritmo presentado es tal como se explica en la parte conceptual

    Se deja como tarea convertir el algoritmo a funcion def-return de Python.

    # 1Eva_IIT2018_T4 Tasa de interés en hipoteca
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    P = 70000.00
    A = 1200.00
    n = 25*12
    fi = lambda i: P - A*(1-((1+i)**-n))/i
    
    # Intervalo de observación
    # e inicio de Bisección
    La = 0.01
    Lb = 0.50
    
    tolera = 0.0001 #grafica
    
    muestras = 21
    
    # PROCEDIMIENTO
    
    # Método de Bisección
    a = La
    b = Lb
    c = (a+b)/2
    tramo = np.abs(b-a)
    while (tramo>tolera):
        fa = fi(a)
        fb = fi(b)
        fc = fi(c)
        cambio = np.sign(fc)*np.sign(fa)
        if (cambio>0):
            a = c
            b = b
        else:   
            b = c
            a = a
        c = (a+b)/2
        tramo = np.abs(b-a)
    
    # Para la gráfica
    tasa = np.linspace(La,Lb,muestras)
    fr = fi(tasa)
    
    # SALIDA
    print('a, f(a):', a,fa)
    print('c, f(c):', c,fc)
    print('b, f(b):', b,fb)
    print('la raiz esta entre: \n',a,b)
    print('con un error de: ', tramo)
    print('raiz es tasa buscada: ', c)
    print('tasas anual buscada: ',c*12)
    
    # Gráfica
    plt.plot(tasa,fr)
    plt.axhline(0, color='green')
    plt.title('tasa de interes mensual')
    plt.show()
    

    la ejecución del algoritmo da como resultado

    >>> 
     RESTART: D:/MATG1052Ejemplos/HipotecaInteres.py 
    a, f(a): 0.016998291015625 -385.52828922150366
    c, f(c): 0.0170281982421875 -145.85350695741363
    b, f(b): 0.01705810546875 92.28034212642524
    la raiz esta entre: 
     0.016998291015625 0.01705810546875
    con un error de:  5.981445312500111e-05
    raiz es tasa buscada:  0.0170281982421875
    tasas anual buscada:  0.20433837890625
    

    y la gráfica obtenida es:

    tasa interes mensual

  • s1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

    Ejercicio: 1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

    El ejercicio tiene dos partes: la interpolación y el integral.

    Literal a

    t [s] 0 2 4 6 8
    v(t) [m/s] 0.0 16.40 27.77 35.64 41.10

    https://www.dreamstime.com/stock-photo-skydiving-formation-group-people-image62015024No se especifica el método a seguir, por lo que se puede seleccionar el de mayor preferencia.

    Por ejemplo. usando el método de Lagrange, con los puntos primero, medio y último, para cubrir todo el intervalo:

    p_2(t) = 0\frac{(t-4)(t-8)}{(0-4)(0-8)} + + 27.77\frac{(t-0)(t-8)}{(4-0)(4-8)} + + 41.10\frac{(t-0)(t-4)}{(8-0)(8-4)} p_2(t) = 0 + 27.77\frac{t(t-8)}{-16}) + + 41.10\frac{t(t-4)}{32} p_2(t) = -1.73(t^2-8t) + 1.28(t^2-4t) p_2(t) = -0.45 t^2 + 8.74t

    2_IIT2018_T1 Interpola Paracaidista 01


    Literal b

    El tema de integración para primera evaluación se realiza de forma analítica.

    Una de las formas, que es independiente si se resolvió el literal a, es usar los datos proporcionados en la tabla el ejercicio:

    t [s] 0 2 4 6 8
    v(t) [m/s] 0.0 16.40 27.77 35.64 41.10

    Se podría usar el método de Simpson de 1/3, puesto que los tamaños de paso en t son equidistantes se puede aplicar: h=2-0=2

    \int_0^8 v(t)dt = \frac{2}{3}\Big( 0+ 4(16.40)+27.77\Big) + \frac{2}{3}\Big( 27.77+ 4(35.64)+41.10\Big) =203.2

    con error del orden de O(h5) que al considerar h=2 no permite hacer una buena estimación del error. Sin embargo la respuesta es bastante cercana si se usa el método el trapecio con el algoritmo:

        valores de fi:  [ 0.   27.77 41.1 ]
    divisores en L(i):  [ 32. -16.  32.]
    
    Polinomio de Lagrange, expresiones
    -1.735625*x*(x - 8.0) + 1.284375*x*(x - 4.0)
    
    Polinomio de Lagrange: 
    -0.45125*x**2 + 8.7475*x
    Método del trapecio
    distancia recorrida:  193.28
    >>> 
    

    El error entre los métodos es |203.2-193.28|= 9.92

    Revisar el resultado usando un método con mayor precisión que el trapecio.


    Algoritmo con Python

    Las instrucciones en Python para el ejercicio son:

    # 1ra Evaluación II Término 2018
    # Tema 1. Interpolar velocidad del paracaidista
    import numpy as np
    import sympy as sym
    import matplotlib.pyplot as plt
    
    # Literal a)
    # Interpolacion de Lagrange
    # divisoresL solo para mostrar valores
    
    # INGRESO
    t = [0.0, 2, 4, 6, 8]
    v = [0.0, 16.40, 27.77, 35.64, 41.10]
    
    cuales = [0,2,4]
    
    # PROCEDIMIENTO
    xi = np.array(t,dtype=float)
    fi = np.array(v,dtype=float)
    
    xi = xi[cuales]
    fi = fi[cuales]
    
    # 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 = 51
    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(t,v,'o', label = 'Puntos')
    plt.plot(xi,fi,'o', label = 'Puntos en polinomio')
    plt.plot(pxi,pfi, label = 'Polinomio')
    plt.legend()
    plt.xlabel('xi')
    plt.ylabel('fi')
    plt.title('Interpolación Lagrange')
    plt.grid()
    plt.show()
    
    # Literal b
    # INGRESO
    # El ingreso es el polinomio en forma lambda
    # se mantienen las muestras
    
    # intervalo de integración
    # a, b seleccionados para la gráfica anterior
    tramos = muestras -1
    
    # PROCEDIMIENTO
    def integratrapecio_fi(xi,fi):
        ''' sobre muestras de fi para cada xi
            integral con método de trapecio
        '''
        n = len(xi)
        suma = 0
        for i in range(0,n-1,1):
            dx = xi[i+1]-xi[i]
            untrapecio = dx*(fi[i+1]+fi[i])/2
            suma = suma + untrapecio
        return(suma)
    
    
    tramos = muestras-1
    # PROCEDIMIENTO
    distancia = integratrapecio_fi(xi,fi)
    
    # SALIDA
    print('Método del trapecio')
    print('distancia recorrida: ', distancia)
    
  • s1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

    Ejercicio: 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

    Los datos del ejercicio proporcionados son:

    i 0 1 2 3 4 5
    x 1.0 1.1 1.3 1.5 1.9 2.1
    y(x) 1.84 1.90 2.10 2.28 2.91 3.28

    Literal a

    El tema es semejante al tema 1, cambiando el método de interpolación.
    Se usan los puntos de las posiciones 0, 3 y 5.

    p_2(x) = b_0 + b_1x + b_2 x^2

    en la fórmula:

    punto x[0] = 1, y[0]= 1.84

    1.84 = b_0 + b_1(1) + b_2 (1)^2 1.84 = b_0 + b_1 + b_2

    punto x[3] = 1.5, y[3]= 2.28

    2.28 = b_0 + b_1(1.5) + b_2 (1.5)^2 2.28 = b_0 + 1.5 b_1 + 2.25 b_2

    punto x[5] = 2.1, y[5]= 3.28

    3.28= b_0 + b_1(2.1) + b_2 (2.1)^2 3.28= b_0 + 2.1 b_1 + 4.41 b_2

    se obtiene el sistema de ecuaciones:

    b_0 + b_1 + b_2 = 1.84 b_0 + 1.5 b_1 + 2.25 b_2 = 2.28 b_0 + 2.1 b_1 + 4.41 b_2 = 3.28

    Con lo que se plantea la forma Ax=B:

    A = \begin{bmatrix} 1 & 1 & 1\\ 1 & 1.5 & 2.25 \\1 & 2.1 & 4.41 \end{bmatrix} B = \begin{bmatrix} 1.84\\ 2.28 \\ 3.28 \end{bmatrix}

    Matriz Aumentada

    AB = \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 1 & 1.5 & 2.25 & 2.28 \\1 & 2.1 & 4.41 &3.28 \end{bmatrix}

    Pivoteo parcial por filas

    Para el primer pivote no se requieren cambio de filas.
    para el segundo pivote de la diagonal se deben intercambiar la fila segunda con la tercera

    \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 1 & 2.1 & 4.41 &3.28 \\ 1 & 1.5 & 2.25 & 2.28 \end{bmatrix}

    Se aplica eliminación hacia adelante:

    fila = 0, columna=0  pivote = AB[0,0]=1

    factor entre las filas es 1/1=1.

    \begin{bmatrix}1 & 1 & 1 & 1.84 \\ 1-1 & 2.1-1 & 4.41 -1 &3.28 -1.84 \\ 1-1 & 1.5 -1 & 2.25 -1 & 2.28 - 1.84 \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0.5 & 1.25 & 0.44 \end{bmatrix}

    fila =1,  columna=1, pivote=AB[1,1] =1.1

    factor entre filas es 0.5/1.1 = 1/2.2

    \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0.5 -\frac{0.5}{1.1}(1.1)& 1.25 -\frac{0.5}{1.1}(3.41)& 0.44-\frac{0.5}{1.1}(1.44) \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0 & -0.3 & -0.214545 \end{bmatrix}

    aplicando sustitución hacia atrás

    b2 = -0.21/(-0.3) = 0.71515 b1= \frac{1.44-3.41 b_2}{1.1} = \frac{1.44-3.41( 0.71515)}{1.1}=-0.9078 b3= \frac{1.84-b_1-b_2}{1} = \frac{1.84-(-0.9078)-(0.71515)}{1} =2.0327

    con lo que el polinomio buscado es:

    p_2(x) = 2.0327 -0.9078 x + 0.71515 x^2

    y se obtiene el resultado de la interpolación.

    E2_IIT2018_T3 Interpola Sistema Ecuaciones 01Observación: En la gráfica se muestra que el polinomio pasa por los puntos seleccionados de la tabla. En los otros puntos hay un error que se puede calcular como la resta del punto y su valor con p(x). Queda como tarea.

    Usando el algoritmo del polinomio de interpolación con la matriz de Vandermonde se obtiene:

    Matriz Vandermonde: 
    [[1.   1.   1.  ]
     [2.25 1.5  1.  ]
     [4.41 2.1  1.  ]]
    los coeficientes del polinomio: 
    [ 0.71515152 -0.90787879  2.03272727]
    Polinomio de interpolación: 
    0.715151515151516*x**2 - 0.907878787878792*x + 2.03272727272728
    
     formato pprint
                       2                                         
    0.715151515151516*x  - 0.907878787878792*x + 2.03272727272728
    suma de columnas:  [3.   4.75 7.51]
    norma D:  7.51
    numero de condicion:  97.03737354737122
    solucion: 
    [ 0.71515152 -0.90787879  2.03272727]
    >>> 
    
    

    Literal b

    Se requiere calcular una norma de suma de filas. es suficiente para demostrar el conocimiento del concepto el usar A.

    Se adjunta el cálculo del número de condición y la solución al sistema de ecuaciones:

    suma de columnas:  [3.   4.75 7.51]
    norma A:  7.51
    numero de condición:  97.03737354737129
    solución: 
    [ 2.03272727 -0.90787879  0.71515152]
    

    El comentario importante corresponde al número de condición, que es un número muy alto para usar un método iterativo, por lo que la solución debe ser un método directo.
    Se puede estimar será un número mucho mayor que 1, pues la matriz no es diagonal dominante.


    Instrucciones en Python

    # 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones
    # El polinomio de interpolación
    import numpy as np
    import sympy as sym
    import matplotlib.pyplot as plt
    
    # INGRESO
    xj = [1.0,  1.1,  1.3,  1.5,  1.9,  2.1 ]
    yj = [1.84, 1.90, 2.10, 2.28, 2.91, 3.28]
    cuales = [0, 3, 5]
    
    # muestras = tramos+1
    muestras = 51
    
    # PROCEDIMIENTO
    # Convierte a arreglos numpy 
    xi = np.array(xj,dtype=float)
    fi = np.array(yj,dtype=float)
    
    xi = xi[cuales]
    fi  = fi[cuales]
    B = fi
    
    n = len(xi)
    
    # Matriz Vandermonde D
    D = np.zeros(shape=(n,n),dtype =float)
    for i in range(0,n,1):
        for j in range(0,n,1):
            potencia = (n-1)-j # Derecha a izquierda
            D[i,j] = xi[i]**potencia
    
    # Aplicar métodos Unidad03. Tarea
    # Resuelve sistema de ecuaciones A.X=B
    coeficiente = np.linalg.solve(D,B)
    
    # Polinomio en forma simbólica
    x = sym.Symbol('x')
    polinomio = 0
    for i in range(0,n,1):
        potencia = (n-1)-i   # Derecha a izquierda
        termino = coeficiente[i]*(x**potencia)
        polinomio = polinomio + termino
    
    # Polinomio a forma Lambda
    # para evaluación con vectores de datos xin
    px = sym.lambdify(x,polinomio)
    
    # Para graficar el polinomio en [a,b]
    a = np.min(xi)
    b = np.max(xi)
    xin = np.linspace(a,b,muestras)
    yin = px(xin)
    
    # Usando evaluación simbólica
    ##yin = np.zeros(muestras,dtype=float)
    ##for j in range(0,muestras,1):
    ##    yin[j] = polinomio.subs(x,xin[j])
        
    # SALIDA
    print('Matriz Vandermonde: ')
    print(D)
    print('los coeficientes del polinomio: ')
    print(coeficiente)
    print('Polinomio de interpolación: ')
    print(polinomio)
    print('\n formato pprint')
    sym.pprint(polinomio)
    
    # Grafica
    plt.plot(xj,yj,'o', label='Puntos')
    plt.plot(xi,B,'o', label='[xi,fi]')
    plt.plot(xin,yin, label='p(x)')
    plt.xlabel('xi')
    plt.ylabel('fi')
    plt.legend()
    plt.title(polinomio)
    plt.show()
    
    
    # literal b
    sumacolumnas = np.sum(D, axis =1)
    norma = np.max(sumacolumnas)
    print('suma de columnas: ', sumacolumnas)
    print('norma D: ', norma)
    
    numerocondicion = np.linalg.cond(D)
    print('numero de condicion: ', numerocondicion)
    
    solucion = np.linalg.solve(D,B)
    print('solucion: ')
    print(solucion)
    

    .

  • s1Eva_IIT2018_T2 Distancia mínima a un punto

    Ejercicio: 1Eva_IIT2018_T2 Distancia mínima a un punto

    Literal a

    Se requiere analizar la distancias entre una trayectoria y el punto = [1,1]

    Al analizar las distancias de ex y el punto [1,1] se trazan lineas paralelas a los ejes desde el punto [1,1], por lo que se determina que el intervalo de x = [a,b] para distancias se encuentra en:

    a > 0, a = 0.1
    b < 1, b = 0.7

    El ejercicio usa la fórmula de distancia entre dos puntos:

    d = \sqrt{(x_2-x_1)^2+(y_2- y_1)^2}

    en los cuales:

    [x1,y1] = [1,1]
    [x2,y2] = [x, ex]

    que al sustituir en la fórmula se convierte en:

    d = \sqrt{(x-1)^2+(e^x- 1)^2}

    que es lo requerido en el literal a


    Literal b

    Para usar un método de búsqueda de raíces, se requiere encontrar el valor cuando f(x) = d' = 0.

    Un método como el de Newton Raphson requiere también f'(x) = d''

    f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}} f'(x)= \frac{(e^x - 1)e^x + e^{2x} + 1 - \frac{(x + (e^x - 1)e^x - 1)^2}{(x - 1)^2 + (e^x - 1)^2}} {\sqrt{(x - 1)^2 + (e^x - 1)^2}}

    expresiones obtenidas usando Sympy

    f(x) :
    (x + (exp(x) - 1)*exp(x) - 1)/sqrt((x - 1)**2 + (exp(x) - 1)**2)
    f'(x) :
    ((exp(x) - 1)*exp(x) + exp(2*x) + 1 - (x + (exp(x) - 1)*exp(x) - 1)**2/((x - 1)**2 + (exp(x) - 1)**2))/sqrt((x - 1)**2 + (exp(x) - 1)**2)
    
    f(x) :
           / x    \  x        
       x + \e  - 1/*e  - 1    
    --------------------------
        ______________________
       /                    2 
      /         2   / x    \  
    \/   (x - 1)  + \e  - 1/  
    f'(x) :
                                                  2
                             /    / x    \  x    \ 
    / x    \  x    2*x       \x + \e  - 1/*e  - 1/ 
    \e  - 1/*e  + e    + 1 - ----------------------
                                                 2 
                                     2   / x    \  
                              (x - 1)  + \e  - 1/  
    -----------------------------------------------
                   ______________________          
                  /                    2           
                 /         2   / x    \            
               \/   (x - 1)  + \e  - 1/            
    
    
    

    lo que permite observar la raíz de f(x) en una gráfica:
    distancia mínima f(x)
    con las siguientes instrucciones:

    # Eva_IIT2018_T2 Distancia mínima a un punto
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    # INGRESO
    x = sym.Symbol('x')
    fx = sym.sqrt((x-1)**2+(sym.exp(x) -1)**2)
    
    a = 0
    b = 1
    muestras = 21
    
    # PROCEDIMIENTO
    dfx = sym.diff(fx,x,1)
    d2fx = sym.diff(fx,x,2)
    
    f = sym.lambdify(x,dfx)
    xi = np.linspace(a,b,muestras)
    fi = f(xi)
    
    
    # SALIDA
    print('f(x) :')
    print(dfx)
    print("f'(x) :")
    print(d2fx)
    print()
    print('f(x) :')
    sym.pprint(dfx)
    print("f'(x) :")
    sym.pprint(d2fx)
    
    # GRAFICA
    plt.plot(xi,fi, label='f(x)')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.legend()
    plt.grid()
    plt.show()
    

    Usando el método de la bisección para el intervalo dado, se tiene:

    f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}}

    itera = 0 , a = 0, b=1

    c= \frac{0+1}{2} = 0.5 f(0) = \frac{0 + (e^0 - 1)e^0 - 1}{\sqrt{(0 - 1)^2 + (e^0 - 1)^2}} = -1 f(1) = \frac{1 + (e^1 - 1)e^1 - 1}{\sqrt{(1 - 1)^2 + (e^1 - 1)^2}} 2.7183 f(0.5) = \frac{(0.5) + (e^(0.5) - 1)e^(0.5) - 1}{\sqrt{((0.5) - 1)^2 + (e^(0.5) - 1)^2}} = 0.6954

    cambio de signo a la izquierda,

    a= 0, b=c=0.5

    tramo = |0.5-0| = 0.5

    itera = 1

    c= \frac{0+0.5}{2} = 0.25 f(0.25) = \frac{(0.25) + (e^(0.25) - 1)e^(0.25) - 1}{\sqrt{((0.25) - 1)^2 + (e^(0.25) - 1)^2}} = -0.4804

    cambio de signo a la derecha,

    a=c= 0.25, b=0.5

    itera = 2

    c= \frac{0.25+0.5}{2} = 0.375 f(0.375) = \frac{(0.375) + (e^(0.375) - 1)e^(0.375) - 1}{\sqrt{((0.375) - 1)^2 + (e^(0.375) - 1)^2}} = 0.0479

    cambio de signo a la izquierda,

    a= 0.25, b=c=0.375

    se continúan las iteraciones con el algoritmo, para encontrar la raíz en 0.364:

    método de Bisección
    i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
       tramo
    0 [0, 0.5, 1] [-1.      0.6954  2.7183]
       0.5
    1 [0, 0.25, 0.5] [-1.     -0.4804  0.6954]
       0.25
    2 [0.25, 0.375, 0.5] [-0.4804  0.0479  0.6954]
       0.125
    3 [0.25, 0.3125, 0.375] [-0.4804 -0.2388  0.0479]
       0.0625
    4 [0.3125, 0.34375, 0.375] [-0.2388 -0.1004  0.0479]
       0.03125
    5 [0.34375, 0.359375, 0.375] [-0.1004 -0.0274  0.0479]
       0.015625
    6 [0.359375, 0.3671875, 0.375] [-0.0274  0.01    0.0479]
       0.0078125
    7 [0.359375, 0.36328125, 0.3671875] [-0.0274 -0.0088  0.01  ]
       0.00390625
    8 [0.36328125, 0.365234375, 0.3671875] [-0.0088  0.0006  0.01  ]
       0.001953125
    9 [0.36328125, 0.3642578125, 0.365234375] [-0.0088 -0.0041  0.0006]
       0.0009765625
    raíz en:  0.3642578125
    

    Al algoritmo anterior se complementa con las instrucciones de la función para la bisección.

    # Eva_IIT2018_T2 Distancia mínima a un punto
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    # INGRESO
    x = sym.Symbol('x')
    fx = sym.sqrt((x-1)**2+(sym.exp(x) -1)**2)
    
    a = 0
    b = 1
    muestras = 21
    
    # PROCEDIMIENTO
    dfx = sym.diff(fx,x,1)
    d2fx = sym.diff(fx,x,2)
    
    f = sym.lambdify(x,dfx)
    xi = np.linspace(a,b,muestras)
    fi = f(xi)
    
    
    # SALIDA
    print('f(x) :')
    print(dfx)
    print("f'(x) :")
    print(d2fx)
    print()
    print('f(x) :')
    sym.pprint(dfx)
    print("f'(x) :")
    sym.pprint(d2fx)
    
    # GRAFICA
    plt.plot(xi,fi, label='f(x)')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.legend()
    plt.grid()
    plt.show()
    
    # Algoritmo de Bisección
    # [a,b] se escogen de la gráfica de la función
    # error = tolera
    import numpy as np
    
    def biseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
        '''
        Algoritmo de Bisección
        Los valores de [a,b] son seleccionados
        desde la gráfica de la función
        error = tolera
        '''
        fa = fx(a)
        fb = fx(b)
        tramo = np.abs(b-a)
        itera = 0
        cambia = np.sign(fa)*np.sign(fb)
        if cambia<0: # existe cambio de signo f(a) vs f(b)
            if vertabla==True:
                print('método de Bisección')
                print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
                print('  ','tramo')
                np.set_printoptions(precision)
                
            while (tramo>=tolera and itera<=iteramax):
                c = (a+b)/2
                fc = fx(c)
                cambia = np.sign(fa)*np.sign(fc)
                if vertabla==True:
                    print(itera,[a,c,b],np.array([fa,fc,fb]))
                if (cambia<0):
                    b = c
                    fb = fc
                else:
                    a = c
                    fa = fc
                tramo = np.abs(b-a)
                if vertabla==True:
                    print('  ',tramo)
                itera = itera + 1
            respuesta = c
            # Valida respuesta
            if (itera>=iteramax):
                respuesta = np.nan
    
        else: 
            print(' No existe cambio de signo entre f(a) y f(b)')
            print(' f(a) =',fa,',  f(b) =',fb) 
            respuesta=np.nan
        return(respuesta)
    
    # INGRESO
    tolera = 0.001
    
    # PROCEDIMIENTO
    respuesta = biseccion(f,a,b,tolera,vertabla=True)
    # SALIDA
    print('raíz en: ', respuesta)
    

    .

  • Protegido: s3Eva_IT2018_T3 EDP Parabólica, temperatura en varilla

    Este contenido está protegido por contraseña. Para verlo introduce tu contraseña a continuación:

  • s3Eva_IT2018_T2 Drenaje de estanque

    Ejercicio: 3Eva_IT2018_T2 Drenaje de estanque

    literal a

    Se usa interpolación para encontrar los polinomios que pasan por los puntos seleccionados.

    El error de A(5) se obtiene como la diferencia entre el valor de la tabla y el polinomio del tramo [4,6] evaluado en el punto.

    ordenado:  [6 5 4 3 2 1 0]
    hi:  [0 1 2 3 4 5 6]
    Ai:  [ 0.02  0.18  0.32  0.45  0.67  0.97  1.17]
    
    puntos seleccionados:
    h1:  [0, 2, 4, 6]
    A1:  [ 0.02  0.32  0.67  1.17]
    
    Polinomios por tramos: 
     x = [0,2]
    0.000416666666666669*x**3 + 0.148333333333333*x + 0.02
     x = [2,4]
    0.00416666666666666*x**3 - 0.0224999999999999*x**2 + 0.193333333333333*x - 0.00999999999999984
     x = [4,6]
    -0.00458333333333333*x**3 + 0.0824999999999999*x**2 - 0.226666666666666*x + 0.549999999999999
    
    error en px(5):  0.0637499999999998
    

    se observa que la evaluación se realiza para el polinomio entre [4,6]

    drenaje estanque 02 Ah

    Desarrollo en Python

    # 3ra Evaluación I Término 2018
    # Tema 2. Drenaje de Estanque
    
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    def traza3natural(xi,yi):
        # Trazador cúbico natural, splines
        # resultado: polinomio en forma simbólica
        n = len(xi)
        # Valores h
        h = np.zeros(n-1, dtype = float)
        for j in range(0,n-1,1):
            h[j] = xi[j+1] - xi[j]
        
        # Sistema de ecuaciones
        A = np.zeros(shape=(n-2,n-2), dtype = float)
        B = np.zeros(n-2, dtype = float)
        S = np.zeros(n, dtype = float)
        A[0,0] = 2*(h[0]+h[1])
        A[0,1] = h[1]
        B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])
        for i in range(1,n-3,1):
            A[i,i-1] = h[i]
            A[i,i] = 2*(h[i]+h[i+1])
            A[i,i+1] = h[i+1]
            B[i] = 6*((yi[i+2]-yi[i+1])/h[i+1] - (yi[i+1]-yi[i])/h[i])
        A[n-3,n-4] = h[n-3]
        A[n-3,n-3] = 2*(h[n-3]+h[n-2])
        B[n-3] = 6*((yi[n-1]-yi[n-2])/h[n-2] - (yi[n-2]-yi[n-3])/h[n-3])
        
        # Resolver sistema de ecuaciones
        r = np.linalg.solve(A,B)
        # S
        for j in range(1,n-1,1):
            S[j] = r[j-1]
        S[0] = 0
        S[n-1] = 0
        
        # Coeficientes
        a = np.zeros(n-1, dtype = float)
        b = np.zeros(n-1, dtype = float)
        c = np.zeros(n-1, dtype = float)
        d = np.zeros(n-1, dtype = float)
        for j in range(0,n-1,1):
            a[j] = (S[j+1]-S[j])/(6*h[j])
            b[j] = S[j]/2
            c[j] = (yi[j+1]-yi[j])/h[j] - (2*h[j]*S[j]+h[j]*S[j+1])/6
            d[j] = yi[j]
        
        # Polinomio trazador
        x = sym.Symbol('x')
        polinomio = []
        for j in range(0,n-1,1):
            ptramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2 + c[j]*(x-xi[j])+ d[j]
            ptramo = ptramo.expand()
            polinomio.append(ptramo)
        
        return(polinomio)
    
    # PROGRAMA -------------------------
    
    hi = np.array([6, 5, 4, 3, 2, 1, 0])
    Ai = np.array([1.17, 0.97, 0.67, 0.45, 0.32, 0.18, 0.02])
    xk = 5
    
    # PROCEDIMIENTO LITERAL a
    # reordena en forma ascendente
    ordenado = np.argsort(hi)
    hi = hi[ordenado]
    Ai = Ai[ordenado]
    
    # Selecciona puntos
    xi = [0,2,4,6]
    fi = Ai[xi]
    n = len(xi)
    
    polinomio = traza3natural(xi,fi)
    
    # literal a, estima error
    px = polinomio[2]
    pxk = px.subs('x',xk)
    errado = np.abs(Ai[xk] - pxk)
    
    # SALIDA
    print('ordenado: ', ordenado)
    print('hi: ', hi)
    print('Ai: ', Ai)
    print('puntos seleccionados:')
    print('h1: ', xi)
    print('A1: ', fi)
    
    print('Polinomios por tramos: ')
    for tramo in range(1,n,1):
        print(' x = ['+str(xi[tramo-1])+','+str(xi[tramo])+']')
        print(str(polinomio[tramo-1]))
    
    print('error en px(5): ', errado)
    
    # GRAFICA
    # Puntos para grafica en cada tramo
    resolucion = 10 # entre cada par de puntos
    xtrazado = np.array([])
    ytrazado = np.array([])
    tramo = 1
    while not(tramo>=n):
        a = xi[tramo-1]
        b = xi[tramo]
        xtramo = np.linspace(a,b,resolucion)
        
        ptramo = polinomio[tramo-1]
        pxtramo = sym.lambdify('x',ptramo)
        ytramo = pxtramo(xtramo)
        
        xtrazado = np.concatenate((xtrazado,xtramo))
        ytrazado = np.concatenate((ytrazado,ytramo))
        tramo = tramo + 1
    
    # GRAFICA
    # puntos originales
    plt.plot(hi,Ai,'o',label = 'Ai')
    # Trazador cúbico
    plt.plot(xtrazado,ytrazado, label = 'p(h)')
    plt.plot(xi,fi,'o', label = 'Apx')
    plt.title('Trazador cúbico natural (splines)')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()
    

    Literal b

    TAREA

  • 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]]
    >>>