tema

31 MATEMÁTICAS

Integración numérica.
Métodos y aplicaciones.
24-13823-13

Temario 1993
 tema 31

matemáticas

1. Resumen de algunos resultados sobre operadores

2. Diferencias de orden negativo. Operadores ∆–1 y δ–1

2.1. Operador ∆–1

2.2. Operador δ–1

3. Fórmula de integración numérica de Newton con diferencias

laterales descendentes

4. Fórmula de integración numérica con diferencias laterales

ascendentes

5. Fórmula mixta de Gregory para diferencias laterales

6. Fórmula de Gauss-Encke para integración con diferencias

centrales

7. Cálculo de la integral de una función entre dos límites que no

figuran en la tabla. Fórmula de Newton-Cotes

8. Otras fórmulas de integración

9. El error cometido en la integración numérica

3
 tema 31

matemáticas

INTRODUCCIÓN

Con el desarrollo del Cálculo Infinitesimal en los últimos siglos, ha ido creciendo el gus-
to por la perfección del razonamiento, en detrimento del desarrollo de los procedimientos
constructivos utilizables en la práctica. Hasta tal punto ha producido esto un cambio en la
situación, que el estudio de los métodos numéricos lleva camino de convertirse en una de las
ramas más importantes de la Matemática y, desde luego, en una de las de más actualidad.
Un método numérico de resolución de un determinado problema, es un conjunto de reglas
que permite la obtención, mediante un número finito de operaciones elementales, de un re-
sultado que se aproxima de alguna manera a la solución de aquel problema. Debido a ello,
los métodos numéricos lo único que suelen proporcionar es soluciones tan próximas como
se quiera (en teoría, puesto que al efectuar los cálculos se cometen casi siempre pequeños
errores debido a que las calculadoras tienen una precisión, aunque grande, limitada) a la
solución exacta, pero en general, no se obtiene ésta. El conjunto de reglas que definen el
método suele indicarse con el nombre de algoritmo.
El Cálculo Numérico ha hecho que problemas resueltos vuelvan a estudiarse por otros
métodos, por no ser aquéllos utilizables más que en demostraciones y razonamientos teó-
ricos.
En los otros temas, hemos estudiado varios métodos matemáticos para calcular la integral
definida de determinadas funciones. Pero en múltiples problemas de Física, Química, Bio-
logía, Ingeniería, y un largo etcétera, interesa a menudo integrar funciones de las que no
se conoce su expresión analítica; simplemente se nos da su tabla de valores, o una curva
representativa, obtenidas ambas experimentalmente.
Hay ocasiones en las que no se sabe, ni siquiera, calcular la primitiva de alguna función, y
otras veces, no interesa obtener con total exactitud dicha primitiva.
El campo de la integración numérica aborda todos estos problemas, y su objetivo estriba en
obtener la integral definida en los casos mencionados.
La integración numérica no consigue resultados exactos, pero sí aproximaciones tan bue-
nas como queramos. Simples fórmulas como la del trapecio, muy fáciles de programar
incluso con calculadoras de bolsillo, resuelven con bastante rapidez y eficacia integrales de
aspecto complicadísimo.
Acabamos esta introducción diciendo que los algoritmos empleados en la integración nu-
mérica, se llaman fórmulas de cuadratura.

5
tema 31

matemáticas

1 Resumen de algunos resultados sobre operadores

Algunos resultados, parte de los cuales ya se estudiaron en el tema de «Interpola-
ción», son los siguientes:
„„ Operador E (de avance)
Eyx = yx+1
Es tal que: Epyx = yx+p, E–pyx = yx–p y se cumple que Ep. E–p = E0 = 1
„„ Operador ∆ (diferencia lateral descendente)
∆ yx = yx+1 – yx, y es tal que ∆2yx = ∆ yx+1 – ∆ yx, y así sucesivamente, es decir:
∆p yx = ∆p–1 yx+1 – ∆p–1 yx
Es conocido que: E = ∆ + 1
„„ Operador ∇ (diferencia lateral ascendente)
∇ yx = yx – yx–1 ⇒ ∇p yx = ∇p–1 yx – ∇p–1 yx–1
„„ Operador D (derivada)
D = LE = ⇒ E = eD, y también D = L(1 + ∆), donde L significa logaritmo ne-
periano.
„„ Operador δ (diferencia central)
δ yx = y 1 −y 1 ⇒ δ p y x = δ p −1 y 1 − δ p −1 y 1 , ...
x+ x− x+ x−
2 2 2 2

1 1
−
y se sabe que: δ = E 2 − E 2 .

6
 tema 31

matemáticas

2 Diferencias de orden negativo. Operadores ∆–1 y δ–1

2.1. Operador ∆–1

Partiendo de una tabla de una función yx:

