Tema 5

Métodos Numéricos Aplicados I

Integración numérica
 Índice
Esquema. . . . . . . . . . . . . . . . . . . . . . . 2

Ideas clave . . . . . . . . . . . . . . . . . . . . . . 3

5.1 Introducción y objetivos . . . . . . . . . . . . . 3

5.2 Cuadratura numérica y polinomio de Lagrange . . . . 5

5.3 Fórmulas cerradas de Newton-Cotes . . . . . . . . 7

5.4 Fórmulas abiertas de Newton-Cotes . . . . . . . . 14

5.5 Cuadratura de Gauss. . . . . . . . . . . . . . . 18

5.6 Integración múltiple . . . . . . . . . . . . . . . 29

 Esquema

INTEGRACIÓN NUMÉRICA

N O D O S E Q U I E S PA C I A D O S N O D O S N O E Q U I E S PA C I A D O S
Newton-Cotes Cuadratura de Gauss

Métodos cerrados Métodos abiertos

Gauss-Legendre Gauss-Hermite
Trapecios Punto medio
Simpson
Gauss-Chebyshev Gauss-Laguerre

I N T E G R A C I Ó N S I M P L E Y M Ú LT I P L E

Métodos Numéricos Aplicados I

2
Tema 5. Esquema
 Ideas clave

5.1 Introducción y objetivos

En ingeniería, el trabajo W necesario para trasladar un cuerpo de una posición a otra

se define como
ED
~ ~
dW = F , dl ,

donde F~ es el vector fuerza, y d~l es el vector desplazamiento, como se ilustra en la

Figura 1.

Figura 1: Representación del trabajo

Se trata de un producto escalar entre dos vectores, por lo que

dW = F dl cos(α),

donde F = |F~ | y dl = |d~l|. Cuando la expresión anterior depende de x, podemos

obtener el trabajo como
Z b
W = F (x) cos(α(x)) dx. (1)
a

En ocasiones, F (x) tiene una expresión analítica que permite obtener el valor de W

Métodos Numéricos Aplicados I

3
Tema 5. Ideas clave
fácilmente; sin embargo, habitualmente nos encontramos con que los valores de F (x)
y α(x) vienen tabulados (como se representa en la Tabla 1) o no dan lugar a una ex-
presión con una integral conocida.

x F (x) α(x)
0.00 0.00 0.50
1.52 40.04 1.40
3.04 57.83 0.75
4.56 62.28 0.90
6.08 46.71 1.30
7.60 53.38 1.48
9.12 22.24 1.50

Tabla 1: Datos tabulados de F (x) y α(x)

La técnica de cuadratura numérica consiste en calcular una integral a partir de unos

valores discretos; estos pueden tener origen en el valor del integrando en los dife-
rentes nodos o en valores tabulados. En general, la cuadratura numérica consiste en
obtener Z b n
X
I= f (x) dx ≈ ai f (xi ).
a i=0

Cuando los nodos están equiespaciados, trabajaremos con las conocidas fórmulas de
Newton-Cotes. Sin embargo, si los nodos no están equiespaciados, hablaremos de las
cuadraturas de Gauss.

A lo largo del tema, trabajaremos principalmente con integrales simples. En la última

parte del tema veremos cómo los conceptos de la integración simple son extrapolables
a la integración múltiple en todos los casos.

Los objetivos que trataremos de alcanzar en este tema serán los siguientes:

I Conocer la expresión general de las técnicas de cuadratura

Métodos Numéricos Aplicados I

4
Tema 5. Ideas clave
I Comprender e implementar las fórmulas derivadas de las expresiones de Newton-
Cotes, cerradas y abiertas

I Comprender e implementar los casos particulares de la cuadratura de Gauss

I Aplicar la integración numérica sobre integración múltiple

Algunas funciones Matlab utilizadas en este tema

I repmat(x,a,b). Repite el vector x en a filas y b columnas.

5.2 Cuadratura numérica y polinomio de

Lagrange

Nuestro objetivo es obtener el valor de

Z b
I= f (x) dx (2)
a

utilizando la técnica de la cuadratura numérica. Los métodos que vamos a ver a lo

largo del tema se basan en los polinomios de Lagrange ln (x). Sabemos que podemos
expresar
f (x) = ln (x) + (x), (3)

donde
n
f (n+1) (ξ(x)) Y
(x) = (x − xi ),
(n + 1)! i=0

y el polinomio de interpolación de Lagrange ln (x) viene dado por

n n
X Y (x − xj )
ln (x) = Li (x)f (xi ), Li (x) = . (4)
i=0 j=0
(x i − x j )
j6=i

Métodos Numéricos Aplicados I

5
Tema 5. Ideas clave
Reemplazando (3) en (2), obtenemos
Z b
I = (ln (x) + (x)) dx =
a

Z b (X
n n
)
f (n+1) (ξ(x)) Y
= Li (x)f (xi ) + (x − xi ) dx = (5)
a i=0
(n + 1)! i=0

n Z b Z b n
X 1 (n+1)
Y
= f (xi ) Li (x) dx + f (ξ(x)) (x − xi ) dx.
i=0 a (n + 1)! a i=0

De este modo,
n
X Z b
I≈ f (xi ) Li (x) dx, (6)
i=0 a

denominándose error de cuadratura al término

Z b n
1 (n+1)
Y
f (ξ(x)) (x − xi ) dx.
(n + 1)! a i=0

En lo sucesivo del tema, veremos los métodos que surgen del uso de la expresión (6).

En función del número de términos, podremos obtener el valor de I como se muestra

en la Figura 2.

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x x
(a) (b)

Figura 2: Diferentes aproximaciones de I

Métodos Numéricos Aplicados I

