Etiqueta: algoritmo Python

  • 6.2 EDO dy/dx, Runge-Kutta con Python

    [ Runge Kutta  dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    1. EDO \frac{\delta y}{\delta x} con Runge-Kutta de 2do Orden

    Referencia: Burden 5.4 p209, Chapra 25.3 p740, Rodríguez 9.1.7 p354, Boyce DiPrima 4Ed 8.4 p450

    Para una ecuación diferencial ordinaria de primera derivada, el método Runge-Kutta de 2do Orden usa una corrección sobre la derivada a partir de los puntos xi y xi+h,  es decir un tamaño de paso h hacia adelante, calculados como términos K1 y K2.

    EDO Runge-Kutta 2do orden primera derivada _animado
    Considere una ecuación diferencial de primera derivada con una condición de inicio se reordena y se escribe como f(x,y) siguiendo los pasos:

    \frac{\delta y}{\delta x} + etc =0 y'(x) = f(x_i,y_i) y(x_0) = y_0

    Los términos K1 y K2 se calculan para predecir el próximo valor en y[i+1], observe que el término K1 es el mismo que el método de Edo con Taylor de dos términos.

    K_1 = h f(x_i,y_i) K_2 = h f(x_i+h, y_i + K_1) y_{i+1} = y_i + \frac{K_1+K_2}{2} x_{i+1} = x_i + h

    Runge-Kutta 2do Orden 02

    Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

    Las iteraciones se repiten para encontrar el siguiente punto en x[i+1] como se muestra en el gráfico animado.

    Los métodos de Runge-Kutta  mejoran la aproximación a la respuesta de la ecuación diferencial ordinaria sin requerir determinar las expresiones de las derivadas de orden superior, como fue necesario en EDO con Taylor.

    Para observar al idea básica, considere observar un combate aéreo simulado en la década de 1940, donde las armas se encuentras fijas en las alas del avión. Observe dos minutos del video sugerido a partir de donde se encuentra marcado el enlace.

    Video Revisar:

    Luego de observar el video de introducción conteste las siguientes preguntas:
    ¿ Que trayectoria siguen los proyectiles al salir del cañón?
    ¿ Que trayectoria siguen los aviones, el perseguido y el que caza?
    ¿ Cuál es la relación entre las trayectorias de los dos aviones?

    Runge-Kutta 2do Orden primera derivada

    [ Runge Kutta  dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    2. Ejercicio

    Referencia: Rodríguez 9.1.1 ejemplo p335. Chapra 25.1.3 p731

    Se pide encontrar puntos de la solución en la ecuación diferencial usando los tres primeros términos de la serie de Taylor con h=0.1 y punto inicial x0=0, y0=1

    \frac{dy}{dx}-y -x +x^2 -1 = 0 y'-y -x +x^2 -1 = 0

    [ Runge Kutta  dy/dx] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    3. Desarrollo Analítico

    Se reordena la expresión haciendo que la derivada se encuentre en el lado izquierdo:

    f(x,y) = y' = y -x^2 +x +1

    Se usa las expresiones de Runge-Kutta en orden, K1 corresponde a una corrección de EDO con Taylor de dos términos (método de Euler). K2 considera el cálculo a un tamaño de paso más adelante. iteración:

    K_1 = h f(x_i,y_i) = 0.1 (y_i -x_i^2 +x_i +1) K_2 = h f(x_i+h, y_i + K_1) K_2 = 0.1 \Big((y_i+K_1) -(x_i+h)^2 +(x_i+h) +1 \Big) y_{i+1} = y_i + \frac{K_1+K_2}{2} x_{i+1} = x_i + h

    Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

    EDO Runge-Kutta 2do orden primera derivada _animado

    para el ejercicio, el tamaño de paso h=0.1, se realizan tres iteraciones en las actividades del curso con lápiz y papel,

    itera = 0 , x0 = 0, y0 = 1

    K_1 = 0.1 f(0,1) = 0.1 \Big( 1 -0^2 +0 +1 \Big) = 0.2 K_2 = 0.1 f(0+0.1, 1+ 0.2) K_2 = 0.1 \Big( (1+ 0.2) - (0+0.1) ^2 +(0+0.1) +1\Big) = 0.229 y_1 = 1 + \frac{0.2+0.229}{2} = 1.2145 x_1 = 0 + 0.1 = 0.1

    itera = 1 , x1 = 0.1, y1 = 1.2145

    K_1 = 0.1 f(0.1,1.2145) = 0.1( 1.2145 -0.1^2 +0.1 +1) K_1 = 0.2304 K_2 = 0.1 f(0.1+0.1, 1.2145 + 0.2304) =0.1 \Big((1.2145 + 0.2304) -(0.1+0.1)^2 +(0.1+0.1) +1\Big) K_2 = 0.2604 y_2 = 1.2145 + \frac{0.2304+0.2604}{2} = 1.4599 x_2 = 0.1 +0.1 = 0.2

    itera = 2 , x2 = 0.2, y2 = 1.4599

    K_1 = 0.1 f(0.2,1.4599) = 0.1( 1.4599 -0.2^2 +0.2 +1) K_1 = 0.2619 K_2 = 0.1 f(0.2+0.1, 1.4599 + 0.2619) =0.1 \Big((1.4599 + 0.2619) -(0.2+0.1)^2 +(0.2+0.1) +1\Big) K_2 = 0.2931 y_2 = 1.4599 + \frac{0.2619+0.2931}{2} = 1.7375 x_2 = 0.2 +0.1 = 0.3

    luego de las 3 iteraciones en papel, se completan los demás puntos con el algoritmo obteniendo la gráfica resultante para y(x) correspondiente.

    EDO dy/dx con Runge-Kutta 2 Orden
    i, [xi,     yi,     K1,    K2]
    0 [0. 1. 0. 0.]
    1 [0.1    1.2145 0.2    0.229 ]
    2 [0.2      1.459973 0.23045  0.260495]
    3 [0.3      1.73757  0.261997 0.293197]
    4 [0.4      2.048564 0.294757 0.327233]
    5 [0.5      2.394364 0.328856 0.362742]
    >>> 
    

    ecuación diferencial ordinaria con Runge-Kutta de 2do orden

    Compare los resultados con Taylor de 2 y 3 términos.

    Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

    [ Runge Kutta  dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    4. Algoritmo en Python

    Para el ejercicio anterior, se adjunta el programa de prueba que usa la función rungekutta2(d1y,x0,y0,h,muestras) .

    Las iteraciones y sus valores se pueden observar usando vertabla=true

    # EDO dy/dx. M todo de RungeKutta 2do Orden 
    # estima la solucion para muestras espaciadas h en eje x
    # valores iniciales x0,y0, entrega tabla[xi,yi,K1,K2]
    import numpy as np
    
    def rungekutta2(d1y,x0,y0,h,muestras,
                    vertabla=False,precision=6):
        '''solucion a EDO dy/dx, con Runge Kutta de 2do orden
        d1y es la expresi n dy/dx, tambien planteada como f(x,y),
        valores iniciales: x0,y0, tama o de paso h.
        muestras es la cantidad de puntos a calcular. 
        '''
        tamano = muestras + 1
        tabla = np.zeros(shape=(tamano,2+2),dtype=float)
        tabla[0] = [x0,y0,0,0] # incluye el punto [x0,y0]
        
        xi = x0 # valores iniciales
        yi = y0
        for i in range(1,tamano,1):
            K1 = h * d1y(xi,yi)
            K2 = h * d1y(xi+h, yi + K1)
    
            yi = yi + (K1+K2)/2
            xi = xi + h
            
            tabla[i] = [xi,yi,K1,K2]
           
        if vertabla==True:
            np.set_printoptions(precision)
            print( 'EDO dy/dx con Runge-Kutta 2 Orden')
            print('i, [xi,     yi,     K1,    K2]')
            for i in range(0,tamano,1):
                print(i,tabla[i])
    
        return(tabla)
    
    # PROGRAMA PRUEBA -------------------
    
    # INGRESO
    # d1y = y' = f
    d1y = lambda x,y: y -x**2 + x + 1
    x0 = 0
    y0 = 1
    h  = 0.1
    muestras = 5
    
    # PROCEDIMIENTO
    tabla = rungekutta2(d1y,x0,y0,h,muestras,
                        vertabla=True)
    xi = tabla[:,0]
    yiRK2 = tabla[:,1]
    
    # SALIDA
    #print(' [xi, yi, K1, K2]')
    #print(tabla)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    plt.plot(xi,yiRK2)
    plt.plot(xi[0],yiRK2[0],
             'o',color='r', label ='[x0,y0]')
    plt.plot(xi[1:],yiRK2[1:],
             'o',color='m',
             label ='[xi,yi] Runge-Kutta')
    plt.title('EDO dy/dx: Runge-Kutta 2do Orden')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()
    

    [ Runge Kutta  dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]


    5. Cálculo de Error con la solución conocida

    La ecuación diferencial ordinaria del ejercicio tiene una solución conocida, lo que permite encontrar el error real en cada punto respecto a la aproximación estimada.

    y = e^x + x + x^2

    EDO RungeKutta 2Orden 01

    Note que el error crece al distanciarse del punto inicial.

    Error máximo estimado:  0.004357584597315167
    entre puntos: 
    [0.         0.00067092 0.00143026 
     0.0022892  0.00326028 0.00435758]
    >>>
    

    Para las siguientes instrucciones, comente la última línea #plt.show() antes de continuar con:

    # ERROR vs solución conocida -----------------
    y_sol = lambda x: ((np.e)**x) + x + x**2
    
    yi_psol  = y_sol(xi)
    errores  = yi_psol - yiRK2
    errormax = np.max(np.abs(errores))
    
    # SALIDA
    print('Error máximo estimado: ',errormax)
    print('entre puntos: ')
    print(errores)
    
    # GRAFICA [a,b+2*h]
    a = x0
    b = h*muestras+2*h
    muestreo = 10*muestras+2
    xis = np.linspace(a,b,muestreo)
    yis = y_sol(xis)
    
    plt.plot(xis,yis, label='y solución conocida',
             linestyle='dashed')
    plt.legend()
    plt.show()
    

    [ Runge Kutta  dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]


    6. Ejercicio de Evaluación

    2Eva_IT2018_T1 EDO Paracaidista wingsuit

    Solución Propuesta: s2Eva_IT2018_T1 EDO Paracaidista wingsuit , Runge-Rutta para primera derivada.

  • 6.1 EDO con Taylor de 3 términos con Python

    [ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    1. Ecuaciones diferenciales ordinarias aproximadas con Taylor

    Referencia: Rodríguez 9.1.1 ejemplo p335. Chapra 25.1.3 p731

    En los métodos con Taylor para Ecuaciones Diferenciales Ordinarias (EDO) se aproxima el resultado a n términos de la serie, para lo cual se ajusta la expresión del problema a cada derivada correspondiente.

    La solución empieza usando la Serie de Taylor para tres términos ajustada a la variable del ejercicio:

    y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h E = \frac{h^3}{3!} y'''(z) = O(h^3)

    Edo Taylor 3 términos GIF animado
    A partir de la expresión de y'(x) y el punto inicial conocido en x[i],se busca obtener el próximo valor en x[i+1] al avanzar un tamaño de paso h. Se repite el proceso en el siguiente punto encontrado y se continua hasta alcanzar el intervalo objetivo.

    EDO Taylor 3 terminos

    En éstos métodos la solución siempre es una tabla de puntos xi,yi que se pueden usar para interpolar y obtener una función polinómica.

    [ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

    ..


    2. Ejercicio

    Referencia: Rodríguez 9.1.1 ejemplo p335. Chapra 25.1.3 p731

    Se requiere encontrar puntos de la solución en la ecuación diferencial usando los tres primeros términos de la serie de Taylor con h=0.1 y punto inicial x0=0, y0=1

    \frac{dy}{dx}-y -x +x^2 -1 = 0

    que con nomenclatura simplificada:

    y'-y -x +x^2 -1 = 0

    [ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    3. Desarrollo Analítico

    Al despejar el valor de  y' de expresión del ejercicio,

    y' = y -x^2 +x +1

    se puede obtener y" al derivar una vez,

    y'' = y' -2x + 1

    para luego combinar las expresiones en

    y'' = (y -x^2 +x +1) -2x + 1

    simplificando:

    y'' = y -x^2 -x +2

    Ecuaciones que permiten estimar nuevos valores yi+1 para nuevos puntos  muestra distanciados en i*h desde el punto inicial siguiendo las siguientes expresiones de iteración:

    y'_i = y_i -x_i^2 + x_i +1 y''_i = y_i -x_i^2 - x_i +2 y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h

    Edo Taylor 3 términos GIF animado

    Se empieza evaluando el nuevo punto a una distancia x1= x0+h del punto de origen con lo que se obtiene y1 , repitiendo el proceso para el siguiente punto en forma sucesiva.

    itera = 0 , x0 = 0, y0 = 1

    y'_0 = 1 -0^2 +0 +1 = 2 y''_0 = 1 -0^2 -0 +2 = 3 y_1 = y_{0} + h y'_0 + \frac{h^2}{2!} y''_0 y_1 = 1 + 0.1 (2) + \frac{0.1^2}{2!} 3 = 1.215 x_1 = 0 + 0.1

    itera = 1 , x = 0.1, y = 1.215

    y'_1 = 1.215 - 0.1^2 + 0.1 +1 = 2.305 y''_1 = 1.215 - 0.1^2 - 0.1 +2 = 3.105 y_2 = 1.215 + 0.1 (2.305) + \frac{0.1^2}{2!} 3.105 = 1.461 x_2 = 0.1 + 0.1 = 0.2

    itera = 2 , x = 0.2, y = 1.461

    y'_2 = 1.461 - 0.2^2 + 0.2 +1 = 2.621 y''_2 = 1.461 - 0.2^2 - 0.2 +2 = 3.221 y_3 = 1.461 + 0.1 (2.621) + \frac{0.1^2}{2!} 3.221 = 1.7392 x_3 = 0.2 + 0.1 = 0.3

    completando los puntos con el algoritmo y realizando la gráfica se obtiene

     EDO con Taylor 3 términos
    [xi, yi, d1yi, d2yi]
    [[0.         1.         0.         0.        ]
     [0.1        1.215      2.         3.        ]
     [0.2        1.461025   2.305      3.105     ]
     [0.3        1.73923262 2.621025   3.221025  ]
     [0.4        2.05090205 2.94923262 3.34923262]
     [0.5        2.39744677 3.29090205 3.49090205]]
    >>>

    Observación, note que los resultados de las derivadas, se encuentran desplazados una fila para cada iteración. Asunto a ser considerado en la gráfica de las derivadas en caso de incluirlas.

    EDO_Taylor_3terminos01

    Nota: Compare luego los pasos del algoritmo con el método de Runge-Kutta de 2do orden.

    [ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
    ..


    4. Algoritmo en Python

    Para simplificar los cálculos se crea una función edo_taylor3t() para encontrar  los valores para una cantidad de muestras distanciadas entre si h veces del punto inicial [x0,y0]

    # EDO. Método de Taylor con3 términos 
    # estima solucion para muestras separadas h en eje x
    # valores iniciales x0,y0
    import numpy as np
    
    def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
        ''' solucion a EDO usando tres términos de Taylor, x0,y0 son valores iniciales
            muestras es la cantidad de puntos a calcular con tamaño de paso h.
        '''
        tamano = muestras + 1
        tabla = np.zeros(shape=(tamano,4),dtype=float)
        # incluye el punto [x0,y0]
        tabla[0] = [x0,y0,0,0]
        x = x0
        y = y0
        for i in range(1,tamano,1):
            d1yi = d1y(x,y)
            d2yi = d2y(x,y)
            y = y + h*d1yi + ((h**2)/2)*d2yi
            x = x + h
            tabla[i] = [x,y,d1yi,d2yi]
        return(tabla)
    
    # PROGRAMA PRUEBA -----------------
    # Ref Rodriguez 9.1.1 p335 ejemplo.
    # prueba y'-y-x+(x**2)-1 =0, y(0)=1
    
    # INGRESO.
    # d1y = y', d2y = y''
    d1y = lambda x,y: y - x**2 + x + 1
    d2y = lambda x,y: y - x**2 - x + 2
    x0 = 0
    y0 = 1
    h = 0.1
    muestras = 5
    
    # PROCEDIMIENTO
    tabla = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
    xi = tabla[:,0]
    yi = tabla[:,1]
    
    # SALIDA
    print(' EDO con Taylor 3 términos')
    print('[xi, yi, d1yi, d2yi]')
    print(tabla)
    
    # Gráfica
    import matplotlib.pyplot as plt
    plt.plot(xi,yi)
    plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
    plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
    plt.title('EDO: Solución con Taylor 3 términos')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show() # plt.show() #comentar para la siguiente gráfica
    

    Tarea: Realizar el ejercicio con más puntos muestra, donde se visualice que el error aumenta al aumentar la distancia del punto inicial [x0,y0]

    [ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]


    5. Cálculo de Error con la solución conocida

    La ecuación diferencial ordinaria del ejercicio tiene una solución conocida, lo que permite encontrar el error real en cada punto respecto a la aproximación estimada.

    y = e^x + x + x^2

    Note que el error crece al distanciarse del punto inicial

    Para las siguientes instrucciones, comente la última línea #plt.show() antes de continuar con:

    # ERROR vs solución conocida
    y_sol = lambda x: ((np.e)**x) + x + x**2
    
    yi_psol = y_sol(xi)
    errores = yi_psol - yi
    errormax = np.max(np.abs(errores))
    
    # SALIDA
    print('Error máximo estimado: ',errormax)
    print('entre puntos: ')
    print(errores)
    
    # GRAFICA [a,b+2*h]
    a = x0
    b = h*muestras+2*h
    muestreo = 10*muestras+2
    xis = np.linspace(a,b,muestreo)
    yis = y_sol(xis)
    
    plt.plot(xis,yis,linestyle='dashed', label='y solución conocida')
    plt.legend()
    plt.show()
    

    Se puede observar los siguientes resultados:

    Error máximo estimado:  0.0012745047595
    entre puntos: 
    [ 0.  0.000170  0.000377  0.000626  0.000922  0.00127 ]
    

    [ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]


    Differential equations, a tourist's guide | DE1. 3Blue1Brown. 31-mayo-2019
    Los Fundamentos de Ecuaciones Diferenciales que Nadie te Explica. 3Blue1Brown Español. 19-Junio-2024

     

  • 5.8 Integral y Derivada de f(x), expresiones con Sympy

    [ Derivadas f(x) ] [ Derivadas Sin Evaluar ] [ Integral definido [a,b] ] [ Integral Indefinida ]


    1. Funciones de prueba

    Para los ejemplos se usan f(x) de variable independiente 'x', y constantes 'a' y 'b'

    f(x) = a \cos(x) f(x) =a e^{-3x}

    [ Derivadas f(x) ] [ Derivadas Sin Evaluar ] [ Integral definido [a,b] ] [ Integral Indefinida ]
    ..


    2. Derivadas de f(x) con Sympy

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

    f(x): a*exp(-3*x)
    f(x) con sym.pprint:
       -3⋅x
    a⋅ℯ    
    
     df/dx: -3*a*exp(-3*x)
    df/dx con sym.pprint:
          -3⋅x
    -3⋅a⋅ℯ    
    
     d2f/dx2: 9*a*exp(-3*x)
    d2f/dx2 con sym.pprint:
         -3⋅x
    9⋅a⋅ℯ    
    

    Instrucciones en Python

    # derivadas de f(x) expresiones Sympy
    import sympy as sym
    # INGRESO
    a = sym.Symbol('a') # constantes sin valor
    b = sym.Symbol('b')
    x = sym.Symbol('x',real=True) # variable independente
    
    #fx = a*sym.cos(x)
    fx = a*sym.exp(-3*x)
    
    #PROCEDIMIENTO
    dfx = fx.diff(x,1)
    d2fx = fx.diff(x,2)
    
    # SALIDA
    print('f(x):',fx)
    print('f(x) con sym.pprint:')
    sym.pprint(fx)
    
    print('\n df/dx:',dfx)
    print('df/dx con sym.pprint:')
    sym.pprint(dfx)
    
    print('\n d2f/dx2:',d2fx)
    print('d2f/dx2 con sym.pprint:')
    sym.pprint(d2fx)
    

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

    Ejemplos:

    Polinomio de Taylor – Ejemplos con Sympy-Python

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

    [ Derivadas f(x) ] [ Derivadas Sin Evaluar ] [ Integral definido [a,b] ] [ Integral Indefinida ]
    ..


    2. Derivadas sin evaluar df(x)/dx con Sympy

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

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

    Para más adelante definir f(x).

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

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

    Ejemplos:

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

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

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

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

    [ Derivadas f(x) ] [ Derivadas Sin Evaluar ] [ Integral definido [a,b] ] [ Integral Indefinida ]
    ..


    3. Integrales definida de f(x) con Sympy

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

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

    [ Derivadas f(x) ] [ Derivadas Sin Evaluar ] [ Integral definido [a,b] ] [ Integral Indefinida ]
    ..


    4. Integrales Indefinida de f(x) con Sympy

    La operación se puede realizar con sym.integrate(fx,x), la expresión obtenida no añade la constante.

    # integral f(x) eindefinido xpresiones Sympy
    import sympy as sym
    # INGRESO
    a = sym.Symbol('a') # constantes sin valor
    b = sym.Symbol('b')
    c = sym.Symbol('c')
    x = sym.Symbol('x',real=True) # variable independente
    
    #fx = a*sym.cos(x)
    fx = a*sym.exp(-3*x)
    
    #PROCEDIMIENTO
    y = sym.integrate(fx,x) + c
    
    # SALIDA
    print('f(x):',fx)
    print('f(x) con sym.pprint:')
    sym.pprint(fx)
    
    print('\n y:',y)
    print('y con sym.pprint:')
    sym.pprint(y)
    

    con el siguiente resultado:

    f(x): a*exp(-3*x)
    f(x) con sym.pprint:
       -3⋅x
    a⋅ℯ    
    
     y: -a*exp(-3*x)/3 + c
    y con sym.pprint:
         -3⋅x    
      a⋅ℯ        
    - ─────── + c
         3       
    
    

    Referencia: https://docs.sympy.org/latest/modules/integrals/integrals.html


    [ Derivadas f(x) ] [ Derivadas Sin Evaluar ] [ Integral definido [a,b] ] [ Integral Indefinida ]

  • 5.7.2 Diferenciación numérica - Tablas con diferencias divididas

    Diferencias Divididas [ hacia adelante ] [ centradas ] [hacia atrás] ..


    Referencia: Chapra Fig.23.1 p669, Burden 4.1 p167, Rodríguez 8.2,3,4,6 p324

    Diferencias divididas hacia adelante

    Primera derivada

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

    Segunda derivada

    f''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+f(x_i)}{h^2} + O(h) f''(x_i) = \frac{-f(x_{i+3})+4f(x_{i+2})-5f(x_{i+1})+2f(x_i)}{h^2} + O(h^2)

    Tercera derivada

    f'''(x_i) = \frac{f(x_{i+3})-3f(x_{i+2})+3f(x_{i+1})-f(x_i)}{h^3} + O(h) f'''(x_i) = \frac{-3f(x_{i+4})+14f(x_{i+3})-24f(x_{i+2})+18f(x_{i+1})-5f(x_i)}{2h^3} + O(h^2)

    Cuarta derivada

    f''''(x_i) = \frac{f(x_{i+4})-4f(x_{i+3})+6f(x_{i+2})-4f(x_{i+1})+f(x_i)}{h^3} + O(h)

    Diferencias Divididas [ hacia adelante ] [ centradas ] [hacia atrás]

    ..


    Diferencias divididas centradas

    Primera derivada

    f'(x_i) = \frac{f(x_{i+1})-f(x_{i-1})}{2h} + O(h^2) f'(x_i) = \frac{-f(x_{i+2})+8f(x_{i+1})-8f(x_{i-1}) +f(x_{i-2})}{12h} + O(h^4)

    Segunda derivada

    f''(x_i) = \frac{f(x_{i+1})-2f(x_{i})+f(x_{i-1})}{h^2} + O(h^2) f''(x_i) = \frac{-f(x_{i+2})+16f(x_{i+1})-30f(x_{i})+16f(x_{i-1})-f(x_{i-2})}{12h^2} + O(h^4)

    Tercera derivada

    f'''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+2f(x_{i-1})-f(x_{i-2})}{2h^3} + O(h^2) f'''(x_i) = \frac{-f(x_{i+3})+8f(x_{i+2})-13f(x_{i+1})+13f(x_{i-1})-8f(x_{i-2})+f(x_{i-3})}{8h^3} + O(h^4)

    Diferencias Divididas [ hacia adelante ] [ centradas ] [hacia atrás]

    ..


    Diferencias divididas hacia atrás

    Primera derivada

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

    Segunda derivada

    f''(x_i) = \frac{f(x_{i})-2f(x_{i-1})+f(x_{i-2})}{h^2} + O(h) f''(x_i) = \frac{2f(x_{i})-5f(x_{i-1})+4f(x_{i-2})-f(x_{i-3})}{h^2} + O(h^2)

    Tercera derivada

    f'''(x_i) = \frac{f(x_{i})-3f(x_{i-1})+3f(x_{i-2})-f(x_{i-3})}{h^3} + O(h) f'''(x_i) = \frac{5f(x_{i})-18f(x_{i-1})+24f(x_{i-2})-14f(x_{i-3})+3f(x_{i-4})}{2h^3} + O(h^2)

    Diferencias Divididas [ hacia adelante ] [ centradas ] [hacia atrás]

  • 5.7.1 Diferenciación numérica - Error con df/dx en intervalo

    [ compara ] [ resultados ] [ algoritmo ]
    ..


    1. Compara derivadas numéricas con analíticas

    Para observar el efecto de disminuir h al aumentar los tramos, se realiza la gráfica de la derivada df/dx analítica versus las derivadas numéricas hacia atrás, hacia centrada y hacia adelante.

    Diferenciacion numérica en un intervalo compara con df/dxSe observa en la animación la reducción del error en cada diferencia dividida.

    [ compara ] [ resultados ] [ algoritmo ]
    ..


    2. Resultados del algoritmo

       tramos: 4 	, h: 0.5
    |errado_atras|:    0.3981193548993356
    |errado_centro|:   0.04939087954655119
    |errado_adelante|: 0.41231110309099034
    |errado_max|:      0.41231110309099034
    [  dfx,        df1_atras, df1_centro, df1_adelante]
    [[ 0.9610378          nan         nan  0.76041177]
     [ 0.49386065  0.76041177  0.44446977  0.12852777]
     [-0.26703531  0.12852777 -0.27540932 -0.67934641]
     [-1.07746577 -0.67934641 -1.04151373 -1.40368104]
     [-1.67397947 -1.40368104         nan         nan]]
       tramos: 8 	, h: 0.25
    |errado_atras|:    0.20757887012431775
    |errado_centro|:   0.016528173718511008
    |errado_adelante|: 0.20828851979214785
    |errado_max|:      0.20828851979214785
       tramos: 32 	, h: 0.0625
    |errado_atras|:    0.05218398011852454
    |errado_centro|:   0.0011877802167017393
    |errado_adelante|: 0.052183627077393824
    |errado_max|:      0.05218398011852454
       tramos: 64 	, h: 0.03125
    |errado_atras|:    0.026094601614190305
    |errado_centro|:   0.000302562296112141
    |errado_adelante|: 0.026094780174240162
    |errado_max|:      094780174240162

    [ compara ] [ resultados ] [ algoritmo ]
    ..


    3. Algoritmo con Python

    Las instrucciones en Python usadas para la gráfica son:

    # Diferenciacion Numerica - 1ra Derivada
    # Compara resultados en un intervalo
    import numpy as np
    
    # INGRESO
    fx = lambda x: np.sqrt(x)*np.sin(x)
    dfx = lambda x: np.sin(x)/(2*np.sqrt(x))+ np.sqrt(x)*np.cos(x)
    a = 1  # intervalo de integracion
    b = 3
    tramos = 64 # >=2, al menos 2 tramos
    
    # PROCEDIMIENTO
    muestras = tramos + 1
    h = (b-a)/tramos
    # puntos de muestras
    xi = np.linspace(a,b,muestras)
    fi = fx(xi)
    
    # tabla para comparar derivada y diferencias divididas
    tabla = np.zeros(shape=(muestras,4),dtype=float)
    for i in range(0,muestras,1): # muestras en intervalo
        df1_atras = np.nan ; df1_centro = np.nan; df1_adelante = np.nan
        
        dfxi = dfx(xi[i]) # derivada analitica, como referencia
        if i>0 and i<muestras: # hacia Atras, diferencias divididas
            df1_atras = (fi[i]-fi[i-1])/h
        if i>0 and i<(muestras-1): # Centro, diferencias divididas
            df1_centro = (fi[i+1]-fi[i-1])/(2*h)
        if i>=0 and i<(muestras -1): # hacia Adelante, diferencias divididas
            df1_adelante = (fi[i+1]-fi[i])/h
        tabla[i]=[dfxi,df1_atras,df1_centro,df1_adelante]
    
    # errado desde dfxi
    errado_atras = np.max(np.abs(tabla[1:,0]-tabla[1:,1]))
    errado_centro = np.max(np.abs(tabla[1:muestras-1,0]-tabla[1:muestras-1,2]))
    errado_adelante = np.max(np.abs(tabla[:muestras-1,0]-tabla[:muestras-1,3]))
    errado_max = np.max([errado_atras,errado_centro,errado_adelante])
    
    # SALIDA
    print('   tramos:', tramos,'\t, h:',h,)
    print('|errado_atras|:   ',errado_atras)
    print('|errado_centro|:  ',errado_centro)
    print('|errado_adelante|:',errado_adelante)
    print('|errado_max|:     ',errado_max)
    print('[  dfx,        df1_atras, df1_centro, df1_adelante]')
    print(tabla)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    
    # fx suave aumentando muestras
    muestrasfxSuave = tramos*10 + 1
    xk = np.linspace(a,b,muestrasfxSuave)
    dfk = dfx(xk)
    
    # Graficar f(x), puntos
    plt.plot(xk,dfk, label ='df(x)',linestyle='dotted')
    if tramos<64:
        plt.plot(xi,tabla[:,0], 'o') # muestras
    
    plt.plot(xi,tabla[:,1],label='df1_atras')
    plt.plot(xi,tabla[:,2],label='df1_centrada')
    plt.plot(xi,tabla[:,3],label='df1_adelante')
    
    plt.xlabel('xi')
    plt.ylabel('df(xi)')
    txt = 'Diferenciaci n Num rica'
    txt = txt + ', tramos:'+str(tramos)
    txt = txt + ', h:'+str(h)
    txt = txt + ', |errado_max|'+str(round(errado_max,4))
    plt.title(txt)
    plt.legend()
    plt.tight_layout()
    
    plt.show()
    

    [ compara ] [ resultados ] [ algoritmo ]

  • 5.7 Diferenciación numérica

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]

    Referencia: Chapra 23.1 p668, Rodríguez 8.2 p324

    La diferenciación numérica aproxima la derivada f'(x) utilizando muestras xi alrededor del punto de interés que se encuentran a distancia h.

    Diferenciacion numérica con varias h

    Se pueden generar fórmulas por diferencias divididas de alta exactitud
    tomando términos adicionales en la expansión de la serie de Taylor. Como referencia, el polinomio de Taylor muestra una aproximación de una función f(x) em el punto x0:

    P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 + ...

    de donde se obtiene la aproximación de las derivadas, al seleccionar la cantidad de términos y despejar, como se mostrará a continuación.

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]

    ..


    1. Primera Derivada con Diferencias divididas hacia adelante

    Una aproximación a primera derivada, usa los primeros dos términos del polinomio de Taylor alrededor de xi en para un punto a la derecha xi+1 a una distancia h = xi+1-xi

    Diferenciacion Adelante 01

    f(x_{i+1}) = f(x_i)+\frac{f'(x_i)}{1!} (h) + \frac{f''(x_i)}{2!}(h)^2 + ...

    se puede simplificar en un polinomio de grado uno y un término de error:

    f_{i+1} = f_i + (h)f'_i + O(h^2) ...

    Despejando la expresión para f'i


    f'_i = \frac{f_{i+1}-f_i}{h} = \frac{\Delta f_i}{h}

    La expresión también es la primera diferencia finita dividida con un error del orden O(h). (tema usado en interpolación).

    Revise que el término de la expresión queda O(h2)/h con lo que se disminuye el exponente en uno.

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]

    ..


    2. Primera derivada con diferencias divididas centradas

    Se realiza el mismo procedimiento que el anterior, usando un punto xi+1 y xi-1 alrededor de xi. En el término xi-1 el valor de h es negativo al invertir el orden de la resta.

    Diferenciacion Centrada 01

    f_{i+1} = f_i+\frac{f'_i}{1!}(h) + \frac{f''_i}{2!}(h)^2 + O(h^3) ... f_{i-1} = f_i-\frac{f'_i}{1!}(h) + \frac{f''_i}{2!}(h)^2

    restando la ecuaciones se tiene que

    f_{i+1} - f_{i-1} = (h)f'_i +(h)f'_i f_{i+1} - f_{i-1} = 2h f'_i

    La expresión de primera derivada usando un punto antes y otro después del punto central queda como:

    f'_i = \frac{f_{i+1} - f_{i-1}}{2h}

    con un error del orden O(h2)

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]
    ..


    3. Algoritmo en Python

    Compara las diferencias divididas hacia adelante, centro y atrás.

    El punto xk para la derivada se encuentra en la posición i , de dónde se toman para el cálculo los puntos xi  hacia la izquierda o derecha .

    Se requieren al menos 2 tramos para usar las fórmulas, para disponer de un valor de h dentro del intervalo.

    Para observar el error en la gráfica, se incorpora f'(x) obtenida de forma analítica como referencia en el ejercicio.

    El resultado del algoritmo es:

    diferenciación numérica:
    en xk: 1.5 	, h: 0.5
       tramos: 4
    [  dfx/dx,  valor,  error]
    ['df_analítica', 0.4938606479864909, 0]
    ['df_atras', 0.7604117685484817, 0.2665511205619908]
    ['df_centro', 0.4444697684399397, -0.04939087954655119]
    ['df_adelante', 0.12852776833139767, -0.3653328796550932]

    Instrucciones en Python

    # Diferenciación Numérica - 1ra Derivada
    # Compara con diferenciación numérica
    import numpy as np
    
    # INGRESO
    fx = lambda x: np.sqrt(x)*np.sin(x)
    dfx = lambda x: np.sin(x)/(2*np.sqrt(x))+ np.sqrt(x)*np.cos(x)
    xk = 1.5 # punto derivada
    
    a = 1  # intervalo de observación
    b = 3
    tramos = 64 # >=2, al menos 2 tramos 
    
    # PROCEDIMIENTO
    # revisa valores, Tarea: incluir en algoritmo
    cond1 = a<b
    cond2 = a<xk and xk<b
    cond3 = tramos>=2
    if not(cond1):
        print('a>b, no es ascendente')
    if not(cond2):
        print('ak fuera de intervalo [a,b]')
    if not(cond3):
        print('tramos deben ser al menos 2')
        
    # tramo menor hacia atras y adelante
    dx_atras = abs(xk-a) # xk hacia atras intervalo
    dx_adelante = abs(xk-b) # xk hacia frente intervalo
    dx_min = min(dx_atras,dx_adelante)
    
    h_tramo = (b-a)/tramos
    h = min(h_tramo,dx_min)
    
    xi_atras = np.sort(np.arange(xk,a-h/2,-h))
    xi_adelante = np.arange(xk,b+h/2,h)
    xi = np.concatenate((xi_atras,xi_adelante[1:]),axis=0)
    
    i = len(xi_atras)-1
    # puntos de muestras
    fi = fx(xi)
    tabla = [] # tangentes a xi[i]
    
    # derivada analítica, como referencia
    dfi = dfx(xi[i])
    px = lambda x: fi[i] + dfi*(x-xi[i])
    pxi = px(xi)
    tabla.append(['df_analítica',dfi,0,pxi,px])
    
    # hacia Atras, diferencias divididas
    df1_atras = (fi[i]-fi[i-1])/h
    px = lambda x: fi[i] + df1_atras*(x-xi[i])
    pxi = px(xi)
    errado = df1_atras - dfi
    tabla.append(['df_atras',df1_atras,errado,pxi,px])
    
    # Centro, diferencias divididas
    df1_centro = (fi[i+1]-fi[i-1])/(2*h)
    px = lambda x: fi[i] + df1_centro*(x-xi[i])
    pxi = px(xi)
    errado = df1_centro - dfi
    tabla.append(['df_centro',df1_centro,errado,pxi,px])
    
    # hacia Adelante, diferencias divididas
    df1_adelante = (fi[i+1]-fi[i])/h
    px = lambda x: fi[i] + df1_adelante*(x-xi[i])
    pxi = px(xi)
    errado = df1_adelante - dfi
    tabla.append(['df_adelante',df1_adelante,errado,pxi,px])
    
    # SALIDA
    print('diferenciación numérica:')
    print('en xk:',xk,'\t, h:',h,)
    print('   tramos:', tramos)
    print('[  dfx/dx,  valor,  error]')
    for unadifdiv in tabla:
        print(unadifdiv[0:3])

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]

    ..


    3.1Gráficas en Python

    diferencias divididas para 4 tramosSe complementa el algoritmo anterior con las siguientes instrucciones.

    # GRAFICA  --------------
    import matplotlib.pyplot as plt
    txt = 'Diferenciación Numérica'
    txt = txt + ', xi:'+str(xk)
    txt = txt + ', tramos:'+str(tramos)
    titulo = txt + ', h:'+str(h)
    
    # fx suave aumentando muestras
    muestrasfxSuave = tramos*10 + 1
    xk = np.linspace(a,b,muestrasfxSuave)
    fk = fx(xk)
    
    # Graficar f(x), puntos
    plt.ylim(min(fk),max(fk)*1.5)
    plt.plot(xk,fk, label ='f(x)',linestyle='dotted')
    if tramos<64:
        plt.plot(xi,fi, 'o') # muestras
    
    for unadifdiv in tabla: # tangentes a xi[i]
        linea = 'dashed'
        if unadifdiv[0]=='df_analítica':
            linea = 'solid'
        plt.plot(xi,unadifdiv[3],label=unadifdiv[0],
                 linestyle=linea)
    
    plt.plot(xi[i],fi[i],'ro') # punto i observado
    
    plt.xlabel('xi')
    plt.ylabel('f(xi)')
    plt.title(titulo)
    plt.legend()
    plt.tight_layout()
    plt.show()

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]


    Segundas derivadas

    Al continuar con el procedimiento mostrado se pueden obtener las fórmulas para segundas derivadas, las que se resumen en las tablas de Fórmulas de diferenciación por diferencias divididas.

    Diferencias Divididas [ hacia adelante ] [ centradas ] || [ algoritmo ] [ gráfica ]

  • 5.6 Integrales dobles sobre área rectangular con Python

    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]
    ..


    1. Integrales Dobles sobre área rectangular

    Referencia: Chapra 21.5 p643, Burden 4.9 p174

    Un ejemplo sencillo de integral doble es una función sobre un área rectangular. Las técnicas revisadas en las secciones anteriores de pueden reutilizar para la aproximación de integrales múltiples. Integra Multiple wireframe

    Considere la integral doble

    \int_R \int f(x,y) \delta A

    siendo

    R = {(x,y) | 
         x0 ≤ x ≤ xn, 
         y0 ≤ y ≤ yn }

    que es una región rectangular en el plano.

    \int_{x_o}^{x_n}\int_{y_o}^{y_n} f(x,y) \delta y \delta x \int_{y_o}^{y_n} \int_{x_o}^{x_n} f(x,y) \delta x \delta y

    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]

    ..


    2. Ejercicio

    Referencia3Eva_IIT2018_T1 Integral doble con Cuadratura de Gauss

    Aproxime el resultado de la integral doble (en área rectangular)

    \displaystyle \int_{0}^{\pi/4} \displaystyle\int_{0}^{\pi/4} \Big( 2y \sin(x) + \cos ^2 (x) \Big) \delta y \delta x

    Considere los tramosx = 5 como tamaños de paso en eje x , tramosy = 4 como tamaños de paso en eje y.

    Integral Multiple wireframe bar 01
    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]
    ..


    3. Desarrollo Analítico

    Use las fórmulas compuestas de Simpson, haciendo primero constante un valor en y, luego propagando el proceso. Al obtener los integrales en un sentido, repetir el proceso con los resultados en el otro eje.

    j: 0 , yj[0]: 0.0
    Fórmulas compuestas, tramos: 5
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 2 0 0]
     Simpson 3/8 : 0.4378989163640264
     Simpson 1/3 : 0.20482799857900125
    j: 1 , yj[1]: 0.19634954084936207
    Fórmulas compuestas, tramos: 5
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 2 0 0]
     Simpson 3/8 : 0.4807008818748735
     Simpson 1/3 : 0.2770455037573468
    j: 2 , yj[2]: 0.39269908169872414
    Fórmulas compuestas, tramos: 5
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 2 0 0]
     Simpson 3/8 : 0.5235028473857204
     Simpson 1/3 : 0.34926300893569256
    j: 3 , yj[3]: 0.5890486225480862
    Fórmulas compuestas, tramos: 5
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 2 0 0]
     Simpson 3/8 : 0.5663048128965673
     Simpson 1/3 : 0.42148051411403814
    j: 4 , yj[4]: 0.7853981633974483
    Fórmulas compuestas, tramos: 5
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 2 0 0]
     Simpson 3/8 : 0.6091067784074142
     Simpson 1/3 : 0.4936980192923838
    Fórmulas compuestas, tramos: 4
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 1 0]
     Simpson 3/8 : 0.4802254950852897
     trapecio : 0.20524320554554915
    xi [0.     0.1571 0.3142 0.4712 0.6283 0.7854]
    yj [0.     0.1963 0.3927 0.589  0.7854]
    f(x,y)
    [[1.     0.9755 0.9045 0.7939 0.6545 0.5   ]
     [1.     1.037  1.0259 0.9722 0.8853 0.7777]
     [1.     1.0984 1.1472 1.1505 1.1162 1.0554]
     [1.     1.1598 1.2686 1.3287 1.347  1.333 ]
     [1.     1.2213 1.3899 1.507  1.5778 1.6107]]
    Integral Volumen: 0.6854687006308389
    

    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]

    ..


    4. Algoritmo en Python

    El punto de partida es usar las Fórmulas compuestas con Δx segmentos desiguales, donde el algoritmo selecciona aplicar Simpson 3/8, Simpson 1/3 y trapecios según se cumplan los requisitos.

    Instrucciones en Python

    # 3Eva_IIT2018_T1 Integral doble con Cuadratura de Gauss
    # Integrales de volumen o dobles
    import numpy as np
    
    # INGRESO
    fxy = lambda x,y: 2*y*np.sin(x)+(np.cos(x))**2 
    x0 = 0    # intervalo x
    xn = np.pi/4
    y0 = 0    # intervalo y
    yn = np.pi/4
    tramosx = 5 # tramos por eje
    tramosy = 4
    titulo = 'f(x,y)'
    precision = 4 # decimales a mostrar
    
    # PROCEDIMIENTO
    xi = np.linspace(x0,xn,tramosx+1)
    yj = np.linspace(y0,yn,tramosy+1)
    n = len(xi)
    m = len(yj)
    hx = (xn-x0)/tramosx # tamaños de paso
    hy = (yn-y0)/tramosy
    Xi,Yj = np.meshgrid(xi,yj)
    
    # Matriz u[xi,yj] para f(x,y)
    F = np.zeros(shape=(n,m),dtype=float)
    F = fxy(Xi,Yj)
    
    def simpson_compuesto(xi,fi,vertabla=False,casicero=1e-7):
        '''Método compuesto Simpson 3/8, Simpson 1/3 y trapecio
        salida: integral,cuenta_h
        '''
        # vectores como arreglo, numeros reales
        xi = np.array(xi,dtype=float)
        fi = np.array(fi,dtype=float)
        n = len(xi)
    
        # Método compuesto Simpson 3/8, Simpson 1/3 y trapecio
        cuenta_h = (-1)*np.ones(n,dtype=int) # sin usar
        suma = 0 # integral total
        suma_parcial =[] # integral por cada método
        i = 0
        while i<(n-1): # i<tramos, al menos un tramo
            tramo_usado = False
            h = xi[i+1]-xi[i] # tamaño de paso, supone constante
            
            # tres tramos iguales
            if (tramo_usado==False) and (i+3)<n:
                d2x = np.diff(xi[i:i+4],2) # diferencias entre tramos
                errado = np.max(np.abs(d2x))
                if errado<casicero: # Simpson 3/8
                    S38 = (3/8)*h*(fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3])
                    suma = suma + S38
                    cuenta_h[i:i+3] = [3,0,0] # Simpson 3/8
                    suma_parcial.append(['Simpson 3/8',S38])
                    i = i+3
                    tramo_usado = True
    
            # dos tramos iguales
            if (tramo_usado==False) and (i+2)<n:
                d2x = np.diff(xi[i:i+3],2) # diferencias entre tramos
                errado = np.max(np.abs(d2x))
                if errado<casicero: # Simpson 1/3
                    S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
                    suma = suma + S13
                    cuenta_h[i:i+2] = [2,0]
                    suma_parcial.append(['Simpson 1/3',S13])
                    i = i+2
                    tramo_usado = True
                    
            # un tramo igual
            if (tramo_usado == False) and (i+1)<n:
                trapecio = (h/2)*(fi[i]+fi[i+1])
                suma = suma + trapecio
                cuenta_h[i:i+1] = [1] # usar trapecio
                suma_parcial.append(['trapecio',trapecio])
                i = i+1
                tramo_usado = True
                
        cuenta_h[n-1] = 0 # ultima casilla
        
        if vertabla==True: #mostrar datos parciales
            print('Fórmulas compuestas, tramos:',n-1)
            print('métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado')
            print('tramos iguales en xi:',cuenta_h)
            for unparcial in suma_parcial:
                print('',unparcial[0],':',unparcial[1])
        
        return(suma,cuenta_h)
    
    # PROGRAMA --------------
    areaj = []
    cuenta_hj = []
    j = 0 # y[0] constante
    for j in range(0,tramosy+1,1): # cada yj
        print('j:',j,', yj['+str(j)+']:',yj[j])
        xj = Xi[j]
        fj = F[j]
        area,cuenta_h=simpson_compuesto(xj,fj,
                                        vertabla=True,
                                        casicero=1e-7)
        areaj.append(area)
        cuenta_hj.append(cuenta_h)
    area,cuenta_h = simpson_compuesto(yj,areaj,
                                      vertabla=True,
                                      casicero=1e-7)
    # SALIDA
    np.set_printoptions(precision)
    print('xi',xi)
    print('yj',yj)
    print(titulo)
    print(F)
    print('Integral Volumen:', area)
    

    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]
    ..


    4. Gráfica de malla

    Permite observar la superficie de f(x,y) sobre la que se calculará el volumen respecto al plano x,y en cero.

    Integra Multiple wireframe

    Instrucciones adicionales al algoritmo

    # GRAFICA 3D de malla ------
    import matplotlib.pyplot as plt
    
    # Malla para cada eje X,Y
    fig3D = plt.figure()
    graf3D = fig3D.add_subplot(111, projection='3d')
    graf3D.plot_wireframe(Xi,Yj,F,label='f(x,y)')
    graf3D.plot(0,0,0,'ro',label='[0,0,0]')
    graf3D.set_title(titulo)
    graf3D.set_xlabel('x')
    graf3D.set_ylabel('y')
    graf3D.set_zlabel('z')
    graf3D.set_title(titulo)
    graf3D.legend()
    graf3D.view_init(35, -45)
    plt.tight_layout()
    #plt.show()
    

    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]
    ..


    5. Gráfica de barras

    Realizado para observar las barras con lados de segmentos dx=hx, dy=hy.

    Integral Multiple bar

    Instrucciones complementarias a la gráfica anterior.

    # Grafica barras 3D ------
    factor = 0.95 # barra proporcion dx, dy
    
    fig_3D = plt.figure()
    graf_3D = fig_3D.add_subplot(111,projection='3d')
    
    F_min = np.min([np.min(F),0])
    F_max = np.max([np.max(F)])
    color_F = F[:m-1,:n-1].ravel()/F_max
    color_ij = plt.cm.rainbow(color_F)
    
    graf_3D.bar3d(Xi[:m-1,:n-1].ravel(),
                  Yj[:m-1,:n-1].ravel(),
                  np.zeros((n-1)*(m-1)),
                  hx*factor,hy*factor,
                  F[:m-1,:n-1].ravel(),
                  color=color_ij)
    graf_3D.plot(0,0,0,'ro',label='[0,0,0]')
    
    # escala eje z
    graf_3D.set_zlim(F_min,F_max)
    # etiqueta
    graf_3D.set_xlabel('xi')
    graf_3D.set_ylabel('yj')
    graf_3D.set_zlabel('z')
    graf_3D.set_title(titulo)
    graf_3D.legend()
    
    # Barra de color
    color_escala = plt.Normalize(F_min,F_max)
    color_barra = plt.cm.ScalarMappable(norm=color_escala,
                    cmap=plt.colormaps["rainbow"])
    fig_3D.colorbar(color_barra,ax=graf_3D,label="F[xi,yj]")
    
    graf_3D.view_init(35, -45)
    plt.tight_layout()
    plt.show()
    

    [ Integral doble ] [ ejercicio ] [ algoritmo ] | gráfica: [ mallas ] [ barras ]

  • 5.5.1 Cuadratura con dos puntos - Experimento con Python

    [ Cuadratura de Gauss ] [ Experimento 2 puntos ] [ Experimento algoritmo ]
    ..


    Cuadratura con dos puntos - Experimento

    Para visualizar el concepto de Cuadratura de Gauss de dos puntos considere lo siguiente:

    Se tiene un corte transversal del un recipiente rectangular lleno de líquido limitado en x entre [-1,1], al que se le ha depositado encima otro recipiente con perfil f(x) = x2 hasta que reposan sus extremos en x[-1,1].

    La altura de ambos recipientes es la misma.

    Cuadratura 2p 01

    La superficie entre f(x) y el eje x es el integral de f(x) en el intervalo.

    Si suponemos que la figura es el corte transversal de una vasija y la parte en amarillo es líquida, la vasija ha desplazado el líquido que ocupa ahora el "área" mostrada en la gráfica que corresponde al integral de f(x)=x2. entre [-1,1].

    Ahora, suponga que se perfora el perfil de f(x) en dos puntos equidistantes cercanos  x=0. Los orificios permitirían el desplazamiento del liquido al interior de f(x) que dejando pasar suficiente tiempo, permitiría tener todo el líquido en el recipiente rectangular entre [-1,1] como una línea horizontal.

    Cuadratura 2p 02

    Podría medir la altura que tiene el líquido y que tiene un equivalente en un punto f(x1). Debería encontrar el valor de x1 que permite disponer del mismo valor entre el área debajo de f(x) y el rectángulo del corte transversal amarillo ahora formado.

    Se usa el resultado analítico del integral restando el área del rectángulo obtenido al evaluar la función f(x) entre [0,1], teniendo un problema de búsqueda de raíces. Obtenemos el valor de x1.

    Se muestra que el área bajo f(x) es equivalente al área del rectángulo conformado.

    Si  utilizamos el desplazamiento horizontal desde el centro para un punto encontrado como un "factor", tendremos que el área del rectángulo se mantendría equivalente, y el desplazamiento proporcional a la mitad del intervalo si se varía el intervalo de observación. Este factor coincide con el factor de Cuadratura de Gauss de dos puntos.

    funcion  fx:   x**2
    Integral Fx:   x**3/3
    I analitico:   0.6666666666666666
    I aproximado:  0.6666666666654123
    desplaza centro:   0.5773502691890826
    factor desplaza:   0.5773502691890826
    Factor CuadGauss:  0.5773502691896258
    erradoFactor:    1.2545520178264269e-12
    error integral:  1.2543299732215019e-12

    El error del integral es del orden de 10-12


    Cambiamos la figura geométrica a un trapecio generado por la recta que pasa por los puntos xi desplazados desde el centro.

    Usamos la función f(x) = x2 + x + 1, observaremos si los resultado son equivalentes.

    La figura al inicio del experimento será:

    Cuadratura 2p 03

    Luego de realizar realizar el mismo cálculo anterior usando un equivalente a trapecio se tiene:

    Cuadratura 2p 04

    con valores numéricos:

    funcion  fx:   x**2 + x + 1
    Integral Fx:   x**3/3 + x**2/2 + x
    I analitico:   2.6666666666666665
    I aproximado:  2.6666666666654124
    desplaza centro:   0.5773502691890826
    factor desplaza:   0.5773502691890826
    Factor CuadGauss:  0.5773502691896258
    erradoFactor:    1.2545520178264269e-12
    error integral:  1.2541079286165768e-12

    El error del integral es también del orden de 10-12, además observe que el factor de cuadratura de Gauss se mantiene.

    [ Cuadratura de Gauss ] [ Experimento 2 puntos ] [ Experimento algoritmo ]


    Tarea

    Realice el experimento usando un polinomio de grado superior y observe los errores para el integral y las diferencia con el coeficiente de Cuadratura de 2 puntos.

    [ Cuadratura de Gauss ] [ Experimento 2 puntos ] [ Experimento algoritmo ]

    ..


    Algoritmo con Python

    Para resumir la cantidad de instrucciones, se usa el método de la bisección desde la librería Scipy y el subgrupo de funciones de optimización.

    Los cálculos para realizar las gráficas se tratan en un bloque luego de mostrar los resultado principales.

    # Integración: Cuadratura de Gauss de dos puntos
    # modelo con varios tramos entre [a,b]
    # para un solo segmento.
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    import scipy.optimize as op
    
    # INGRESO
    x = sym.Symbol('x')
    
    # fx = (x)**2
    fx = x**2 + x + 1
    # fx = 0.2 + 25.0*x-200*(x**2)+675.0*(x**3)-900.0*(x**4)+400.0*(x**5)
    
    a = -1
    b = 1
    muestras = 51
    tolera = 1e-12
    iteramax = 100
    
    # PROCEDIMIENTO
    # Desarrollo analítico con Sympy
    Fx = sym.integrate(fx,x)
    Fxn = sym.lambdify('x',Fx,'numpy')
    Fiab = Fxn(b)-Fxn(a)
    
    # Busca igualar trapecio con Integral analitico
    fxn = sym.lambdify('x',fx,'numpy')
    base = b-a
    mitad = base/2
    xc = (a+b)/2  # centro
    diferencia = lambda x: Fiab-base*(fxn(xc-x)+fxn(xc+x))/2
    desplazado =  op.bisect(diferencia,0,mitad,
                            xtol=tolera,maxiter=iteramax)
    factor = desplazado/mitad
    
    # Integral aproximando con trapecio
    x0 = xc - factor*mitad
    x1 = xc + factor*mitad
    Faprox = base*(fxn(x0)+fxn(x1))/2
    
    # Integral cuadratura Gauss
    xa = xc + mitad/np.sqrt(3)
    xb = xc - mitad/np.sqrt(3)
    FcuadG = base*(fxn(xa)+fxn(xb))/2
    erradofactor = np.abs(FcuadG - Faprox)
    erradoIntegral = np.abs(Fiab-Faprox)
    # SALIDA
    print('funcion  fx:  ', fx)
    print('Integral Fx:  ', Fx)
    print('I analitico:  ', Fiab)
    print('I aproximado: ', Faprox)
    print('desplaza centro:  ', desplazado)
    print('factor desplaza:  ', factor)
    print('Factor CuadGauss: ', 1/np.sqrt(3))
    print('erradoFactor:   ', erradofactor)
    print('error integral: ', erradoIntegral) 
    
    # Grafica
    # Para GRAFICAR 
    # Para gráfica f(x)
    xi = np.linspace(a,b,muestras)
    fi = fxn(xi)
    
    # Para gráfica Trapecio
    m = (fxn(x1)-fxn(x0))/(x1-x0)
    trapeciof = lambda x: fxn(x0)+m*(x-x0)
    trapecioi = trapeciof(xi)
    
    # Areas Trapecio para cada punto que busca
    k = int(muestras/2)
    xicg = xi[k:muestras-1]
    Fcg = [base*(fxn(xi[k+0])+fxn(xi[k-0]))/2]
    for i in range(1,k,1):
        untrapecio = base*(fxn(xi[k+i])+fxn(xi[k-i]))/2
        Fcg.append(untrapecio)
    
    # Punto buscado
    Fiaprox = base*(fxn(x1)+fxn(x0))/2
    
    Fi = Fxn(xi)-Fxn(a)
    
    # Areas de curvas y trapecio
    
    plt.subplot(211) # Grafica superior
    plt.xlim(a,b)
    plt.plot(xi, fi, label='f(x)')
    # Solo fi
    # plt.fill_between(xi,0, fi,
    #                label='integral fi',
    #                 color='yellow')
    # usando cuadratura
    plt.fill_between(xi,0, trapecioi,
                     label='Cuadratura 2 puntos',
                     color='yellow')
    plt.axvline(x0,color='white')
    plt.axvline(x1,color='white')
    plt.plot([x0,x1],[fxn(x0),fxn(x1)],
             'ro', label='x0,x1')
    plt.axvline(0,color='black')
    plt.xlabel('x')
    plt.ylabel('f(x) y Cuadratura de 2 puntos')
    plt.legend()
    
    # Valores de integrales
    plt.subplot(212) # Grafica inferior
    plt.xlim(a,b)
    plt.axhline(Fiab, label='F[a,b]')
    # plt.plot(xi,Fi,label='F(x)')
    plt.plot(xicg,Fcg,color='orange',label='Aprox Trapecio')
    
    plt.axvline(x1,color='yellow')
    plt.axvline((1/np.sqrt(3))*(b-a)/2 + xc ,color='magenta')
    plt.plot(x1,Fiaprox,'ro', label='x0,x1')
    plt.axvline(0,color='black')
    
    plt.xlabel('x')
    plt.legend()
    plt.ylabel('Integrando')
    plt.show()
    

    [ Cuadratura de Gauss ] [ Experimento 2 puntos ] [ Experimento algoritmo ]

  • 5.5 Cuadratura de Gauss con Python

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]
    ..


    1. Cuadratura de Gauss

    Referencia: Chapra 22.3 p655, Rodríguez 7.3 p294, Burden 4.7 p168

    La cuadratura de Gauss aproxima el integral de una función en un intervalo [a,b] centrado en cero mediante un cálculo numérico con menos operaciones y evaluaciones de la función. Se representa como una suma ponderada:

    I \cong \Big( \frac{b-a}{2} \Big)\Big(c_0f(x_a) + c_1f(x_b)\Big)

    regla Cuadratura Gauss 01

    para la fórmula de dos puntos con referencia a una función centrada en cero y ancho unitario a cada lado, se tiene:

    c_0 = c_1 = 1, x_0 = -\frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}

    Para un intervalo de evaluación desplazado en el eje x se requiere convertir los puntos al nuevo intervalo. Se desplaza el punto cero al centro del intervalo [a,b] y se corrige el desplazamiento hacia la izquierda y derecha del centro con x0 y x1.

    Cuadratura Gauss integración xa xb animado

    x_a = \frac{b+a}{2} - \frac{b-a}{2}\Big(\frac{1}{\sqrt{3}} \Big) x_b = \frac{b+a}{2} + \frac{b-a}{2}\Big(\frac{1}{\sqrt{3}} \Big)

    con lo que el resultado aproximado del integral se convierte en:

    I \cong \frac{b-a}{2}(f(x_a) + f(x_b))

    cuya fórmula es semejante a una mejor aproximación de un trapecio, cuyos promedios de alturas son puntos internos de [a,b], concepto mostrado en la gráfica.

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]
    ..


    2. Ejercicio

    Para el ejercicio y comparación de resultado con los otros métodos, se realiza el cálculo para un tramo en el intervalo [a,b].

    \int_1^3 \sqrt{x} \sin(x) dx x_a = \frac{3+1}{2} - \frac{3-1}{2}\Big(\frac{1}{\sqrt{3}}\Big)= 1.4226 x_b = \frac{3+1}{2} + \frac{3-1}{2}\Big(\frac{1}{\sqrt{3}} \Big) = 2.5773 f(1.4226) = \sqrt{1.4226} \sin(1.4226) = 1.1796 f(2.5773) = \sqrt{2.5773} \sin(2.5773) = 0.8585

    con lo que el resultado aproximado del integral se convierte en:

    I \cong \frac{3-1}{2}(1.1796 + 0.8585) = 2.0382

    que usando instrucciones de Python para obtener los valores:

    >>> import numpy as np
    >>> (3+1)/2-(3-1)/2*(1/np.sqrt(3))
    1.4226497308103743
    >>> (3+1)/2+(3-1)/2*(1/np.sqrt(3))
    2.5773502691896257
    >>> fx = lambda x: np.sqrt(x)*np.sin(x)
    >>> fx(1.4226)
    1.1796544827404145
    >>> fx(2.5773)
    0.8585957175067221
    >>> ((3-1)/2)*(fx(1.4226)+fx(2.5773))
    2.0382
    

    el resultado se puede mejorar aumentando el número de tramos en el intervalo [a,b]. Por ejemplo, el resultado usando 4 tramos el resultado es semejante al usar el método del trapecio con 128 tramos, lo que muestra el ahorro en cálculos entre los métodos

    Integral:  2.05357719003
    >>> 
    

    Concepto:

    Ejercicio:

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]
    ..


    3. Algoritmo con Python

    Para el ejercicio anterior, usando al 4 segmentos y en cada uno aplicando Cuadratura de Gauss, el integral resultante es:

    Factores Gauss-Legendre puntos: 2
    xgl: [-0.5773502691896258, 0.5773502691896258]
    cgl: [1.0, 1.0]
    tabla por intervalos [a,b]
    [a,b] : [1.  1.5]
    xi : [1.10566 1.39434]
    fi : [0.93979 1.16248]
    area : 0.5255697290782242
    [a,b] : [1.5 2. ]
    xi : [1.60566 1.89434]
    fi : [1.26638 1.30494]
    area : 0.642828854689653
    [a,b] : [2.  2.5]
    xi : [2.10566 2.39434]
    fi : [1.24843 1.05163]
    area : 0.5750145898962924
    [a,b] : [2.5 3. ]
    xi : [2.60566 2.89434]
    fi : [0.82428 0.41638]
    area : 0.31016401636449004
    Integral:  2.0535771900286597
    

    Instrucciones en Python usando la Cuadratura de Gauss de dos puntos para una función f(x):

    # Integración: Cuadratura de Gauss de 2 puntos
    # modelo con varios tramos entre [a,b]
    import numpy as np
    
    # INGRESO
    fx = lambda x: np.sqrt(x)*np.sin(x)
    a = 1 # intervalo de integración
    b = 3
    tramos = 4  # subintervalos a integrar
    precision = 5 # decimales en tabla
    
    # PROCEDIMIENTO
    # cuadratura de 2 puntos
    n_puntos = 2
    xgl = [-1/np.sqrt(3),1/np.sqrt(3)]
    cgl = [1.,1.]
    # cuadratura de n_puntos, fórmulas Gauss-Legendre
    #xgl, cgl = np.polynomial.legendre.leggauss(2)
    
    x_h = np.linspace(a,b,tramos+1)
    tabla = {}
    suma = 0
    for k in range(0,tramos,1):
        a = x_h[k]
        b = x_h[k+1]
        centro = (a+b)/2
        mitad = (b-a)/2
        xa = centro + xgl[0]*mitad
        xb = centro + xgl[1]*mitad
    
        area = ((b-a)/2)*(cgl[0]*fx(xa) + cgl[1]*fx(xb))
        tabla[k]= {'[a,b]': np.array([a,b]),
                   'xi': np.array([xa,xb]),
                   'fi': np.array([fx(xa),fx(xb)]),
                   'area':area
                   }
        suma = suma + area
    
    # SALIDA
    np.set_printoptions(precision)
    print('Factores Gauss-Legendre puntos:',n_puntos)
    print('xgl:',xgl)
    print('cgl:',cgl)
    print('tabla por intervalos [a,b]')
    for k in range(0,tramos,1):
        for entrada in tabla[k]:
            print(entrada,':',tabla[k][entrada])
    print('Integral: ', suma)
    

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]
    ..


    4. Gráfica por tramos

    La gráfica que complementa el resultado anterior, se realiza añadiendo las instrucciones presentadas a continuación.

    Considere que la gráfica es útil con pocos tramos en el intervalo[a,b]

    Cuadratura Gauss 2 puntos, varios tramos

    # GRAFICA -------------
    # concepto con 'pocos' tramos/segmentos
    import matplotlib.pyplot as plt
    
    titulo = 'Cuadratura Gauss-Legendre puntos:'+str(n_puntos)
    
    a = np.min(x_h); b = np.max(x_h)
    muestras_k = tramos*10+1 # fx suave en grafico
    xk = np.linspace(a,b,muestras_k)
    fk = fx(xk)
    
    # tramos marcas de división
    for i in range(0,tramos+1,1):
        plt.axvline(x_h[i],linestyle='dashed',
                    color='tab:gray')
    for k in range(0,tramos,1):
        xi = tabla[k]['xi']
        fi = tabla[k]['fi']
        for i in range(0,n_puntos,1):
            plt.vlines(xi[i],0,fi[i],
                       linestyle='dotted',
                       color='gray')
            plt.plot(xi[i],fi[i],'o',color='gray')
    if n_puntos>2:
        plt.fill_between(xk,0,fk,color='tab:olive')
    if n_puntos==2: # recta trapecio
        for k in range(0,tramos,1):
            xi = tabla[k]['xi']
            fi = tabla[k]['fi']
            df = (fi[1]-fi[0])/(xi[1]-xi[0]) # pendiente
            f0 = fi[0] + df*(x_h[k]-xi[0])   # en xi[i]
            f1 = fi[0] + df*(x_h[k+1]-xi[0]) # en xi[i+1]
            plt.fill_between([x_h[k],x_h[k+1]],
                             0,[f0,f1],
                             color='tab:olive')
            plt.plot([x_h[k],x_h[k+1]],[f0,f1],
                     linestyle='dashed',color='orange')
    # linea fx
    plt.plot(xk,fk, label='f(x)', color = 'blue')
    plt.title(titulo)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.legend()
    plt.tight_layout()
    plt.show()
    

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]
    ..


    5. Tabla Cuadratura de Gauss-Legendre para n puntos

    Referencia: Chapra Tabla 22.1  p661, Burden 10Ed Tabla 4.12 p172

    Para un integral con intervalo centrado en el origen.

    I \cong c_0 f(x_0) + c_1f(x_1) + … + c_{n-1}f(x_{n-1})
    Factores usados en las fórmulas de Gauss-Legendre.
    Puntos Factor de ponderación Argumentos de la función Error de truncamiento
    2 c0 = 1.0000000
    c1 = 1.0000000
    x0 = - 0.577350269
    x1 = 0.577350269
    ≅ f(4)(x )
    3 c0 = 0.5555556
    c1 = 0.8888889
    c2 = 0.5555556
    x0 = - 0.774596669
    x1 = 0.0
    x2 = 0.774596669
    ≅ f(6)(x )
    4 c0 = 0.3478548
    c1 = 0.6521452
    c2 = 0.6521452
    c3 = 0.3478548
    x0 = - 0.861136312
    x1 = - 0.339981044
    x2 = 0.339981044
    x3 = 0.861136312
    ≅f (8)(x )

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]


    6. Función en librería np.polynomial.legendre.leggauss para n_puntos

    Referencianp.polynomial.legendre.leggauss

    Para obtener los valores de la tabla, se puede usar la librería Numpy.

    import numpy as np
    # n_puntos Gauss-Legendre
    n_puntos = 2
    xg, cg = np.polynomial.legendre.leggauss(n_puntos)
    print("xg",xg)
    print("cg",cg)
    

    que presenta como resultado:

    xg [-0.57735027  0.57735027]
    cg [1. 1.]

    donde se puede obtener los valores de la tabla con mas puntos.

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]

    ..


    7. Algoritmo de Gauss-Legendre n_puntos

    La gráfica se realiza con la sección gráfica anterior.

    # Integración: Cuadratura de Gauss de n_puntos
    # modelo con varios tramos entre [a,b]
    import numpy as np
    
    # INGRESO
    fx = lambda x: np.sqrt(x)*np.sin(x)
    a = 1 # intervalo de integración
    b = 3
    n_puntos = 5 # n_puntos>=2
    tramos = 3  # subintervalos a integrar
    precision = 5 # decimales en tabla
    
    # PROCEDIMIENTO
    # cuadratura de n_puntos, fórmulas Gauss-Legendre
    xgl, cgl = np.polynomial.legendre.leggauss(n_puntos)
    
    x_h = np.linspace(a,b,tramos+1) # tramos
    tabla = {}
    suma = 0
    for k in range(0,tramos,1):
        a = x_h[k] # intervalo[a,b]
        b = x_h[k+1]
        centro = (a+b)/2
        mitad = (b-a)/2
        xi=[]; fi=[]; suma_gl=0
        for j in range(0,n_puntos,1):
            xj = centro + xgl[j]*mitad
            fj = fx(xj)
            suma_gl = suma_gl + cgl[j]*fj
            xi.append(xj)
            fi.append(fj)
        area = ((b-a)/2)*suma_gl
        tabla[k]= {'[a,b]': np.array([a,b]),
                   'xi': np.array(xi),
                   'fi': np.array(fi),
                   'area':area
                   }
        suma = suma + area
    
    # SALIDA
    np.set_printoptions(precision)
    print('Factores Gauss-Legendre puntos:',n_puntos)
    print('xgl:',xgl)
    print('cgl:',cgl)
    print('tabla por intervalos [a,b]')
    for k in range(0,tramos,1):
        for entrada in tabla[k]:
            print(entrada,':',tabla[k][entrada])
    print('Integral: ', suma)
    

    resultado del algoritmo:

    Factores Gauss-Legendre puntos: 5
    xgl: [-0.90618 -0.53847  0.       0.53847  0.90618]
    cgl: [0.23693 0.47863 0.56889 0.47863 0.23693]
    tabla por intervalos [a,b]
    [a,b] : [1.      1.66667]
    xi : [1.03127 1.15384 1.33333 1.51282 1.63539]
    fi : [0.87127 0.98214 1.1223  1.2279  1.27616]
    area : 0.7350121364546989
    [a,b] : [1.66667 2.33333]
    xi : [1.69794 1.82051 2.      2.17949 2.30206]
    fi : [1.29253 1.30741 1.28594 1.21116 1.12934]
    area : 0.8369414195585647
    [a,b] : [2.33333 3.     ]
    xi : [2.36461 2.48718 2.66667 2.84616 2.96873]
    fi : [1.07815 0.95996 0.74672 0.4912  0.29637]
    area : 0.4816765248057005
    Integral:  2.053630080818964

    La gráfica se realiza con la sección gráfica anterior.

    Cuadratura Gauss Legendre 5 puntos tramos 3

    [ Cuadratura Gauss ] [ Ejercicio ] [ Algoritmo 2puntos ] [ gráfica ]
    [ tabla puntos ] [ algoritmo n_puntos]

  • 5.4 Fórmulas compuestas con Δx segmentos desiguales con Python

    [ Concepto ] [ Algoritmo /función ] [ Gráfica ] [ diferencias entre tramos ]

    ..


    Referencia: Chapra 21.3 p640

    En la práctica, existen muchas situaciones con segmentos o tramos de tamaños desiguales. Por ejemplo, los datos obtenidos en experimentos frecuentemente son de este tipo. Para integración numérica, se busca identificar los segmentos de igual tamaño y aplicar: Simpson 3/8, Simpson 1/3  o trapecio.

    Integral Fórrmula Compuesta Simpson 3/8, Simpson 1/3, Trapecio

    El algoritmo se desarrolla al incorporar cada método como integral numérico de "fórmulas compuestas".

    [ Concepto ] [ Algoritmo /función ] [ Gráfica ] [ diferencias entre tramos ]


    1. Ejercicio

    Los datos de ejemplo son:

    xi = [1.        , 1.22222222, 1.44444444, 1.66666667,
          1.88888889, 2.11111111, 2+1/3, 2+1/3+0.2,
          2+1/3+0.4, 3.        ]
    fi = [0.84147098, 1.03905509, 1.19226953, 1.28506615,
          1.30542157, 1.24598661, 1.10453193, 0.90952929,
          0.65637234, 0.24442702]

    [ Concepto ] [ Algoritmo /función ] [ Gráfica ] [ diferencias entre tramos ]
    ..


    2. Algoritmo como función en Python

    Los resultados del algoritmo muestran los detalles parciales al aplicar cada método acorde a los requisitos de cada uno en el siguiente orden: Tres tramos iguales permiten aplicar Simpson de 3/8, Dos tramos iguales permiten aplicar Simpson de 1/3, Un tramo desigual le aplica trapecio.

    Los métodos usados de identifican por el arreglo de tramos iguales:

    Fórmulas compuestas, tramos: 9
    métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado
    tramos iguales en xi: [3 0 0 3 0 0 2 0 1 0]
     Simpson 3/8 : 0.7350425751495739
     Simpson 3/8 : 0.8369852099634817
     Simpson 1/3 : 0.3599347620000003
     trapecio : 0.1201065813333333
    tramos iguales en xi: [3 0 0 3 0 0 2 0 1 0]
    Integral: 2.0520691284463894

    Instrucciones en Python

    # Integración Simpson 3/8, Simpson 1/3 y trapecio
    # tramos no iguales, formulas compuestas
    import numpy as np
    
    # INGRESO
    xi = [1.        , 1.22222222, 1.44444444, 1.66666667,
          1.88888889, 2.11111111, 2+1/3, 2+1/3+0.2,
          2+1/3+0.4, 3.        ]
    fi = [0.84147098, 1.03905509, 1.19226953, 1.28506615,
          1.30542157, 1.24598661, 1.10453193, 0.90952929,
          0.65637234, 0.24442702]
    
    # PROCEDIMIENTO
    def simpson_compuesto(xi,fi,vertabla=False,casicero=1e-7):
        '''Método compuesto Simpson 3/8, Simpson 1/3 y trapecio
        salida: integral,cuenta_h
        '''
        # vectores como arreglo, numeros reales
        xi = np.array(xi,dtype=float)
        fi = np.array(fi,dtype=float)
        n = len(xi)
    
        # Método compuesto Simpson 3/8, Simpson 1/3 y trapecio
        cuenta_h = (-1)*np.ones(n,dtype=int) # sin usar
        suma = 0 # integral total
        suma_parcial =[] # integral por cada método
        i = 0
        while i<(n-1): # i<tramos, al menos un tramo
            tramo_usado = False
            h = xi[i+1]-xi[i] # tamaño de paso, supone constante
            
            # tres tramos iguales
            if (tramo_usado==False) and (i+3)<n:
                d2x = np.diff(xi[i:i+4],2) # diferencias entre tramos
                errado = np.max(np.abs(d2x))
                if errado<casicero: # Simpson 3/8
                    S38 = (3/8)*h*(fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3])
                    suma = suma + S38
                    cuenta_h[i:i+3] = [3,0,0] # Simpson 3/8
                    suma_parcial.append(['Simpson 3/8',S38])
                    i = i+3
                    tramo_usado = True
    
            # dos tramos iguales
            if (tramo_usado==False) and (i+2)<n:
                d2x = np.diff(xi[i:i+3],2) # diferencias entre tramos
                errado = np.max(np.abs(d2x))
                if errado<casicero: # Simpson 1/3
                    S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
                    suma = suma + S13
                    cuenta_h[i:i+2] = [2,0]
                    suma_parcial.append(['Simpson 1/3',S13])
                    i = i+2
                    tramo_usado = True
                    
            # un tramo igual
            if (tramo_usado == False) and (i+1)<n:
                trapecio = (h/2)*(fi[i]+fi[i+1])
                suma = suma + trapecio
                cuenta_h[i:i+1] = [1] # usar trapecio
                suma_parcial.append(['trapecio',trapecio])
                i = i+1
                tramo_usado = True
                
        cuenta_h[n-1] = 0 # ultima casilla
        
        if vertabla==True: #mostrar datos parciales
            print('Fórmulas compuestas, tramos:',n-1)
            print('métodos 3:Simpson3/8, 2:Simpson1/3, 1:Trapecio, 0:usado')
            print('tramos iguales en xi:',cuenta_h)
            for unparcial in suma_parcial:
                print('',unparcial[0],':',unparcial[1])
        
        return(suma,cuenta_h)
    
    # usa función
    area,cuenta_h = simpson_compuesto(xi,fi,vertabla=True)
    
    # SALIDA
    print('tramos iguales en xi:',cuenta_h)
    print('Integral:',area)
    

    [ Concepto ] [ Algoritmo /función ] [ Gráfica ] [ diferencias entre tramos ]
    ..


    3. Gráfica de segmentos desiguales

    En la gráfica se diferencian los métodos aplicados por colores.

    Se añaden al algoritmo anterior las siguientes instrucciones,

    # GRAFICA ---------------------
    import matplotlib.pyplot as plt
    n = len(xi) # muestras
    tramos = n-1
    metodotipo = [['na','grey'],
                  ['Trapecio','lightblue'],
                  ['Simpson 1/3','lightgreen'],
                  ['Simpson 3/8','cyan']]
    # etiquetas linea
    plt.plot(xi[0],fi[0],'o', label=metodotipo[1][0],
             color=metodotipo[1][1])
    plt.plot(xi[0],fi[0],'o', label=metodotipo[2][0],
             color=metodotipo[2][1])
    plt.plot(xi[0],fi[0],'o', label=metodotipo[3][0],
             color=metodotipo[3][1])
    
    # relleno y bordes
    tipo = 0
    for i in range(0,tramos,1):
        if cuenta_h[i]>0:
            tipo = cuenta_h[i]
        plt.fill_between([xi[i],xi[i+1]],
                         [fi[i],fi[i+1]],[0,0],
                         color=metodotipo[tipo][1])        
    # Divisiones entre tramos
    for i in range(0,n,1):
        tipolinea = 'dashed' # inicia un método
        if cuenta_h[i]==0: 
            tipolinea = 'dotted' # dentro de un método
        plt.vlines(xi[i],0,fi[i],linestyle=tipolinea,
                   color='orange')
    
    # puntos de xi,fi
    plt.plot(xi,fi,marker='o',linestyle='dashed',
             color='orange',label='f(xi)')
    
    plt.axhline(0,color='black') # eje x
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.title('Integral numérico xi,fi con Fórmulas Compuestas')
    plt.legend()
    plt.tight_layout()
    plt.show()
    

    [ Concepto ] [ Algoritmo /función ] [ Gráfica ] [ diferencias entre tramos ]
    ..


    4. Diferencias entre tramos xi

    Las diferencias entre tramos pueden ser analizadas con una gráfica, usando operaciones de diferencias finitas.

    Diferencias Tramos XiResultado del algoritmo:

    d1xi: [0.22222222 0.22222222 0.22222223 0.22222222 0.22222222 0.22222222
     0.2        0.2        0.26666667]
    d2xi: [ 2.22044605e-16  9.99999972e-09 -9.99999972e-09 -2.22044605e-16
      3.33333361e-09 -2.22222233e-02 -4.44089210e-16  6.66666667e-02]
    diferencia máxima: 4.440892098500626e-16  en i: 7
    tramos iguales en xi: [3 0 0 2 0 1 2 0 1 0]

    Instrucciones en Python

    # revisa xi por
    # tramos equidistantes y no iguales
    import numpy as np
    
    # INGRESO
    xi = [1.        , 1.22222222, 1.44444444, 1.66666667,
          1.88888889, 2.11111111, 2+1/3, 2+1/3+0.2,
          2+1/3+0.4, 3.        ]
    
    # PROCEDIMIENTO
    casicero=1e-7
    # vectores como arreglo, números reales
    xi = np.array(xi,dtype=float)
    n = len(xi)
    
    # diferencias finitas en xi
    d1xi = np.diff(xi,1) # magnitud de tramos
    d2xi = np.diff(xi,2) # diferencia entre tramos
    errado = np.max(np.abs(d2xi)) # error mayor
    donde = np.argmax(np.abs(d2xi)) # donde es error mayor
    
    # revisa tramos iguales
    cuenta_h = -1*np.ones(n,dtype=int)
    i = 0
    while i<(n-1): # i<tramos, al menos un tramo
        tramo_usado = False
        
        # tres tramos iguales
        if (tramo_usado==False) and (i+3)<n:
            d2x = np.diff(xi[i:i+4],2) # diferencias entre tramos
            errado = np.max(np.abs(d2x))
            if errado<casicero: # Simpson 3/8
                cuenta_h[i:i+3] = [3,0,0]
                i = i+3
                tramo_usado==True
                
        # dos tramos iguales
        if (tramo_usado==False) and (i+2)<n:
            d2x = np.diff(xi[i:i+3],2) # diferencias entre tramos
            errado = np.max(np.abs(d2x))
            if errado<casicero: # Simpson 1/3
                cuenta_h[i:i+2] = [2,0]
                i = i+2
                tramo_usado = True
                
        # un tramo igual
        if (tramo_usado == False) and (i+1)<n:
            cuenta_h[i:i+1] = [1] # usar trapecio
            i = i+1
            tramo_usado = True
    
        cuenta_h[n-1] = 0 # ultima casilla
    
    # SALIDA
    print('d1xi:',d1xi)
    print('d2xi:',d2xi)
    print('diferencia máxima:',errado,' en i:',donde)
    
    print('tramos iguales en xi:',cuenta_h)
    
    # GRAFICA
    import matplotlib.pyplot as plt
    m = len(xi)-1
    
    plt.subplot(311) # diferencia entre tramos
    plt.stem(d2xi,linefmt=':',markerfmt ='C03')
    plt.plot(np.ones(n)*casicero,color='red',linestyle='dotted')
    plt.xlim(-0.1,m+0.1)
    plt.ylabel('diferencia tramos')
    plt.title('Diferencia entre tramos, muestras:'+str(n))
    
    plt.subplot(312) # tramos
    plt.stem(d1xi,linefmt=':',markerfmt ='C02')
    plt.xlim(-0.1,m+0.1)
    plt.ylabel('|tramo|')
    
    plt.subplot(313) # muestras
    plt.stem(xi,linefmt=':',markerfmt ='C01')
    plt.xlim(-0.1,m+0.1)
    plt.ylabel('xi')
    plt.xlabel('muestra i')
    
    plt.tight_layout()
    plt.show()
    

    [ Concepto ] [ Algoritmo /función ] [ Gráfica ] [ diferencias entre tramos ]