x y
. .
. .
. .
–1 y–1
0 y0

1 y1
. .
. .
. .
hemos definido, en el tema ya citado, las diferencias p-ésimas ∆p yn, siendo p un
entero positivo por las fórmulas de recurrencia
∆ yn = yn+1 – yn,   ∆p yn = ∆p–1 yn+1 – ∆p–1 yn
Conviene remarcar que la primera de estas fórmulas viene de un caso particular de
la segunda si ponemos ∆0 yn = yn, teniendo en cuenta que ∆0 = 1.
Vamos a definir el operador ∆–1 por una condición idéntica a las precedentes, este
operador deberá verificar la condición yn = ∆0 yn = ∆–1 yn+1 – ∆–1 yn. ∆–1 está defini-
do salvo una constante, porque sólo las diferencias de dos valores consecutivos de
∆–1 son conocidas. Esta constante es indiferente en las aplicaciones, pues ∆–1 no
intervendrá más que de forma tal, que dicha constante se elimine.
La expresión de ∆–1 yn es fácil de formar por recurrencia.
Si n es positivo, se encuentra: ∆–1 yn = ∆–1 y0 + y0 + y1 + ... + yn–1;
y, por otro lado, si n es negativo: ∆–1 y–n = A– 1 y0 – y–1 – y–2 – ... – y–n.
Siendo n > p, estos enteros pueden tomar signos cualesquiera y considerando las
tres posibilidades para esos signos, se verifica que:
∆–1 yn – ∆–1 yp = yp + yp+1 + ... + yn–1, n > p
Vamos a mostrar que la notación ∆–1 está justificada por el hecho de que ∆–1 es el
inverso de ∆ a derecha e izquierda, es decir que
∆ · ∆–1 = 1   y   ∆–1 · ∆ = 1
La verificación es inmediata: por ejemplo si n > 0
∆ · ∆–1 yn = ∆–1 yn+1 – ∆– 1 yn = (∆–1 y0 + y0 + ... + yn) – (∆–1 y0 + y0 + ... + yn – 1) = yn

7
tema 31

matemáticas

Recíprocamente, el operador considerado es el único inverso de ∆. Busquemos

por ejemplo un operador A, tal que A∆ = 1, y apliquemos los dos miembros de
esta identidad a yn:
A∆ yn = yn, es decir, Ayn+1 – Ayn = yn
Esta es la condición que define a ∆–1.
Hemos visto anteriormente que el operador ∆–1 está definido, salvo una cons-
tante aditiva. En adelante para fijar las ideas, supondremos que ∆–1 y0 = 0.
De tal forma:
∆–1 yn = y0 + y1 + ... + yn–1,    si n > 0
∆–1 y–p = –y–1 – y–2 – ... – y–p,    si p > 0
Si bien las fórmulas que vamos a establecer son válidas cualquiera que sea el
valor adoptado para ∆–1 y0.
La disposición que se utiliza para trabajar con el operador ∆–1 es la siguiente:
x ∆–1 y
. . .
. . .
–2 ∆–1 y–2 y–2
–1 ∆–1 y–1 y–1

0 ∆–1 y0 y0
+1 ∆–1 y1 y1
+2 ∆–1 y2 y2

. . .
. . .
. . .
La formación de cada columna es sencilla, se pone ∆–1 y0 = 0, y se forma cada
∆–1 yi, de índice positivo, por adición del ∆–1 colocado inmediatamente por enci-
ma, y de la yi colocada inmediatamente a la derecha de ∆–1 yi; para las ∆–1 yi de
índice i negativo el proceso es muy similar, si bien se disminuye del ∆–1 colo-
cado inmediatamente por debajo, la yi colocada inmediatamente a la derecha de
∆–1 yi+1.
Si 0 es el principio de la tabla y n el final, la columna de las ∆–1 yi comenzará en
∆–1 y0 = 0 y acabará en ∆–1 yn+1.
De idéntica manera se van formando las ∆–2 yi, después las ∆–3 yi, y así sucesiva-
mente hasta las ∆–k yi.
Poniendo ∆–p y0 = 0, ∆–1 puede expresarse como suma finita de y, pues:
∆–1 yn – ∆–1 yp = yp + yp–1 + ... + yn–1,   n > p

8
 tema 31

matemáticas

Se verifica que este operador de suma finita, llamado ∆–1, satisface las reglas clá-
sicas del cálculo algebraico. Basta para ello tener presente que:
1
∆–1 es igual a
∆

2.2. Operador δ–1

Se le define como el precedente por la condición:

yi = δ−1 y 1 − δ−1 y 1
i+ i−
2 2
Tomando
δ−1 y 1 =0
−
2