6
Tema 5. Ideas clave
Cuando trabajamos con nodos equiespaciados estaremos en el caso de fórmulas de
Newton-Cotes, que desarrollaremos en los apartados 5.3 y 5.4. En el caso de nodos
equiespaciados, en el apartado 5.5 nos adentraremos en las fórmulas de cuadratura
de Gauss.

5.3 Fórmulas cerradas de Newton-Cotes

En este apartado nos centraremos en el caso en que los nodos están equiespaciados
y, además, se toman todos los nodos del intervalo x ∈ [a, b], con la situación de la
Figura 3, con 𝑥𝑖 𝑥𝑖+1
b−a
𝑎 𝑥 h= .
n

𝑎 𝑎+ℎ 𝑎 + 2ℎ 𝑏−ℎ 𝑏
Figura 3: Esquema de nodos para fórmulas cerradas de Newton-Cotes
𝑥𝑖−1 𝑥𝑖 𝑥𝑖+1 𝑎 𝑎+ℎ 𝑎 + 2ℎ

Para desarrollar cada una de las fórmulas de integración, recurriremos a particulari-

zar la expresión (6) para algunos valores particulares de n. La Tabla 2 recoge estas
expresiones en su versión simple, algunos de cuyos desarrollos veremos a lo largo del
apartado. En todos los casos es de aplicación el Teorema 1.
Rb
Método Aproximación de a
f (x) dx 
3
Trapecios h
2
(f (a) + f (b)) − h12 f 00 (ξ)
h 5
f (a) + 4f a+b − h90 f IV (ξ)



Simpson 3 2
+ f (b)
3h 2a+b a+2b 5
− 3h




Simpson 3/8 8
f (a) + 3f 3
+ 3f 3
+ f (b) 80
f IV (ξ)
2h 7
7f (a) + 32f 4 + 12f 2 + 32f a+3b
3a+b a+b
− 8h





Milne 45 4
+ 7f (b) 945
f V I (ξ)

Tabla 2: Fórmulas cerradas simples de Newton-Cotes

Métodos Numéricos Aplicados I

7
Tema 5. Ideas clave
 Teorema 1: Error de las fórmulas cerradas de Newton-Cotes

Pn
Sea i=1 ai f (xi ) la fórmula de Newton-Cotes cerrada de n + 1 puntos donde
b−a
x0 = a, xn = b y h = n
. Entonces existe un ξ ∈ (a, b) tal que

I para valores pares de n y si f ∈ C n+2 [a, b] se cumple que

b n (n+2) n
hn+3
Z X Z
f (x) dx = ai f (xi )+ (ξ)(n + 2)! τ 2 (τ −1) · · · (τ −n) dτ.
a i=0
f 0

I para valores impares de n y si f ∈ C n+1 [a, b] se cumple que

b n (n+1) n
hn+2
Z X Z
f (x) dx = ai f (xi ) + (ξ)(n + 1)! τ (τ − 1) · · · (τ − n) dτ.
a i=0
f 0

Método de Trapecios

Se denomina método de Trapecios al caso en que (6) tiene n = 2. En este sentido,

Z b Z b
I ≈ f (a) L0 (x) dx + f (b) L1 (x) dx =
a a

Z b Z b
f (a) f (b) (7)
= − (x − b) dx + (x − a) dx =
h a h a

h
= (f (a) + f (b)) .
2

El nombre del método de Trapecios viene dado porque (7) es la expresión del área de
un trapecio.

Se puede demostrar que el error es

h3 00
=− f (ξ(x)).
12

Métodos Numéricos Aplicados I

8
Tema 5. Ideas clave
Para disminuir el error, se divide el intervalo en subintervalos y se aplica el método de
trapecios (7) sobre cada subintervalo, como se ilustra en la Figura 4.

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 1 2 0 1 2
x x
(a) (b)

Figura 4: Método de trapecios (a) con el paso original, (b) dividiendo el intervalo en
subintervalos

De este modo,

