Cuadraturas: Métodos de Integración numéricos E. Galíndezα ∗ Electromagnetismo computacional ∗ Departamento de Física. Universidad Nacional de Colombia. Sede Bogotá. 16 de febrero de 2017
Resumen En el siguiente reporte, se abordó de una forma breve la teoría de los métods de integración de Romberg y Gauss, resaltando las ideas y ecuaciones básicas que permitieron poder construir el algoritmo en FORTRAN de ambos métodos. Comparando el método de Gauss-Hermite con el de Romberg, se concluyó que ambos métodos tienen sus ventajas particulares, y tienen una diferencia tajante en integrales de funciones que no son bien comportadas. Palabras Clave: Gauss-Hermite.
1.
Fortran, Include, Interface, cuadraturas, Romberg, Gauss-Legendre,
Introducción
Existe una amplia gama de algoritmos para calcular el valor numérico de una integral definida, el término cuadratura numérica, se usa a menudo como sinónimo de integración numérica, en especial cuando se trata de integrales unidimensionales. El problema común de todos los métodos es calcular una solución aproximada a la integral definida: Z b A= f (x)dx. (1)
pecio, cuya primera iteración resultaría: (b − a) (f (a) + f (b)). (2) 2 Es decir, el área de un trapecio de altura el ancho, h = b−a, del intervalo, y de bases, los valores extremales de la función. Sin embargo, dependiendo de la variación de f(x) en el intervalo, la aproximación puede ser bastante inexacta. Un nueva iteración, dividiendo (a, b) en N = 2n−1 intervalos equidistantes, puede traer una mejor estimación: A≈
a
La necesidad de tener una solución numérica, surge de la imposibilidad de llegar a una solución analítica, sobre todo, en integrales que necesitan amplios conocimientos en matemática avanzada pueden ser solucionadas de una manera más breve, a costo de un pequeño error relativo. Los algoritmos más usados, están basados en funciones de interpolación, en donde la función a integrar se sustituye por una de primitiva conocida, tal que compartan una serie de puntos, incluidos los extremos.
An,1 =
Siempre se puede disminuir el error relativo, tomando como ancho de la partición h0 = 0,5h; la nueva iteración, que conserva el valor de In,1 para disminuir el número de evaluaciones de la función, sería: n−1
An+1,1
Cuadratura de Romberg
N −1 X h [f (a) + f (b)] + h f (a + ih). (3) 2 i=1
2X 1 f (a + [2i − 1]h0 ). (4) = In,1 + h0 2 i=1
Ejecutando esta refinación n veces, es decir La integral (1) puede ser aproximada, en primera instancia, por el conocido método del tra- con hn = 0,5hn−1 , se obtiene un error que deα
[email protected]
Cuadraturas
2
ALGORITMO
Algoritmo pende solo de las potencias pares de N1 , resul- 2. tando un error de 0,25 en la iteración n + 1 con respecto a la iteración anterior en el término más 2.1. Include significativo; es posible hacer una mejor aproLa sentencia INCLUDE, realiza una copia del ximación, al eliminar este término, cosa que se archivo llamado, y remplaza las declaraciones consigue al sumar ambas iteraciones, dejando soque se encuentran específicadas en el archivo por lamente el error causado por la siguiente potenla sentencia.(Ver Anexo 4.1.1) cia de la expansión: 1
A
≈ In+1 +
K N2
INCLUDE ’ f i l e . i n c ’
(5)
Si el nombre de archivo está encerrado entre co(6) millas, tal como en el caso anterior, no se añade ninguna extensión al nombre, así el compilador K K 3A ≈ 4In+1,1 − In+1 + 2 − 2 (7) busca el archivo en el directorio actual, y de no N N encontrarlo, busca sucesivamente en las subcar4 1 A ≈ In+1,1 − In+1 . (8) petas. En resumen INCLUDE, crea una copia 3 3 del código fuente en el archivo, en la línea en Siguiendo esta dinámica, es posible llegar a que es llamado.(Ver Anexos 4.1) disminuir el error tanto como se quiera ([1, Apendix A]): 4A
≈ 4In+1,1 +
K N2
2.2.
Interface
Este bloque es un sentencia segura que per1 4m I − I . (9) mite al programa principal y a los subprogran+1,m n,m 4m − 1 4m − 1 mas externos interconectar adecuadamente. Un bloque de interfaz asegura que el programa de Cuadratura de Gauss llamada y el subprograma tienen el número y el Es posible pensar, a diferencia del método an- tipo de argumentos correctos; esto ayuda a que el terior, en intervalos (a, b) no equipartidos, como compilador detecte el uso incorrecto de un subplantea el método de Gauss, permitiendo esco- programa en el tiempo de compilación. Un bloger el ancho de las particiones de tal forma que que de interfaz se compone de (Ver Anexo 4.1): optimiza el algoritmo, aumentando la precisión El numero de argumentos y disminuyendo el número de iteraciones: In+1,m+1 =
A≈
N X
El tipo de cada argumento wi f (xi ).
(10)
el tipo del valor o los valores retornados por el subprograma.
i=1
En la mayoría de los casos, se tiene wi 6= wk , reduciendo el, problema a encontrar los valores 2.3. Integración numérica por el de wi , exigiendo que la ecuación (??) coincida método de Romberg con el valor de (1) para un polinomio de grado (2N − 1), con N el número de puntos a evaluar, (Ver Anexo 4.1.2, Función externa) en el intervalo (−1, 1), mediante las raices de los Primero, es necesario definir los argumentos a polinomios de Legendre (Ver Burden R. [2]), es utilizar (líneas 31,32): la función a integrar, los posible determinar estos anchos: puntos que definen el intervalo, el error relatiZ 1 Y n vo y absoluto, y el número máximo y mínimo x − xj wi = dx, (11) de iteraciones a realizar, además de las variables −1 j=1,j6=i xi − xj internas del programa(líneas 33,36). Si los dadonde x1 , x2 , ..., xn son las raíces del n-ésimo tos son coherentes, es posible continuar hacia la polinomio de Legendre; así como para este primera iteración, siendo sum el área del trapemétodo, conocido como cuadratura de Gauss- cio(líneas 53-56). Legendre, se consideró la expanción de f (x) en Luego, es necesario refinar mediante la segunpolinomios de Legendre, la base puede cmabiar, da iteración (Ecuación (4), líneas 60-67), para modificando así los wi , permitiendo añadir tra- luego, iterar n veces según (9), permitiendo, reatamiento de singularidades. lizar tantas operaciones como el mínimo se ellas 2
E. Galíndez
Cuadraturas
3
CONCLUSIONES
se haya impuesto (línea 78), para luego, compro- integral un poco más complicada: bar la presición mediante el error relativo, y que Z 2π ((1 + 2 cos(x))2 cos(2 ∗ x)) no se intente aproximar una integral más pequedx = 1,639850651, 3 + 2 cos(x) ña que el error absoluto(líneas 85,86). 0 (14) el método de Romberg, se aproxima a 1,63985 2.4. Integración numérica por el en la séptima iteración, generando un error relativo porcentual al método ‘Exacto’ del orden método de Gauss-Hermite de 10− 5 %, al contrario del método de Gauss, (Ver Anexo 4.1.3, Función externa) que genera un error del 3 %. Sin embargo, cuanComo se dijo con anterioridad, las funciones do el integrando tiene una discontinuidad en el pueden ser expandidas en cualquier base, en es- dominio, pero se sabe que la integral converge: Z π te caso, la función será expandida en sistemas x sin(x) = 4,355172181, (15) extendidos de Hermite, que son combinaciones −π 2 − 2 cos(x) de polinomios y logaritmos naturales, en la literatura [3], es fácil encontrar los valores para los el método de Romberg no converge, mientras que anchos de las divisiones y los puntos de evalua- el método de Gauss converge con un error relatición en términos de fracciones de intervalo. Para vo porcentual del orden de 10−4 . En conclusión, implementar el algoritmo, se definen las varia- si el integrando tiene singularidades y se sospebles utilizadas por la función (líneas 10-13), se- cha que el valor de la integral converje, la mejor guido de los vectores con los valores encontrados opción es usar el método de las cuadraturas de en J. Ma[3] (líneas 15-30), para declarar el ancho Gauss-Hermite, debido a que en el método de del intervalo y el complejo 0 − j0 (líneas 32-33), Gauss-Hermite, es bastante probable que en las concluyendo con n = 15 iteraciones (Ecuación evaluaciones, el integrando sea evaluado en la (10)). singularidad; el otro caso es, que en la mayoriá de integrandos bien comportados, Romberg converge con un mayor tiempo de ejecución, pero 3. Conclusiones con menos puntos que Gauss, además de que, en este caso, él método puede hacerse tan extenso Los métodos de integración numéricos son como se quiera. de gran utilidad en la resolución de ecuaciones integro-diferenciales; de los métodos estudiados, ambos pueden llegar a simplificar un cálculo analítico muy extenso en cuestión de segundos, aunque a la hora de tiempo de ejecución e iteraciones, realizar la integral por el método de Gauss, disminuye el numero de iteraciones sin afectar la precisión. Por ejemplo, calculando la conocida integral: Z
π
sin(x)dx = 2,
(12)
0
El método de Gauss-Hermite, con las 15 iteraciones programadas, resulta un valor de 2.00000, que el método de Romberg, consigue con igual precisión en 4 iteraciones. Utilizando el método de integración de Mathematica, se compararon ambos métodos con la integral: Z
1
exp(x2 )dx = 0,74682413,
(13)
0
Ambos resultados, tienen como resolución 0,746824, pero el método de Romberg lo alcanzó a la cuarta iteración. Realizando ahora, una 3
E. Galíndez
Cuadraturas
4. 4.1.1.
3
5
7
9
11
13
15
17
19
! La ! La ! La real
ANEXOS
Anexos
4.1.
1
4
Códigos de los programas Funciones
f u n c i ón f ( x ) r e t o r n a e l v a l o r de s i n ( x ) f u n c i ón g ( x ) r e t o r n a e l v a l o r de e ( x ) f u n c i ón z f u n ( x ) r e t o r n a e l v a l o r de s i n ( x ) FUNCTION f ( x ) IMPLICIT NONE REAL : : x f= SIN ( x )
END ! REAL FUNCTION g ( x ) IMPLICIT NONE REAL : : x g=EXP( x ) END ! COMPLEX FUNCTION z f u n ( x ) IMPLICIT NONE REAL : : x z f u n= SIN ( x ) END
file.inc 4.1.2.
Cuadratura de Romberg
Archivo Principal 2
4
6
8
10
12
! INCLUDE ’ f i l e . i n c ’ ! PROGRAM CuadraturaR IMPLICIT NONE ! INTERFACE REAL FUNCTION romberg ( f , a , b , r e r r , a e r r , maxp , minp ) REAL, INTENT( IN ) : : a , b , r e r r , a e r r INTEGER, INTENT( IN ) : : maxp , minp EXTERNAL f END FUNCTION romberg END INTERFACE
14
16
18
20
22
24
26
28
REAL : : a , b , r e r r , a e r r , ans INTEGER : : maxp , minp EXTERNAL f WRITE( ∗ , ∗ ) ’ I n i c i o d e l i n t e r v a l o \ t ’ READ( ∗ , ∗ ) a WRITE( ∗ , ∗ ) ’ F i n a l d e l i n t e r v a l o \ t ’ READ( ∗ , ∗ ) b WRITE( ∗ , ∗ ) ’ E r r o r r e l a t i v o \ t ’ READ( ∗ , ∗ ) r e r r WRITE( ∗ , ∗ ) ’ E r r o r a b s o l u t o \ t ’ READ( ∗ , ∗ ) a e r r WRITE( ∗ , ∗ ) ’ Maximo de e v a l u a c i o n e s \ t ’ READ( ∗ , ∗ ) maxp WRITE( ∗ , ∗ ) ’ Minimo de e v a l u a c i o n e s \ t ’ READ( ∗ , ∗ ) minp
30
32
ans=romberg ( f , a , b , r e r r , a e r r , maxp , minp ) WRITE( ∗ , ∗ ) ’ El v a l o r de l a i n t e g r a l e s= ’ , ans END PROGRAM CuadraturaR
MainRomberg.f90
4
E. Galíndez
Cuadraturas
4
ANEXOS
Función externa 1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
c c c c c c c c c c c c c c c c c c c c c c c c c c c c
31
33
35
37
c
39
41
43
c
45
47
49
51
c c c
55
59
f rerr aerr maxp minp
r
ier
e s una f u n c i ón r e a l de l a forma f ( x ) , donde x e s una v a r i a b l e r e a l ( f debe s e r d e c l a r a d a como e x t e r n a l en e l programa que l a l l a m e ) . e s e l e r r o r r e l a t i v o d e s e a d o ( r e r r debe e s t a r comprendido e n t r e 0.00001 y 1.0) . e s e l e r r o r a b s o l u t o d e s e a d o ( s i e l v a l o r e s t i m a d o c a e por d e b a j o de e s t a c a n t i d a d , e l a l g o r i t m o s e t e r m i n a ) . e s t a b l e c e e l número má ximo de e v a l u a c i o n e s de l a f u n c i ón como 1+2∗∗(maxp−1) . e s t a b l e c e e l número mí nimo de e v a l u a c i o n e s de l a f u n c i ón como 1+2∗∗(minp−1) ( podr í a s e r n e c e s a r i o en c a s o de que e l integrando o s c i l e v a r i a s veces dentro del i n t e r v a l o ) . e s un a r r e g l o r e a l de l o n g i t u d i n i c i a l ’ maxp ’ ( p r o p o r c i o n a d o d e s d e e l programa p r i n c i p a l ) . De v u e l t a , e s t e a r r e g l o c o n t i e n e l a ú l t i m a f i l a de e x t r a p o l a c i o n e s d e l mé todo de Romberg . e s 0 s i l a p r e c i s i ón d e s e a d a f u e a l c a n z a d a e s 1 s i r e r r e s t a f u e r a de rango e s 2 s i maxp < 1 ó maxp < minp e s 3 s i l a p r e c i s i ón d e s e a d a no f u e a l c a n z a d a
Autor : J . D. Baena Fecha : 21 e n e r o 2017 D e c l a r a c i ón de v a r i a b l e s i m p l i c i t none real f , a , b , rerr , aerr i n t e g e r maxp , minp r e a l r ( maxp ) integer ier r e a l sum , d e l , x , r l a s t , w integer i , j ,k ,n e r r o r 2 : s i e l maximo de i t e r a c i o n e s e s menor a 1 o e l minimo de e s mayor a e l maximo i f ( ( maxp . l t . 1 ) . o r . ( minp . g t . maxp ) ) then i e r =2 romberg = ( 0 . , 0 . ) go t o 30 endif E r r o r 1 : El e r r o r r e l a t i v o e s menor a 0 . 0 0 0 0 1 o mayor que 1 i f ( ( r e r r . l t . 0 . 0 0 0 0 1 ) . o r . ( r e r r . g t . 1 . ) ) then i e r =1 romberg = ( 0 . , 0 . ) go t o 30 endif i e r =0
c c c
iteraciones
Primera e s t i m a c i ón de l a i n t e g r a l , á r e a d e l p r i m e r t r a p e c i o sum =0 .5∗( f ( a )+f ( b ) ) r ( 1 )=sum ∗ ( b−a ) w r i t e ( ∗ , ∗ ) ’ Primera e s t i m a c i ón de l a i n t e g r a l : n=1
53
57
r e a l f u n c t i o n romberg ( f , a , b , r e r r , a e r r , maxp , minp ) −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− Devuelve e l v a l o r e s t i m a d o de l a i n t e g r a l de una f u n c i ón r e a l f s o b r e e l i n t e r v a l o r e a l ( a , b ) , u t i l i z a n d o e l metodo de Romberg .
’ , r (1)
R e f i n a m i e n t o de l a e s t i m a c i ón u t i l i z a n d o l a r e g l a de l o s t r a p e c i o s do 20 k=2 ,maxp r l a s t=r ( 1 ) d e l =(b−a ) / f l o a t ( n ) x=a −0.5∗ d e l do 10 i =1 ,n 10 sum=sum+f ( x+d e l ∗ f l o a t ( i ) ) n=n∗2 r ( k ) =0.5∗ d e l ∗sum
61
63
65
67
c
5
E. Galíndez
Cuadraturas
69
c c
4
ANEXOS
update Romberg e x t r a p o l a t i o n a r r a y w=4. do 15 j =2 ,k r ( k+1− j ) =(w∗ r ( k+2− j )−r ( k+1− j ) ) / (w−1) w=4.∗w 15 c o n t i n u e
71
73
75
c romberg=r ( 1 ) i f ( k . l t . minp ) go t o 20
77
79
81
83
c c c c c c
87
91
i f change i n e s t i m a t e i s w i t h i n d e s i r e d range , o r i f t h e e s t i m a t e f a l l s below t h e a b s o l u t e e r r o r d e s i r e d , r e t u r n i f ( abs ( romberg ) . l t . a e r r ) r e t u r n i f ( abs ( ( romberg−r l a s t ) / romberg ) . l t . r e r r ) r e t u r n 20 c o n t i n u e i e r =3
85
89
check i n t e g r a l estim ate :
c c c
error 30 w r i t e ( ∗ , ∗ ) return end
93
’ E r r o r en romberg . i e r = ’ , i e r
romberg.f
4.1.3.
Cuadratura de Gauss
Archivo Principal INCLUDE ’ f i l e . i n c ’ 2
4
6
8
10
PROGRAM CuadraturaGauss IMPLICIT NONE INTERFACE COMPLEX FUNCTION QGLL15( zfun , a , b ) REAL, INTENT( IN ) : : a , b EXTERNAL z f u n END FUNCTION QGLL15 END INTERFACE REAL : : a , b , ans
12
14
16
18
20
EXTERNAL z f u n WRITE( ∗ , ∗ ) ’ I n i c i o d e l i n t e r v a l o \ t ’ READ( ∗ , ∗ ) a WRITE( ∗ , ∗ ) ’ F i n a l d e l i n t e r v a l o \ t ’ READ( ∗ , ∗ ) b ans=QGLL15( zfun , a , b ) WRITE( ∗ , ∗ ) ’ El v a l o r de l a i n t e g r a l e s= ’ , ans END PROGRAM CuadraturaGauss
MainCGauss.f90 Función externa complex f u n c t i o n QGLL15( zfun , a , b )
1
3
5
7
9
C C C C C C C C
This f u n c t i o n computes t h e i n t e g r a l o f t h e f u n c t i o n ’ f u n ’ from a t o b u s i n g o r d e r 15 Gaussian q u a d r a t u r e but with p o l y n o m i a l s and n a t u r a l l o g a r i t h m s a s b a s i s f u n c t i o n s . See Ma, Rokhlin , and Wandzura , " G e n e r a l i z e d Gaussian Quadrature R u l e s f o r Systems o f A r b i t r a r y F u n c t i o n s , " R e s e a r c h Report , Yale U n i v e r s i t y . complex asum , z f u n
6
E. Galíndez
Cuadraturas
4
ANEXOS
r e a l x ( 1 5 ) , w( 1 5 ) integer n e x t e r n a l zfun
11
13
C 15
+ + + + + + +
17
19
21
23
+ + + + + + +
25
27
29
31
C x l e n g t h = b−a asum = ( 0 . , 0 . )
33
35
37
39
data x / . 1 0 5 7 8 4 5 4 8 4 5 8 6 2 9 e − 3 , . 1 5 6 6 2 4 3 8 3 6 1 6 7 8 2 e −2 , . 7 5 9 5 2 1 8 9 0 3 2 0 7 0 9 e − 2 , . 2 2 8 3 1 0 6 7 3 9 3 9 8 6 2 e −1 , . 5 2 3 8 8 6 3 0 1 5 6 8 2 0 0 e − 1 , . 1 0 0 7 5 8 6 8 5 2 0 1 2 1 3 e0 , . 1 7 0 7 4 0 7 6 8 8 4 9 9 4 3 e0 , . 2 6 2 5 9 1 2 0 6 1 1 8 9 9 3 e0 , . 3 7 3 5 3 6 5 0 5 1 8 4 5 5 8 e0 , . 4 9 7 7 4 6 3 5 8 4 1 4 5 3 3 e0 , . 6 2 6 7 8 9 0 3 1 3 9 2 3 7 3 e0 , . 7 5 0 5 1 6 1 0 3 4 6 1 4 0 8 e0 , . 8 5 8 2 5 5 3 3 5 2 0 7 8 6 1 e0 , . 9 4 0 1 4 1 2 9 1 2 1 2 3 4 6 e0 , . 9 8 8 4 0 1 5 9 5 9 8 6 3 4 2 e0 / data w/ . 4 0 3 2 1 7 7 2 4 6 4 8 4 6 0 e − 3 , . 3 0 6 2 9 7 8 4 3 4 7 8 7 e −2 , . 9 7 8 4 2 1 2 1 1 8 7 6 6 1 5 e − 2 , . 2 1 5 5 8 7 5 2 2 2 5 5 8 1 3 e −1 , . 3 8 3 2 3 0 6 7 3 7 0 8 8 9 2 e − 1 , . 5 8 8 9 8 1 9 9 0 2 6 3 0 0 4 e −1 , . 8 1 1 1 7 0 2 9 9 3 9 2 5 9 5 e − 1 , . 1 0 2 1 2 2 1 0 1 9 7 2 0 6 9 e0 , . 1 1 8 7 8 9 0 5 9 0 3 0 4 0 1 e0 , . 1 2 8 2 1 0 3 1 6 4 4 6 6 9 4 e0 , . 1 2 8 1 6 3 3 2 7 4 1 7 0 9 3 e0 , . 1 1 7 4 8 9 4 6 5 8 8 8 4 9 2 e0 , . 9 6 3 2 3 0 1 8 5 6 9 5 9 0 4 e − 1 , . 6 6 1 3 4 5 3 9 8 3 1 8 9 3 4 e −1 , . 2 9 6 2 0 7 1 4 0 0 3 5 3 5 5 e −1/
10
do 10 n = 1 , 1 5 asum = w( n ) ∗ z f u n ( x l e n g t h ∗x ( n )+a )+asum continue QGLL15 = x l e n g t h ∗asum return end
cuadratura_gauss.f
7
E. Galíndez
Cuadraturas
REFERENCIAS
Referencias
gage Learning, 9 edition. pág. 231.
[1] R. Mittra . F. Peterson, S. L. Ray. Computational Methods for Electromagne- [3] J. Ma, V. Rokhlin, and S. Wandzura. Generalized gaussian quadrature rules for systems tics. IEEE, 1997. of arbitrary functions. SIAM Journal on Nu[2] Faires J. Burden R. Numerical analysis. Cen- merical Analysis, 33(3):971–996, 1996.
8
E. Galíndez