las δ–1 son lo mismo que los ∆–1, y se les dispone en los mismos lugares, siendo
sólo diferentes por la notación.
. . . . .
. . . . .
. y−1 . δ y−1 .
2

δ−1 y 1 δy 1 δ3 y 1
− − −
2 2 2

y0 δ2 y0
δ−1 y 1 δy 1 δ3 y 1
2 2 2

y1 δ y12

−1
δ y3 δy 3 δ3 y 3
2 2 2

y2 δ y22

−1
δ y5 δy 5 δ3 y 5
2 2 2

y3 δ y32

y análogamente se definen las δ–p y.

Las diferencias de orden negativo que acabamos de ver tienen la ventaja de que
con ellas podemos efectuar desarrollos de cualquier operador.

9
tema 31

matemáticas

3 Fórmula de integración numérica de Newton con

diferencias laterales descendentes
El problema que se nos plantea es el siguiente: Tenemos una tabla de valores de
una función:
x   x0   x1   x2 ......  xn
y    y0   y1   y2 ......  yn

∫
b
y pretendemos calcular una integral del tipo y x dx, en la cual los límites de
a
integración supondremos que pertenecen a la tabla. Supondremos asimismo que
la tabla es de intervalo constante.

XX El operador D–1 de integración

Está definido:

∫ ∫
n p
D −1 yn = y x dx ⇒ y x dx = D −1 y p − D −1 yn
0 n

El problema que tenemos planteado es desarrollar D–1 según las potencias de ∆.

Recordemos que D = L(1 + ∆), y consideremos el operador P = L–1 (1 + ∆) que
verifica la ecuación PD = 1.
Representa, pues, una primitiva de la función yx:
D–1 yn = Pyn + Cl

∫
n
Cl es una constante que como hemos tomado por convenio D −1 yn = y x dx,
entonces D–1 y0 = 0, por lo tanto C1 = –Py0 0

Por otro lado recordemos que:

1 1 1
L(1 + ∆) = ∆ – ∆2 + ∆3 – ∆4 + ...
2 3 4
P = L–1 (1 + ∆) = H –1
–1
∆–1 + H 0–1 ∆0 + H 1–1 ∆1 + H 2–1 ∆2 + ...
i
donde los H–1 son unos coeficientes que tenemos que determinar.
Hemos de obligar a que PD sea igual a 1, es decir, a que:
L–1(1 + ∆) · L(1 + ∆) = 1
Multiplicando miembro a miembro los dos desarrollos e identificando coeficien-
tes, se obtiene:
Coeficiente de ∆0: H –1
–1
=1
1 –1 1
Coeficiente de ∆: H 0–1 – H –1 = 0 ⇒ H 0–1 =
2 2
1 1 1
Coeficiente de ∆2: H 1–1 – H 0–1 + H –1 –1
= 0 ⇒ H 1–1 = –
2 3 12

10
 tema 31

matemáticas

Análogamente se calcularán:
1 19 3 863 275
H 2–1 = , H 3–1 = – , H 4–1 = , H 5–1 = – , H 6–1 = ,
24 720 160 60480 24192
3393
H 7–1 – , etc.
3628800
Pyn = H –1
–1
∆–1 yn + H 0–1 yn + H 1–1 ∆ yn + H 2–1 ∆2 yn + ...

Py0 = H –1
–1
∆–1 y0 + H 0–1 y0 + H 1–1 ∆ y0 + H 2–1 ∆2 y0 + ... = –C1

∫
n
luego: y x dx = Pyn – Py0 =
0

= [H –1
–1
∆–1 yn + H 0–1 yn + H 1–1 ∆ yn + H 2–1 ∆2 yn + ...]
– [H –1
–1
∆–1 y0 + H 0–1 y0 + H 1–1 ∆ y0 + H 2–1 ∆2 y0 + ...] =
= H –1
–1
(∆–1 yn – ∆–1 y0) + H 0–1 (yn – y0) + H 1–1 (∆ yn – ∆ y0) +
+ H 2–1 (∆2 yn – ∆2 y0) + ...
  ∆–1 y0
  .   y0
  .   . ∆ y0
  .   . . ∆2 y0 ......
  ∆–1 yn
    yn
  .   . ∆ yn
  .   . . ∆2 yn ......
Vamos a prestar nuestra atención al término H –1
–1
(∆–1 yn – ∆–1 y0).
Recordemos que:
∆h yn = ∆h–1yn+1 – ∆h–1yn
∆0 yn = yn = ∆–1 yn+1 – ∆–1 yn
∆0 yn–1 = yn–1 = ∆–1 yn – ∆–1 yn–1
luego,
∆–1 yn = ∆–1 yn–1 + yn–1 = ∆–1 yn–2 + yn–2 + yn–1
y así sucesivamente llegaríamos a:
∆–1 yn = ∆–1 y0 + y0 + y1 + ...... + yn–1 ⇒ ∆–1 yn – ∆–1 y0 = y0 + y1 + ...... + yn–1
con lo cual,
H –1
–1
(∆–1 yn – ∆–1 y0) = y0 + y1 + ...... + yn–1