n−1
hX
I ≈ (f (xi ) + f (xi+1 )) =
2 i=0
(8)
h
= (f (x0 ) + 2f (x1 ) + 2f (x2 ) + · · · + 2f (xn−1 + f (xn )) .
2

A este método se le conoce como fórmula de trapecios compuesta. Su expresión del

error es
h2
 = − (b − a)f 00 (ξ(x)).
12

A continuación se muestra la implementación del método de trapecios en Matlab.

Métodos Numéricos Aplicados I

9
Tema 5. Ideas clave
 Trapecios.m

function I= trapecios (f ,a,b,n)

Métodos Numéricos Aplicados I

10
Tema 5. Ideas clave
 % I= trapecios (f,a,b,n) obtiene la integral de f(x)
% con la fórmula de trapecios compuesta.
h = (b-a)/n;
x = a:h:b;
pesos = [1 2* ones (1 ,n -1) 1];
I=h /2* sum ( pesos. *f(x));
end

Método de Simpson

Cuando tomamos (6) con n = 3, tenemos el método de Simpson.

Z b  Z b Z b
a+b
I ≈ f (a) L0 (x) dx + f L1 (x) dx + f (b) L2 (x) dx =
a 2 a a
(9)
   
h a+b
= f (a) + 4f + f (b) ,
3 2

siendo su error
h5 IV
=− f (ξ(x)).
90

Al igual que en el método de trapecios, podemos disminuir el error si tomamos subin-

tervalos. Por tanto,

n−2
hX
I ≈ (f (xi ) + 4f (xi+1 ) + f (xi+2 )) =
3 i=0

(10)
h
= (f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + · · · +
3
+2f (xn−2 ) + 4f (xn−1 ) + f (xn )) ,

dando lugar a la fórmula de Simpson compuesta. En este caso, el error pasa a ser

h4
=− (b − a)f IV (ξ(x)).
180

Métodos Numéricos Aplicados I

11
Tema 5. Ideas clave
A continuación se muestra la implementación del método de Simpson en Matlab.

Simpson.m

function I= simpson (f,a,b,n)

% I= simpson (f,a,b,n) obtiene la integral de f(x)
% con la fórmula de Simpson compuesta.
h = (b-a)/n;
x = a:h:b;
pesos = ones (1 ,n +1) ;
pesos ([Link] n) = 4; pesos ([Link]n -1) = 2;
I=h /3* sum ( pesos. *f(x));
end

Ejemplo 1. Calcula Z π
2
I= sin(x)e−x dx,
0

utilizando los métodos de Trapecios y Simpson tomando 4 y 8 subintervalos.

Para cada caso, calcula el error cometido sabiendo que el resultado analítico es
π
1 − e− 2
.
2

Interpreta los resultados.

Ejecutamos en Matlab

>> f=@(x) sin (x).* exp (-x);

>> a =0; b=pi /2; n1 =4; n2 =8;
>> IT1 = trapecios (f,a,b,n1)
IT1 =
0 .380590604382816
>> IT2 = trapecios (f,a,b,n2)
IT2 =

Métodos Numéricos Aplicados I

12
Tema 5. Ideas clave
 0 .392182862002726
>> IS1 = simpson (f,a,b,n1)
IS1 =
0 .395839444235324
>> IS2 = simpson (f,a,b,n2)
IS2 =
0 .396046947876029
>> error = abs ([ IT1 IT2 IS1 IS2 ] -(1 - exp (-pi /2) ) /2)
error =
0 .015469607441803 0 .003877349821893 ...
0 .000220767589295 0 .000013263948590

Vemos que para un mismo método, se obtiene menor error cuando disminuimos
el tamaño de h. Además, comparando ambos métodos, las fórmulas de Simpson
tienen menor error.

Ejemplo 2. Con los datos de la Tabla 1, calcula el trabajo realizado utilizando la

expresión (1) utilizando los métodos de Trapecios y Simpson. Toma como paso
la diferencia entre los valores de x.
En los programas Trapecios.m y Simpson.m hemos tomado como valor de en-
trada la función f (x). Cuando trabajamos con datos tabulados, son estos datos
los que tienen que ser el parámetro de entrada, y hacer ligeras modificaciones
sobre los programas, que llamaremos TrapeciosDatos.m y SimpsonDatos.m.

Para obtener el trabajo realizado, ejecutamos en Matlab

>> x =0:1 .52 :9 .12 ;

>> F =[0 40 .04 57 .83 62 .28 46 .71 53 .38 22 .24 ];
>> alpha =[ .5 1.4 .75 .9 1.3 1 .48 1.5 ];
>> datos =F.* cos ( alpha );
>> IT= trapeciosDatos ( datos ,x)

Métodos Numéricos Aplicados I

13
Tema 5. Ideas clave
 IT =
1 .610507484870318e +02
>> IS= simpsonDatos ( datos ,x)
IS =
1 .583980289247424e +02

5.4 Fórmulas abiertas de Newton-Cotes

+1
Si en el apartado anterior desarrollamos algunas fórmulas de Newton-Cotes cerradas,
en este apartado nos centraremos en las fórmulas abiertas. En este caso, los nodos
primero y último no se toman, y se utiliza la nomenclatura de la Figura 5. Ahora,

ℎ 𝑎 + 2ℎ 𝑏−ℎ 𝑏 b−a
h= .
n−2

𝑖 𝑥𝑖+1
𝑥−1 𝑥0 𝑥1 𝑥𝑛 𝑥𝑛+1

𝑎 𝑎+ℎ 𝑎 + 2ℎ 𝑏−ℎ 𝑏
Figura 5: Esquema de nodos para fórmulas abiertas de Newton-Cotes

De nuevo partiremos de (6) para desarrollar las fórmulas de integración. La Tabla 3 re-
coge estas expresiones en su versión simple, desarrollando el método de punto medio
en este apartado. Para todos los casos se puede aplicar el Teorema 2.
Rb
Método Aproximación de a
f (x) dx 
h3
a+b
f 00 (ξ)


n = 0: Punto medio (b − a)f 2 3
b−a 2a+b a+2b 3h3 00




n=1 2
f 3
+ f 3 4
f (ξ)
b−a 14h5 IV
2f 4 + f 2 + 2f a+3b
3a+b a+b





n=2 3 4 45
f (ξ)

Tabla 3: Fórmulas abiertas simples de Newton-Cotes

Métodos Numéricos Aplicados I

14
Tema 5. Ideas clave
 Teorema 2: Error de las fórmulas abiertas de Newton-Cotes

Pn
Sea i=1 ai f (xi ) la fórmula de Newton-Cotes cerrada de n + 1 puntos donde
b−a
x−1 = a, xn+1 = b y h = n+2
. Entonces existe un ξ ∈ (a, b) tal que

I para valores pares de n y si f ∈ C n+2 [a, b] se cumple que

b n (n+2) n+1
hn+3
Z X Z
f (x) dx = ai f (xi )+ (ξ)(n + 2)! τ 2 (τ −1) · · · (τ −n) dτ.
a i=0
f −1

I para valores impares de n y si f ∈ C n+1 [a, b] se cumple que

b n (n+1) n+1
hn+2
Z X Z
f (x) dx = ai f (xi )+ (ξ)(n + 1)! τ (τ −1) · · · (τ −n) dτ.
a i=0
f −1

Método de Punto Medio

Tomando (6) con n = 0, tendremos los nodos

a+b
x−1 = a, x0 = , x1 = b.
2

Por tanto,
Z b  
a+b
I ≈ f (x0 ) dx = (b − a)f (x0 ) = (b − a)f . (11)
a 2

siendo su error
h3 00
=− f (ξ(x)).
3

Para disminuir el error, tomamos subintervalos. Como el método de punto medio abar-
ca los puntos anterior y posterior, tomaremos solo los puntos pares para no solapar;

Métodos Numéricos Aplicados I

15
Tema 5. Ideas clave
además, el número de subintervalos n deberá ser par. Por tanto,

n/2
X
I ≈ (x2i+1 − x2i−1 )f (x2i ) =
i=0
(12)
= 2h (f (x0 ) + f (x2 ) + · · · + f (xn−2 ) + f (xn )) .

dando lugar a la fórmula de punto medio compuesta, siendo su error

b − a 2 00
= h f (ξ).
6

A continuación se muestra la implementación del método de Punto Medio en Matlab.

PuntoMedio.m

function I= puntoMedio (f,a,b,n)

% I= puntoMedio (f,a,b,n) obtiene la integral de f(x)
% con la fórmula de punto medio compuesta.
h = (b-a)/(n +2) ;
x = a+h:h:b-h;
I =2* h* sum (f(x (1:2 :end)));
end

Ejemplo 3. Calcula Z π
2
I= sin(x)e−x dx,
0

utilizando el método de Punto Medio tomando 4 y 8 subintervalos. Para cada

caso, calcula el error cometido sabiendo que el resultado analítico es
π
1 − e− 2
.
2

Interpreta los resultados.

Ejecutamos en Matlab

Métodos Numéricos Aplicados I

16
Tema 5. Ideas clave
>> f=@(x) sin (x).* exp (-x);
>> a =0; b=pi /2; n1 =4; n2 =8;
>> IPM1 = puntoMedio (f,a,b,n1)
IPM1 =
0 .409710138362011
>> IPM2 = puntoMedio (f,a,b,n2)
IPM2 =
0 .401008515076530
>> error = abs ([ IPM1 IPM2 ] -(1 - exp (-pi /2) ) /2)
error =
0 .013649926537392 0 .004948303251911

Se evidencia la reducción del error tomando un mayor número de subintervalos.

Ejemplo 4. Con los datos de la Tabla 1, calcula el trabajo realizado utilizando la

expresión (1) utilizando el método de Punto Medio.
Debemos interpretar este problema a partir de la comprensión del método de
Punto Medio. Este método calcula el área del rectángulo que tiene como base la
diferencia entre tres puntos consecutivos y de altura el valor de la función en el
punto intermedio. Si queremos utilizar el valor de la función en todos los puntos,
podemos determinar unos nodos intermedios por los que no se calcula la función.
Así, los nodos serán

1.52
xi = i, i = −1, 0, . . . , 13,
2

1.52
siendo h = 2
, como refleja la siguiente figura.

Métodos Numéricos Aplicados I

17
Tema 5. Ideas clave
 40

30

20

10

0
-1.52 0 1.52 3.04 4.56 6.08 7.6 9.12 10.64
x

Aplicando (12), obtenemos

I = 2h (f (x0 ) + f (x2 ) + · · · + f (x10 ) + f (x12 )) = 162.24638.

5.5 Cuadratura de Gauss

La cuadratura de Gauss es un tipo de cuadratura para nodos no equiespaciados. Esta

cuadratura pasa por la obtención de los nodos x1 , x2 , · · · , xn y de los coeficientes
c1 , c2 , . . . , cn que minimizan el error obtenido en la aproximación

Z b n
X
w(x)f (x) dx ≈ ci f (xi ),
a i=1

donde a la función w(x) se la denomina función peso, y cumple w(x) > 0, x ∈ [a, b].
Para la obtención de estos coeficientes, se deben verificar las siguientes condiciones:

1. Hay que determinar unos polinomios p1 (x), . . . , pn (x) tales que

hpi , pj i = 0, i 6= j.

2. Los nodos x1 , . . . , xn son las raíces del polinomio pn (x).

A partir de estas dos condiciones, se pueden obtener los coeficientes c1 , . . . , cn y, por

Métodos Numéricos Aplicados I

18
Tema 5. Ideas clave
tanto, la fórmula de cuadratura.

El error cometido al realizar la aproximación de la integral por las diferentes cuadra-

turas de Gauss es

b
f (2n) (ξ(x))
Z
= p2n (x)w(x) dx, ξ ∈ (a, b).
(2n)! a

Las diferentes cuadraturas vienen definidas por la función w(x).

Cuadratura de Gauss-Legendre

La fórmula de cuadratura de Gauss-Legendre utiliza la función peso w(x) = 1. El in-

tervalo [a, b] = [−1, 1], pero veremos más adelante que con un cambio de variable
sencillo se puede extender a cualquier intervalo [a, b]. Por tanto, las integrales que
aproximaremos serán de la forma
Z 1
f (x) dx.
−1

Los polinomios pi (x) de la cuadratura de Gauss-Legendre son

p0 (x) = 1,
p1 (x) = x,
1
pk+1 (x) = k+1
[(2k + 1)xpk (x) − kpk−1 (x)] , k = 1, 2, . . . , n − 1.

Los polinomios anteriores cumplen con la condición de ortogonalidad. Además, las

raíces de pn (x) son
   
1 1 4k − 1
xi = 1− 2 + 3 cos π .
8n 8n 4n + 2

Los coeficientes ci se obtienen a partir de la expresión

2
ci = .
(1 − x2i )(p0n (xi ))2

Métodos Numéricos Aplicados I

19
Tema 5. Ideas clave
La Tabla 4 recoge los valores de los nodos y de los coeficientes para diferentes valores
de n.

i 1 2 3 4 5
xi -0.577350 0.577350 - - -
n=2
ci 1.000000 1.000000 - - -
xi 0.000000 -0.774597 0.774597 - -
n=3
ci 0.888889 0.555556 0.555556 - -
xi -0.339981 -0.861136 0.339981 0.861136 -
n=4
ci 0.652145 0.347855 0.652145 0.347855 -
xi 0.000000 -0.538469 -0.906180 0.538469 0.906180
n=5
ci 0.568889 0.478629 0.236927 0.478629 0.236927

Tabla 4: Nodos y coeficientes de la cuadratura de Gauss-Legendre

Por ejemplo, si tuviéramos que obtener la expresión de la fórmula de cuadratura de

Gauss-Legendre para n = 3, tendríamos que
Z 1
f (x) dx = c1 f (x1 ) + c2 f (x2 ) + c3 f (x3 ) =
−1

= 0.888889f (0) + 0.555556f (−0.774597) + 0.555556f (0.774597).

� Accede al vídeo: Obtención de nodos y coeficientes de la cuadratura de Gauss-

Legendre en Matlab

Para transformar el intervalo [−1, 1] en el intervalo genérico [a, b] utilizamos el cambio

de variable
b−a b+a b−a
x= y+ , dx = dy, (13)
2 2 2
por lo que podremos aproximar

b n  
b−aX b−a
Z
b+1
I= f (x) dx ≈ ci f xi + .
a 2 i=1 2 2

Métodos Numéricos Aplicados I

20
Tema 5. Ideas clave
 Ejemplo 5. Obtén el valor de
Z 1.5
2
I= e−x dx
1

con la cuadratura de Gauss-Legendre utilizando n = 2 y n = 3.

Transformamos el intervalo [1, 1.5] en el [−1, 1] utilizando (13) como

1.5 − 1 1.5 + 1 y+5 y

x= y+ = , dx = ,
2 2 4 4

de modo que Z 1.5 Z 1

−x2 1 (y+5)2
I= e dx = e− 16 dy.
1 4 −1

Con n = 2 tenemos que

Z 1  
1 −
(y+5)2 1 −
(−0.577350+5)2
−
(0.577350+5)2
e 16 dy ≈ e 16 +e 16 = 0.109400.
4 −1 4

Para n = 3,
Z 1
(y+5)2
1
4
e− 16 dy ≈
−1 
(0+5)2 (−0.774597+5)2 (0.774597+5)2

≈ 14 0.888889e− 16 + 0.555556e− 16 + 0.555556e− 16 =
= 0.109364.

Cuadratura de Gauss-Chebyshev

La fórmula de cuadratura de Gauss-Chebyshev utiliza la función peso w(x) = √ 1 .

1−x2

De nuevo, utilizamos el intervalo [a, b] = [−1, 1], pero se puede generalizar para cual-
quier intervalo [a, b] con el cambio de variable (13). Las integrales que aproximaremos
serán de la forma Z 1
f (x)
√ dx.
−1 1 − x2

Métodos Numéricos Aplicados I

21
Tema 5. Ideas clave
Los polinomios Ti (x) de la cuadratura de Gauss-Chebyshev son

T0 (x) = 1,
T1 (x) = x,
Tk (x) = 2xTk−1 (x) − Tk−2 , k = 2, 3, . . . , n.

Los polinomios anteriores cumplen con la condición de ortogonalidad. Las raíces de

Tn (x) son  
2k − 1
xi = cos π ,
2n
mientras que los coeficientes ci tiene la expresión

π
ci = .
n

De este modo, la expresión genérica de la fórmula de cuadratura de Gauss-Chebyshev

será Z 1 n
f (x) πX
√ dx ≈ f (xi ).
−1 1 − x2 n i=1

La Tabla 5 recoge los valores de los nodos para diferentes valores de n.

i 1 2 3 4 5
n=2 xi -0.707107 0.707107 - - -
n=3 xi -0.866025 0.000000 0.866025 - -
n=4 xi -0.923880 -0.382683 0.382683 0.923880 -
n=5 xi -0.951057 -0.587785 0.000000 0.587785 0.951057

Tabla 5: Nodos de la cuadratura de Gauss-Chebyshev

Ejemplo 6. Calcular la integral

1
ex
Z
I= √ dx
−1 1 − x2

Métodos Numéricos Aplicados I

22
Tema 5. Ideas clave
con seis decimales exactos. Es decir, encontrar el valor de n tal que

2π
= eξ < 10−6 , ξ ∈ (−1, 1).
22n (2n)!

La única variable que aparece en la expresión de  es ξ, así que vamos a acotarla.

2π 2π
|eξ | ≤ e ↔ || = eξ ≤ e .
22n (2n)! 22n (2n)!

Calculamos para diferentes valores de n, obteniendo:

I n = 1 : || < 2.134933555

I n = 2 : || < 0.044477782

I n = 3 : || < 3.70648 · 10−4

I n = 4 : || < 1.65468 · 10−6

I n = 5 : || < 4.59633 · 10−9

Por tanto, para n = 5 se garantiza que el error es menor a 10−6 . Así, la integral se
obtiene como

1
ex
Z
I = √ dx ≈
−1 1 − x2

π −0.951057
+ e−0.587785 + e0 + e0.587785 + e0.951057 = 3.977463.


≈ e
5

Métodos Numéricos Aplicados I

23
Tema 5. Ideas clave
Cuadratura de Gauss-Laguerre

La fórmula de cuadratura de Gauss-Laguerre tiene como función peso a w(x) = e−x

en el intervalo [a, b] = [0, +∞), por lo que las integrales a aproximar serán de la forma
Z +∞
e−x f (x) dx.
0

Los polinomios Li (x) de la cuadratura del Gauss-Laguerre son

L0 (x) = 1,
L1 (x) = 1 − x,
Lk+2 (x) = (2k + 3 − x)Lk+1 (x) − (k + 1)2 Lk (x), k = 0, 1, . . . , n − 2.

Los polinomios anteriores son ortogonales respecto del producto escalar. Las raíces de
Ln (x) son
j2 2
 
−2 + j0k
xi = 0i 1+ ,
4hn 48h2n
siendo hn = n+ 12 y j0i la raíz i-ésima de la función J0 (x), es decir, la función de Bessel
de primera especie y orden. Los coeficientes ci se calculan como

(n!)2 xi
ci = .
L2n+1 (xi )

Ejemplo 7. Calcular la integral

Z +∞
I= e−10x sin(x) dx
0

utilizando la cuadratura de Gauss-Laguerre para n = 2 y n = 4.

Vemos que el dominio de la integral están en [0, +∞). Sin embargo, no es de la
forma Z +∞
I= e−x f (x) dx,
0

así que hagamos un cambio de variable. Si tomamos y = 10x → dy = 10 dx,

Métodos Numéricos Aplicados I

24
Tema 5. Ideas clave
entonces
Z +∞ Z +∞
−10x 1 −y
y
I= e sin(x) dx = e sin dy.
0 10 0 10

Esta última integral ya tiene la forma esperada. Calculamos los polinomios de La-
guerre hasta n = 5.

L0 (x) = 1,
L1 (x) = 1 − x,
L2 (x) = x2 − 4x + 2,
L3 (x) = −x3 + 9x2 − 18x + 6,
`L4 (x) = x4 − 16x3 + 72x2 − 96x + 24,
L5 (x) = −x5 + 25x4 − 200x3 + 600x2 + 120.

Calculamos las raíces de L2 (x) y L4 (x) ejecutando en Matlab

>> raicesL2 = roots ([1 -4 2])

raicesL2 =
3 .414213562373095
0 .585786437626905
>> raicesL4 = roots ([1 -16 72 -96 24])
raicesL4 =
9 .395070912301119
4 .536620296921135
1 .745761101158346
0 .322547689619392

Obtenemos los coeficientes ejecutando en Matlab

>> syms x
>> L3=-x ^3+9* x ^2 -18* x +6;

Métodos Numéricos Aplicados I

25
Tema 5. Ideas clave
>> c2 =( factorial (2) ) ^2* raicesL2. / ...
( double ( subs (L3 ,x, raicesL2 ))).^2
c2 =
0 .146446609406726
0 .853553390593274
>> L5=-x ^5+25* x ^4 -200* x ^3+600* x ^2+120;
>> c4 =( factorial (4) ) ^2* raicesL4. / ...
( double ( subs (L5 ,x, raicesL4 ))).^2
c4 =
0 .000539294705561
0 .038887908515005
0 .357418692437800
0 .603154104341635

Por tanto, el valor de la integral para n = 2 y n = 4 será

>> I2 =( c2 (1) * sin ( raicesL2 (1) /10) + ...

c2 (2) * sin ( raicesL2 (2) /10) ) /10
I2 =
0 .009900565097779
>> I4 =( c4 (1) * sin ( raicesL4 (1) /10) + ...
c4 (2) * sin ( raicesL4 (2) /10) + ...
c4 (3) * sin ( raicesL4 (3) /10) + ...
c4 (4) * sin ( raicesL4 (4) /10) ) /10
I4 =
0 .009900990092799

Métodos Numéricos Aplicados I

26
Tema 5. Ideas clave
Cuadratura de Gauss-Hermite

2
La fórmula de cuadratura de Gauss-Hermite tiene como función peso a w(x) = e−x
en el intervalo [a, b] = (−∞, +∞), por lo que las integrales a aproximar serán de la
forma Z +∞
2
e−x f (x) dx.
−∞

Los polinomios Hi (x) de la cuadratura del Gauss-Hermite son

H0 (x) = 1,
H1 (x) = 2x,
Hk+2 (x) = 2xHk+1 (x) − 2(k + 1)Hk (x), k = 0, 1, . . . , n − 2.

Los polinomios anteriores son ortogonales respecto del producto escalar. A partir de
las raíces xi de Hn (x), obtenemos los coeficientes
√
2n−1 n! π
ci = 2 2 .
n Hn−1 (xi )

Ejemplo 8. Calcular la integral

Z +∞
2
I= e−4x |x| dx
−∞

utilizando la cuadratura de Gauss-Hermite para n = 2 y n = 4.

Hagamos un cambio de variable para que el integrando tenga la forma sobre la
que podemos aplicar la cuadratura de Gauss-Hermite. Si tomamos y = 2x →
dy = 2 dx, entonces
Z +∞ Z +∞
−4x2 1 2
I= e |x| dx = e−y |y| dy.
−∞ 4 −∞

Métodos Numéricos Aplicados I

27
Tema 5. Ideas clave
Calculamos los polinomios de Hermite hasta n = 4.

H0 (x) = 1,
H1 (x) = 2x,
H2 (x) = 4x2 − 2,
H3 (x) = 8x3 − 12x,
H4 (x) = 16x4 − 48x2 + 12.

Calculamos las raíces de H2 (x) y H4 (x) ejecutando en Matlab

>> raicesH2 = roots ([4 0 -2])

raicesH2 =
0 .707106781186548
-0 .707106781186547
>> raicesH4 = roots ([16 0 -48 0 12])
raicesH4 =
-1 .650680123885785
1 .650680123885786
-0 .524647623275290
0 .524647623275290

Obtenemos los coeficientes ejecutando en Matlab

>> syms x
>> H1 =2* x;
>> c2 =2^(2 -1) * factorial (2) * sqrt (pi)/ ...
2^2 ./( double ( subs (H1 ,x, raicesH2 )). ^2)
c2 =
0 .886226925452758
0 .886226925452758
>> H3 =8* x ^3 -12* x;

Métodos Numéricos Aplicados I

28
Tema 5. Ideas clave
 >> c4 =2^(4 -1) * factorial (4) * sqrt (pi) ...
/4^2 ./( double ( subs (H3 ,x, raicesH4 )). ^2)
c4 =
0 .081312835447245
0 .081312835447245
0 .804914090005513
0 .804914090005513

Por tanto, el valor de la integral para n = 2 y n = 4 será

>> I2 =( c2 (1) * abs ( raicesH2 (1) )+ ...

c2 (2) * abs ( raicesH2 (2) ))/4
I2 =
0 .313328534328875
>> I4 =( c4 (1) * abs ( raicesH4 (1) )+ ...
c4 (2) * abs ( raicesH4 (2) )+ ...
c4 (3) * abs ( raicesH4 (3) )+ ...
c4 (4) * abs ( raicesH4 (4) ))/4
I4 =
0 .278258872775874

� Accede al vídeo: Interpretación geométrica de la integración numérica

5.6 Integración múltiple

La integración múltiple consiste en integrar sobre más de una variable. Es decir, que
f : Rn → R. En este apartado nos centraremos en el caso de f : R2 → R en el
dominio (x, y) = [a, b] × [c, d], pero es extrapolable a otras dimensiones. La Figura 6

Métodos Numéricos Aplicados I

29
Tema 5. Ideas clave
muestra un ejemplo de función con estas características.

1
0.9

0.8
0.5

0.7
0
2
0.5 0.6
0
0
y -2 -0.5
x

Figura 6: Ejemplo de f : R2 → R en recinto rectangular

Nodos equiespaciados

La técnica que se va a utilizar en esta sección es el uso del método de Trapecios para
la aproximación de una integral doble. Se podría utilizar perfectamente otro método.

Podemos escribir la integral doble como

ZZ Z b Z d 
I= f (x, y) dA = f (x, y) dy dx.
R a c

Si nombramos Z d
g(x) = f (x, y) dy, (14)
c

entonces Z b
I= g(x) dx.
a

Apliquemos ahora el método de trapecios compuesto (8). Nombrando

b−a d−c
h= , k= ,
n m

Métodos Numéricos Aplicados I

30
Tema 5. Ideas clave
podemos decir que

Z b n−1
hX
I= g(x) dx ≈ g(xi ) + g(xi+1 ),
a 2 i=0

pero como g(x) viene definida en (14),

n−1 Z d Z d 
hX
I≈ f (xi , y) dy + f (xi+1 , y) dy . (15)
2 i=0 c c

Desarrollemos la primera integral fijando xi y utilizando trapecios sobre la variable y.

Z d m−1
kX
f (xi , y) dy ≈ f (xi , yj ) + f (xi , yj+1 ).
c 2 j=0

Procediendo de forma similar sobre la segunda integral y reemplazando en (15), ob-

tenemos

n−1
(m−1 )
hk X X
I≈ f (xi , yj ) + f (xi , yj+1 ) + f (xi+1 , yj ) + f (xi+1 , yj+1 ) . (16)
4 i=0 j=0

Supongamos que n = m = 2. Entonces (16) pasa a ser

hk
Pn−1
I ≈ 4 i=0 f (xi , y0 ) + 2f (xi , y1 ) + f (xi , y2 ) + f (xi+1 , y0 ) + 2f (xi+1 , y1 ) + f (xi+1 , y2 ) =
hk
= 4
[f (x0 , y0 ) + 2f (x0 , y1 ) + f (x0 , y2 )+
+2f (x1 , y0 ) + 4f (x1 , y1 ) + 2f (x1 , y2 )+
+ f (x2 , y0 ) + 2f (x2 , y1 ) + f (x2 , y2 )] .

La Figura 7 muestra la distribución de los nodos.

𝑦2 1 2 1

𝑦1 2 4 2

𝑦0 1 2 1

𝑥0 𝑥1 𝑥2

Figura 7: Distribución de los nodos para n = m = 2

Métodos Numéricos Aplicados I

31
Tema 5. Ideas clave
Para cualesquiera valores de n y m, la distribución de nodos es la que se representa
en la Figura 8.

𝑦𝑚 1 2 2 2 1

𝑦𝑚−1 2 4 4 4 2

𝑦2 2 4 4 4 2

2 4 4 4 2
𝑦1

1 2 2 2 1
𝑦0
𝑥0 𝑥1 𝑥2 𝑥𝑛−1 𝑥𝑛
Figura 8: Distribución de los nodos para cualesquiera n y m

Nodos no equiespaciados

Cuando trabajamos con nodos no equiespaciados, mencionamos la diferentes cuadra-

turas de Gauss. De nuevo, deberemos transformar el recinto de integración [a, b] ×
[c, d] en el que corresponda a la cuadratura.

En el Ejemplo 9 vemos una adaptación de la cuadratura de Gauss-Legendre para el

caso de integrales dobles.

Ejemplo 9. Calcular la integral

Z 2 Z 1.5
ln(x + 2y) dy dx,
1.4 1

utilizando la cuadratura de Gauss-Legendre con n = m = 3.

Como 1.4 ≤ x ≤ 2 y 1 ≤ y ≤ 1.5, aplicamos los cambios de variables para que

Métodos Numéricos Aplicados I

32
Tema 5. Ideas clave
tengan el recinto [−1, 1] × [−1, 1], de modo que

2x − 3.4 2 10 2y − 2.5 2
u= → du = dx = dx, v= → dv = dy = 4 dy,
0.6 0.6 3 0.5 0.5

es decir,

x = 0.3u + 1.7 → dx = 0.3 du, y = 0.25v + 1.25 → dy = 0.25 dv.

Por tanto,
Z 2 Z 1.5 Z 1 Z 1
ln(x + 2y) dy dx = 0.075 ln(0.3u + 0.5v + 4.2) dv du.
1.4 1 −1 −1

Con los valores de la Tabla 4, los coeficientes y los nodos de la cuadratura de

Gauss-Legendre son

u1 = v1 = 0, c1 = 0.888889;
u2 = v2 = −0.774597, c2 = 0.555556;
u3 = v3 = 0.774597, c3 = 0.555556.

De este modo,
Z 1 Z 1
0.075 ln(0.3u + 0.5v + 4.2) dv du ≈
−1 −1

3 X
X 3
≈ 0.075 ci cj ln(0.3ui + 0.5vj + 4.2).
i=1 j=1

Resolvemos este sumatorio en Matlab como

>> c =[0 .888889 0 .555556 0 .555556 ];

>> C=c '*c
C =
0 .7901 0 .4938 0 .4938
0 .4938 0 .3086 0 .3086
0 .4938 0 .3086 0 .3086

Métodos Numéricos Aplicados I

33
Tema 5. Ideas clave
>> u =[0 -0 .774597 0 .774597 ] ';
>> U= repmat (u ,1 ,3)
U =
0 0 0
-0 .7746 -0 .7746 -0 .7746
0 .7746 0 .7746 0 .7746
>> v =[0 -0 .774597 0 .774597 ];
>> V= repmat (v ,3 ,1)
V =
0 -0 .7746 0 .7746
0 -0 .7746 0 .7746
0 -0 .7746 0 .7746
>> I= sum ( sum (0 .075 *C.* log (0 .3*U+0 .5*V+4 .2)))
I =
0 .4296

Así pues,
Z 2 Z 1.5
I= ln(x + 2y) dy dx = 0.429554959579526.
1.4 1

Ejemplo 10. Calcular la integral

s 2  2
ZZ
∂f ∂f
I= (x, y) + (x, y) dA,
R ∂x ∂y

donde f (x, y) es la superficie de la semiesfera

x2 + y 2 + f 2 (x, y) = 9,

y el recinto de integración es

R = {(x, y) ∈ R2 : 0 ≤ x ≤ 1, 0 ≤ y ≤ 1}.

Métodos Numéricos Aplicados I

34
Tema 5. Ideas clave
Utiliza los métodos de Simpson con n = m = 8 y de Gauss-Legendre con n = 4.
La función f que aparece en el integrando es

p
f (x, y) = 9 − x2 − y 2 ,

y sus derivadas parciales

∂f x ∂f y
(x, y) = − p , (x, y) = − p ,
∂x 9 − x2 − y 2 ∂y 9 − x2 − y 2

de modo que el integrando es

s 2  2 s
∂f ∂f x2 + y 2
(x, y) + (x, y) = .
∂x ∂y 9 − x2 − y 2

El recinto de integración es un cuadrado [0, 1] × [0, 1], por lo que la integral a

resolver es s
1 1
x2 + y 2
Z Z
I= dy dx.
0 0 9 − x2 − y 2

Procediendo de forma análoga al apartado de «Nodos equiespaciados» de la Sec-

ción 5.6 pero con el método de Simpson, la distribución de nodos es la expuesta
en la siguiente tabla para n = m = 8.

1 4 2 4 2 4 2 4 1
4 16 8 16 8 16 8 16 4
2 8 4 8 4 8 4 8 2
4 16 8 16 8 16 8 16 4
2 8 4 8 4 8 4 8 2
4 16 8 16 8 16 8 16 4
2 8 4 8 4 8 4 8 2
4 16 8 16 8 16 8 16 4
1 4 2 4 2 4 2 4 1

Métodos Numéricos Aplicados I

35
Tema 5. Ideas clave
El resultado es
ISimpson = 0.267814255559730.

Para aplicar el método de Gauss-Legendre, aplicamos los cambios de variables

u = 2x − 1 → du = 2 dx, v = 2y − 1 → dv = 2 dy.

Por tanto,
s r
1 1 1 1
x2 + y 2 u2 + 2u + v 2 + 2v + 2
Z Z Z Z
1
I= 2 2
dy dx = dv du
0 0 9−x −y 4 −1 −1 34 − u2 − 2u − v 2 − 2v

Con los valores de la Tabla 4, los coeficientes y los nodos de la cuadratura de

Gauss-Legendre para n = 4 son

u1 = v1 = −0.339981, c1 = 0.652145;
u2 = v2 = −0.861136, c2 = 0.347855;
u3 = v3 = −u1 , c3 = c1 ;
u4 = v4 = −u2 , c4 = c2 .

De este modo,

1 1
r
u2 + 2u + v 2 + 2v + 2
Z Z
1
4
dv du ≈
−1 −1 34 − u2 − 2u − v 2 − 2v

s
4 X
4
1
X u2i + 2ui + vj2 + 2vj + 2
≈ ci cj .
4
i=1 j=1
34 − u2i − 2ui − vj2 − 2vj

Procediendo de forma análoga al Ejemplo 9,

IGauss-Legendre = 0.267770529696778.

Métodos Numéricos Aplicados I

36
Tema 5. Ideas clave