1 1
El siguiente sería (yn – y0), el siguiente – (∆yn – ∆y0), y así sucesivamente...
2 12

11
tema 31

matemáticas

Sustituyendo nos queda:

∫
n 1 1
y x dx = y0 + y1 +  + yn −1 + (∆yn − ∆y0 ) +
0 2 2
1 19 3
+ ( ∆ 2 y n − ∆ 2 y0 ) − (∆ 3 yn − ∆ 3 yn ) + ( ∆ 4 y n − ∆ 4 y0 ) − 
24 720 160

llamada fórmula de Newton.

Dicha fórmula debe aplicarse cuando los límites de integración están cerca del
principio de la tabla.

∫ ∫
b n
Si nos piden calcular , calcularemos y la multiplicaremos por el número de
intervalos.
a 0

12
 tema 31

matemáticas

4 Fórmula de integración numérica con diferencias

laterales ascendentes
Su deducción es totalmente similar a la anterior.
Escribiendo 1 + ∆ en la forma:
−1
 ∆
1 + ∆ = 1 − 
 E

y como P = L–1 (1 + ∆), entonces, se obtiene:

∫
n
y x dx = Pyn − Py0 = H −−11 ( ∆ −1 yn+1 − ∆ −1 y1 ) − H −01 ( yn − y0 ) + H −11 ( ∆yn−1 − ∆yy−1 ) −
0
1 1
− H −21 ( ∆ 2 yn−2 − ∆ 2 y−2 ) +  = y0 + y1 +  + yn−1 + yn + H −11 ( ∆yn−1 − ∆y−1 ) −
2 2
− H −21 ( ∆ 2 yn−2 − ∆ 2 y−2 ) + 

Este método se aplica cuando los límites de integración estén al final de la tabla.
∆2 y–2
∆ y–1 .
y 0 . .
∆–1 y1 . . .
. . . ∆2 yn–2
. . ∆ yn–1
. yn
∆–1 yn+1

13
tema 31

matemáticas

5 Fórmula mixta de Gregory para diferencias

laterales
En este caso se usan las diferencias descendentes para y0, y ascendentes para yn,
y siguiendo un proceso similar al del apartado 4, se obtiene lo siguiente:

∫
n
1 1 1
y x dx = y0 + y1 +  + yn−1 + yn − ( ∆yn−1 − ∆y0 ) −
0 2 2 12
1 19 3
− ( ∆ y n− 2 + ∆ y0 ) −
2 2
( ∆ y n−3 − ∆ y0 ) −
3 3
( ∆ 4 y n− 4 + ∆ 4 y0 ) − 
24 720 160

que nos proporciona la integral por medio de las diferencias descendentes de y0,
y ascendentes de yn:
∆–1 y0
. y0
. . ∆ y0
. . . ∆2 y0 ...
. . . ∆2 yn–2 ...
. . ∆ yn–1
. yn
∆–1 yn+1
Remarcamos el hecho de que el primer término se escribe:
∆–1 yn+1 – ∆–1 y0 = y0 + y1 + ... + yn.
Esta fórmula que acabamos de ver, que es debida a Gregory, sólo es aplicable si 0
es el principio de la tabla y n el extremo, es decir, para calcular la integral definida
de yx para todo el intervalo de la tabla.

14
 tema 31

matemáticas

6 Fórmula de Gauss-Encke para integración con

diferencias centrales

Al igual que sucede con el problema de la interpolación en una tabla, hay aquí
también dos fórmulas de integración utilizando las diferencias centrales. La pri-
mera de ellas nos da la integral entre dos límites que figuran en la tabla, por ejem-

∫
1
plo y x dx.
0

La segunda, por el contrario, nos da la integral entre dos límites iguales a un ente-
1

∫
1
ro más , por ejemplo 2
1
y x dx.
2 −
2

Para obtener la primera, debemos desarrollar:

 δ2 δ2 
L−1 1 + + δ 1 + 
P  2 4 
=
µ δ2
1+
4

según las potencias de δ.

Recordemos que µ es el operador media aritmética que vale

y 1 +y 1 
µy x =  
x+ x−
2 2
2 
 
 

El desarrollo se simplifica teniendo en cuenta que:

2
δ2 δ2  δ δ2 
1+ + δ 1+ =  + 1+ 
2 4  2 4 

Este desarrollo no contiene más que potencias impares de δ, como puede com-
probarse cambiando δ por –δ en el argumento del logaritmo y teniendo en cuenta
que:

 δ2 δ2   δ δ2 
2 + δ + − + + =1
4   2 4 
1 1


P
= N −−11 δ−1 + N −11 δ + N −31 δ3 + 
µ

15
tema 31

matemáticas

Se obtienen los valores numéricos de los primeros coeficientes por los métodos
clásicos de los desarrollos limitados, y se obtiene:

1 11 191 2497
N −−11 = 1, N −11 = − , N −31 = , N −51 = − , N −71 =
12 720 60480 3628800

∫
n
y así: y x dx = Pyn + C1; con C1 = –Py0
0

 δ−1 y 1 + δ−1 y 1 δ−1 y 1 + δ−1 y 1 

2
y x dx = N −−11 
n+ n− −

∫
n
− +
2 2 2
0 2 2
 

 δy 1 + δy 1 δy 1 + δy 1 
2
+N 
n+ n− −
1 2 2
− 2 +
 −1
2 2 
 

Esta fórmula se llama de Gauss-Encke.

Si n es entero, las diferencias que figuran en el segundo miembro son las dife-
1
rencias de orden impar para las abscisas enteras más , y como consecuencia se
2
encuentran en la tabla.

∫
n 1 1
La integral y x dx se obtiene de la fórmula anterior cambiando por p +
p 2 2
1 1
y – por p – .
2 2
La fórmula de Gauss-Encke puede ser puesta bajo una forma más simple teniendo
en cuenta que:

δ−1 y 1 = y0 + y1 +  + yn , δ−1 y 1 = y0
n+
2 2

δ−1 y 1 = y0 + y1 +  + yn−1, δ−1 y 1 =0

n− −
2 2

el primer corchete se escribe pues:

1 1
y0 + y1 + y2 + ... + yn–1 + yn
2 2
Análogamente:
δy 1 = yn+1 − yn , δy 1 = yn − yn−1
n+ n−
2 2

δy 1 = y1 − y0 , δy 1 = y0 − y−1
−
2 2

1
y para el segundo corchete: [ yn+1 – yn–1 – y1 + y–1]
2

16
 tema 31

matemáticas

Todos los términos se transforman de la misma manera y se obtiene finalmente:

1 1  1
∫
n
y x dx = N −−11  y0 + y1 +  + yn−1 + yn  + N −11 ( yn+1 − yn−1 −  − y1 − y−1 ) +
0 2 2  2

1 3 2 1
+ N −1 ( δ yn+1 − δ2 yn−1 − δ2 y1 + δ2 y−1 ) + N −51 ( δ4 yn+1 − δ4 yn−1 − δ4 y1 + δ4 y−1 ) + 
2 2

que hace intervenir las diferencias de orden par para las abscisas de la tabla de
orden menos elevado que los de la primera fórmula.
Para n = 1, la fórmula de Gauss-Encke se simplifica:

1 −1  −1  1  
∫
1
y x dx = N −1  δ y 3 − δ−1 y 1  + N −11  δy 3 − δy 1  +
0 2  2
−
2
2  2 −
2

1 3  3 
+ N −1  δ y 3 − δ3 y 1  + 
2  2
− 
2

1
Se recuerda que el primer término es igual a (y0 + y1). Esta fórmula utiliza las
diferencias: 2

x
−1 δ−1 y 1 . δy 1 . δ3 y 1 ...
− − −
2 2 2

0 .
+1 δ−1 y 3 δy 3 δ3 y 3 ...
2 2 2

+2

Bajo la segunda forma, esta fórmula se escribe:

∫
1
1 1 1
y x dx = N −1 ( y0 + y1 ) + N −11 ( y2 − y0 − y1 + y−1 ) +
0 2 2
1 3 2
+ N −1 (δ y2 − δ2 y0 − δ2 y1 + δ2 y−1 ) + 
2

Para obtener la otra fórmula, desarrollaremos P siguiendo las potencias de δ, y las

mismas observaciones hechas nos muestran que
P = N –1–1 δ–1 + N 1–1 δ + N –1
3
δ3 + ...

1 17 367 27.859
con N –1–1 =1, N 1–1 = 3
, N –1 =– , N 5–1 = , N 7–1 = – .
24 5.760 967.680 464.486.400

∫
n
y x dx = Pyn + C1, pero si tomamos:
0
Pyn = N –1–1 δ–1 yn + N 1–1 δ yn + ...

17
tema 31

matemáticas

C1 = –Py0 = –N –1–1 δ–1 y0 – N 1–1 δ y0 – ..., se obtiene una expresión que depende de
las diferencias δ–1 y0, δy0, ..., δ–1 yn, δ yn, que no figuran en la tabla.
Se utiliza, sin embargo, esta fórmula para el cálculo de

∫
p
y x dx = Py p − Py− p
−p

1
y es necesario tener para p un entero más :
2
 
1
n+

∫ y x dx = N −11  δ−1 y 1 − δ−1 y 1  +

2

 2
1 + − −
− n− n n
2 2

   
+ N −11  δy 1 − δy 1  + N −31  δ3 y 1 − δ3 y 1  + 
 n+ 2 − n− 
2  n+ 2 − n− 
2

Se recuerda que el primer paréntesis es igual a:

y–n + y–(n–1) + ... + y–1 + y0 + y1 + ... + yn
1

En particular, si n = 0:
∫ y x dx = δ−1 y 1 − δ−1 y +
2
−1 −
1
2 2 2

   
+ N −11  δy 1 − δy 1  + N −31  δ3 y 1 − δ3 y 1  + 
 2 − 
2  2
− 
2

y como δ p y 1 − δ p y 1 = δ p−1 y0 , se obtiene:

−
2 2
1

∫ y x dx = y0 + N −11δ2 y0 + N −31δ4 y0 + 
2
−1
2

que utiliza las diferencias centrales situadas sobre la línea de y0.

18
 tema 31

matemáticas

7 Cálculo de la integral de una función entre dos

límites que no figuran en la tabla. Fórmula de
Newton-Cotes
∫
n
El problema que se nos plantea es integrar y x dx = Pyx + C1, donde C1 = –Py0,
0

es decir, reemplazando yx por yx = Exy0, y Pyx por

Pyx = D–1 yx = PExy0 + C1

∫
n
y x dx = Pyx + C1 = PExy0 – Py0
0

Consideremos que Ex = e x L (1+∆) =

x x2 2 xn n
= 1+ L (1 + ∆ ) + L (1 + ∆ ) +  + L (1 + ∆ ) + 
1! 2! n!

tenemos que desarrollar PEx.

2 x3 2 xn
Como P = L–1 (1 + ∆), PEx = L–1 (1 + ∆) + x L(1 + ∆) + L (1 + ∆) + + ... +
Ln–1 (1 + ∆) + ... 2! 3! n!

Considerando: Lp (1 + ∆) = Hpp ∆p + Hpp+1 ∆p+1 + ... + Hpp+k ∆p+k + ..., tenemos:

PEx = H –1
–1
(x) ∆–1 + H 0–1 (x) + H 1–1 (x) ∆ + H 2–1 (x) ∆2 + ...
x 0
donde: H –1
–1
(x) = 1, H 0–1 (x) = H 0–1 + H 0,
1!
x 1 x2 1
H 1–1 (x) = H 1–1 + H + H , ...
1! 0 2 ! 1
y por tanto,

∫
n
y x dx = C1 + ∆–1 y0 + H 0–1 (x) y0 + H 1–1 (x) ∆y0 + H 2–1 (x) ∆2y0 + ...
0

Es decir, reemplazando C1, por su expresión en diferencias descendentes de y0:

∫
n
y x dx = [H 0–1 (x) – H 0–1] y0 + [H 1–1 (x) – H 1–1] ∆y0 + [H 2–1 (x) – H 2–1] ∆2 y0 + ...
0

llamada fórmula de Newton-Cotes.

Un cálculo análogo nos proporciona la fórmula en función de las diferencias as-
cendentes de y0:
(H –1–1)’ (x) = 1
(H 0–1)’ (x) = –H 0–1 + H 00 x
x x2
(H 1–1)’ (x) = H 1–1 – H 10 + H 11
1! 2!
x x2 x3
(H 2–1)’ (x) = –H 2–1 + H 20 – H 21 + H 22 + ...
1! 2! 3!

∫
n
y x = [H 0–1 (x) – H 0–1] y0 + [H 1–1 (x) – H 1–1] ∇y–1 + [H 2–1 (x) + H 2–1] ∇2 y–2 + ...
0

y también podríamos obtenerla, lógicamente, en función de las diferencias cen-

trales.

19
tema 31

matemáticas

8 Otras fórmulas de integración

Están basadas en las de Newton-Cotes que nos permite calcular la integral

∫
n
y x dx por medio de diferencias laterales descendentes.
0

Recordemos la fórmula de interpolación de Newton.

x x( x − 1) 2 x( x − 1) ( x − 2) 3
y x = y0 + ∆y0 + ∆ y0 + ∆ y0 + 
1! 2! 3!
entonces:

x2  x 3 x 2  ∆ 2 y0
∫
n
y x dx = xy0 + ∆y 0 +  −  +
0 2  3 2  2

 x4  ∆ 3 y0  x 5 3 x 4 11 3  ∆ 4 y0
+  − x3 + x2  + − + x + 3x 2  +
 4  6  5 2 3  24

En las aplicaciones, si el intervalo de la tabla es h y no 1 como hemos supuesto en

la fórmula precedente, escribiremos:

x2
∫
1 nh
f (u ) du = xy0 + ∆y 0 + 
h 0 2

Una sustitución análoga deberá ser hecha en las fórmulas que siguen.
„„ Fórmula del trapecio
Despreciando diferencias a partir de ∆2y0 y haciendo n = 1 obtenemos:

∫
1
1
y x dx = ( y0 + y1 )
0 2

Aplicándola sucesivamente a (0, 1), (1, 2), ..., (n – 1, n) esta fórmula nos da:

∫
n
1 1
y x dx = y0 + y1 + y2 +  + yn−1 + yn
0 2 2

„„ Fórmula de Simpson
Despreciando las diferencias a partir de ∆3y0 y adoptando x = 2:

∫
2 1
y x dx = ( y0 + 4 y1 + y2 )
0 3
Aplicándola n veces:

∫
n 1
y x dx = ( y0 + 4 y1 + 2 y2 + 4 y3 + 2 y4 + 4 y5 + )
0 3

Omitimos, llegado a este punto, otra serie de fórmulas de integración numérica

para evitar que el tema se alargue en exceso.

20
 tema 31

matemáticas

9 El error cometido en la integración numérica

El cálculo del límite superior del error cometido al adoptar para la integral de una
función el valor dado por una de las fórmulas que hemos expuesto, está basado en
el método siguiente:
Recordemos que en el tema de interpolación se comentó que había un cierto error
al aproximar una función mediante un polinomio.
R (x) = f (x) – P(x) ⇒ f(x) = P(x) + R(x)
integrando:

∫ ∫ ∫
b b b
f ( x ) dx = P( x ) dx + R( x ) dx
a a a

∫
b
ε= R( x ) dx es el error cometido en la integración numérica.
a
Recordemos también que R(x) se podrá poner:
( x − x0 ) ⋅ ( x − x1 ) ⋅  ⋅ ( x − xn ) ( n+1)
R( x ) = f (ξx )
( n + 1) !

donde ξx es un valor de la variable comprendido entre el más pequeño y el más

grande de los x0, x1, ..., xn, y dependiente de x.

∫
1 b
ε = ( x − x0 ) ⋅ ( x − x1 ) ⋅  ⋅ ( x − xn ) f ( n +1) (ξ x ) dx
(n + 1) ! a

Para transformar esta integral, recordemos el teorema de la media.

Si f(x) y ϕ(x) son dos funciones continuas, y si ϕ(x) mantiene signo constante en
el intervalo (α, β), existe un valor ϑ interior a (α, β) tal que:
β β

∫ α
f ( x ) ϕ( x ) dx = f ( ϑ ) ∫ α
ϕ( x ) dx

Para transformar la expresión de ε, descomponemos el intervalo (a, b) en subin-

tervalos (αi, αi+1) en los cuales
ϕ(x) = (x – x0) · (x – x1) · ... · (x – xn) mantiene signo constante.

f ( n+1) ( ξi )
( xi , xi +1 ) 
→ε= ∑ i
Ai
( n + 1) !
αi+1
Ai = ∫ αi
(x – x0) · (x – x1) · ... · (x – xn) dx

Los ξi están comprendidos entre el más grande y el más pequeño de los x0, x1, ..., xn.
Los Ai son fácilmente calculables. Si se sabe calcular f (n+1) (x) se puede deducir de
la fórmula precedente un límite superior de ε.
En el desarrollo del tema hemos obtenido una serie de fórmulas para el cálculo de
∫0 yx dx. En el caso de que nos pidan calcular ∫ y dx, calcularemos ∫ y dx, y después
n b n
x x
a 0

la multiplicaremos por el intervalo h.

21
tema 31

matemáticas

BIBLIOGRAFÍA

COHEN, A., et al.: Análisis Numérico. Ed. Reverté, 1977 [agotado].

DAHLQUIST y BJÖRCK: Numerical Methodos. Ed. Prentice-Hall.
DEMIDOVICH, B. P., y MARON, I. A.: Cálculo Numérico Fundamental. Ed. Paraninfo, 1993.
FERNÁNDEZ, P., et al.: Métodos de integración numérica. Universidad Politécnica de Valencia,
1997.
GASCA GONZÁLEZ, M.: Cálculo numérico II. U.N.E.D., 1991.
SCHEID, F.: Análisis numérico. McGraw-Hill, 1985.
STOER y BURLIRSCH: Introduction to Numerical Analysis. Ed. Springer Verlay.
WHITAKER AND ROBINSON: The Calculus of observations. Ed. Blackie and Sons.

22
 tema 31

matemáticas

RESUMEN

Integración numérica.
Métodos y aplicaciones.

1.
1 Resumen de algunos resultados sobre operadores
Un método numérico de resolución es un conjunto de reglas que permite la obtención de
un resultado aproximado de la solución. El conjunto de reglas son el algoritmo (fórmulas
de cuadratura). La integración numérica tiene por objetivo obtener integrales definidas en
casos en los que no se sabe o puede calcular la primitiva o no interesa un valor exacto.

2.
2 Diferencias de orden negativo. Operadores ∆–1 y δ–1

2.1. Operador ∆–1

∆–1 yn = ∆–1 y0 + y0 + y1 + ... yn
∆–1 y–n = ∆–1 y0 + y–1 + y–2 + ... y–n
∆∆–1 = ∆–1∆ = 1

2.2. Operador δ–1

Se le define como el precedente por la condición:
yi = δ −1 y 1 − δ −1 y 1
i+ i−
2 2
δδ = δ δ = 1
–1 –1

3.
3 Fórmula de integración numérica de Newton con
diferencias laterales descendentes

∫
n 1 1
y x dx = y0 + y1 +  + yn −1 + (∆yn − ∆y0 ) +
0 2 2
1 19 3
+ ( ∆ 2 y n − ∆ 2 y0 ) − (∆ 3 yn − ∆ 3 yn ) + ( ∆ 4 y n − ∆ 4 y0 ) − 
24 720 160

Se aplica cuando los límites de integración están cerca del principio de la tabla de valor
de la función.

23
tema 31

matemáticas

4.
4 Fórmula de integración numérica con diferencias
laterales ascendentes

∫
n
y x dx = Pyn − Py0 = H −−11 (∆ −1 yn +1 − ∆ −1 y1 ) − H −01 ( yn − y0 ) + H −11 (∆yn −1 − ∆yy−1 ) −
0
1 1
− H −21 (∆ 2 yn − 2 − ∆ 2 y−2 ) +  = y0 + y1 +  + yn −1 + yn + H −11 (∆yn −1 − ∆y−1 ) −
2 2
− H −21 (∆ 2 yn − 2 − ∆ 2 y−2 ) + 

Se aplica cuando los límites de integración están cerca del final de la tabla de valore de la
función.

5.
5 Fórmula mixta de Gregory para diferencias laterales

∫
n 1 1 1
y x dx = y0 + y1 +  + yn −1 + yn − (∆yn −1 − ∆y0 ) −
0 2 2 12
1 19 3
− ( ∆ 2 y n − 2 + ∆ 2 y0 ) − ( ∆ 3 y n −3 − ∆ 3 y0 ) − ( ∆ 4 y n − 4 + ∆ 4 y0 ) − 
24 720 160

Se aplica cuando los límites de integración son el principio y el final de la tabla.

6.
6 Fórmula de Gauss-Encke para integración con
diferencias centrales
1 1  1
∫
n
y x dx = N −−11  y0 + y1 +  + yn −1 + yn  + N −11 ( yn +1 − yn −1 −  − y1 − y−1 ) +
0 2 2  2
1 3 2 1
+ N −1 (δ yn +1 − δ 2 yn −1 − δ 2 y1 + δ 2 y−1 ) + N −51 (δ 4 yn +1 − δ 4 yn−1 − δ 4 y1 + δ 4 y−1 ) + 
2 2

7.
7 Cálculo de la integral de una función entre dos límites
que no figuran en la tabla. Fórmula de Newton-Cotes
∫
x
y x = [H 0 (x) – H 0 ] y + [H 1 (x) – H 1 ] ∇y + [H 2 (x) + H 2 ] ∇2 y + ...
0 –1 –1 0 –1 –1 –1 –1 –1 –2

8.
8 Otras fórmulas de integración
„„ Fórmula del trapecio

∫
n 1 1
y x dx = y0 + y1 + y2 +  + yn −1 + yn
0 2 2

„„ Fórmula de Simpson

∫
n 1
y x dx = ( y0 + 4 y1 + 2 y2 + 4 y3 + 2 y4 + 4 y5 + )
0 3

24
 tema 31

matemáticas

9.
9 El error cometido en la integración numérica

∫
b
ε = R( x) dx tal que R (x) = f (x) – P(x) donde P(x) es el polinomio interpolador.
a

∫
1 b
ε = ( x − x0 ) ⋅ ( x − x1 ) ⋅  ⋅ ( x − xn ) f ( n +1) (ξ x ) dx con ξx [x0, xn] o equivalentemente.
(n + 1) ! a

f ( n +1) (ξ i )
∑A
α i +1
→ε =
( xi , xi +1 ) 
i
i
(n + 1) !
tal que Ai = ∫αi
(x – x0) · (x – x1) · ... · (x – xn) dx

25