Metodo Simpson 711v1u

Ing Yamil Armando Cerquera Rojas

[email protected]

INTEGRACION NUMERICA Método de Simpson Ing Yamil Armando Cerquera Rojas – [email protected] Especialista en Sistemas Universidad Nacional Docente Universidad Surcolombiana Neiva – Huila

Contenido DEFINICIÒN ............................................................................................................................ 2 INTRODUCCIÒN ..................................................................................................................... 2 OBJETIVOS............................................................................................................................. 3 GENERAL............................................................................................................................ 3 ESPECÍFICOS ...................................................................................................................... 3 OBSERVACIONES PRELIMINARES ......................................................................................... 3 CÁLCULO DE ÁREAS.............................................................................................................. 5 EL MÉTODO DE SIMPSON ................................................................................................. 5 DESARROLLO DEL MODELO DE SIMPSON:...................................................................... 6 Deducción del modelo a partir de la Ecuación de la parábola........................................ 6 Deducción del modelo a partir del polinomio de lagrange............................................. 8 Generalización del modelo para n subintervalos.......................................................... 11 Ejemplos.......................................................................................................................... 13 Programación del método de Simpson en lenguaje C. ........................................... 20 La jerarquía de clases para C...................................................................................... 21

Universidad Surcolombiana – Neiva – Huila - Colombia

1 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

DEFINICIÒN De acuerdo a la definición del diccionario, integrar significa “unir todas las partes en un todo; unificar, indicar la cantidad total,…”.1 . Matemáticamente, la b

integración se representa por

∫ f ( x)dx . En los primeros años de ingeniería, se ven a

apartes de cálculo integral. Se aprenden técnicas que obtienen soluciones analíticas o soluciones exactas de integrales definidas e indefinidas. En esta parte se trata de solucionar integrales definidas, o sea integrar una función entre un par de límites dados [a, b] . Integral en la cual el intervalo de integración [a, b] , es finito, y f es una función de una variable real y valor real continua en [a, b] . Una integral definida se define geométricamente como el área bajo la curva f (x) en el intervalo [a, b] . De acuerdo al teorema fundamental del calculo b

integral la ecuación se evalúa como

∫ f ( x)dx = F ( x) a b

. En donde F(x) es la

a

dF ( x) = F ' ( x) = f ( x) . Es dx decir F(x) es una antiderivada de f (x) . La nomenclatura de F ( x) ba es integral de f (x) , esto es, cualquier función tal que

F (b) − F (a) .

INTRODUCCIÒN Desafortunadamente en la mayoría de los casos prácticos es muy difícil o aun imposible hallar una antiderivada de f(x). En estos casos el valor de la integral debe de aproximarse. Esto puede lograrse de las siguientes maneras: Serie de potencias. Método gráfico. Métodos numéricos. Para realizar el cálculo de una integral definida por modelos ó métodos numéricos, además de aplicar la regla Trapezoidal o Rectangular con segmentos cada vez más pequeños, otra manera de obtener una estimación más exacta de una integral, es la de usar polinomios de orden superior para conectar los puntos. Por ejemplo si hay un punto medio extra entre f (a ) y f (b) , entonces se puede conectar los tres puntos con una parábola. A las formulas resultantes para calcular la integral bajo estos polinomios se llama Reglas de SIMPSON. 1

CHAPRA Steven. Canale Raymond. Metodos numericos para ingenieros. Pag 415.McGraw - Hill

Universidad Surcolombiana – Neiva – Huila - Colombia

2 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

Dicho de otra manera se dice que para cada aplicación de la regla de SIMPSON se requieren dos subintervalos, a fin de aplicarla n número de veces, deberá dividirse el intervalo (a,b) en un número de subintervalos o segmentos. Cada subintervalo sucesivo se aproxima por un polinomio de segundo grado (parábola) y se integra de tal manera que la suma de las áreas de cada segmento de la parábola sea la aproximación a la integración deseada.

OBJETIVOS GENERAL Resolver el problema de cálculo del área bajo la curva entre dos límites conocidos, dividiendo en N sub áreas para calcular su valor, asumiendo cada sub área como un pequeño arco de parábola. 1. Comprender las bases conceptuales de la integración aproximada. 2. Comprender los rasgos generales de la integración aproximada utilizando el método de Simpson. 3. Comprender la aproximación del error por truncamiento de la integración aproximada utilizando el método de Simpson, frente al valor exacto. 4. Resolver problemas de integración numérica y apreciar su aplicación en la solución de problemas de ingeniería, utilizando el método de Simpson.

ESPECÍFICOS 1. Conocer la interpretación geométrica de la integral definida. 2. Reconocer que el método de Simpson representa, geométricamente, el área bajo una función polinomial de segundo orden (Cuadrática o Parabólica). 3. Deducir la fórmula de Simpson a partir de la interpretación geométrica de la integral definida. 4. Acotar el error cometido en la integración numérica por el método de Simpson. 5. Explicar la obtención de fórmulas más precisas para calcular, numéricamente, integrales definidas. 6. Aplicar el método de Simpson, para calcular numéricamente, las aproximaciones de algunas integrales definidas.

OBSERVACIONES PRELIMINARES Cuando se realiza un experimento, generalmente, se obtiene una tabla de valores que se espera, tengan un comportamiento funcional.

Universidad Surcolombiana – Neiva – Huila - Colombia

3 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

Sin embargo, no se obtiene la representación explícita de la función que representa la regla de correspondencia entre las variables involucradas. En estos casos, la realización de cualquier operación matemática sobre la nube de puntos que pretenda tratarla como una relación funcional, tropezará con dificultades considerables al no conocerse la expresión explícita de dicha relación. Entre estas operaciones se encuentra la integración de funciones. Además, es conocido que existen relativamente pocas fórmulas y técnicas de integración, frente a la cantidad existente de funciones que se pueden integrar. Es decir, un gran número de integrales de funciones elementales no puede ser expresada en términos de ellas. Entre estos casos singulares se tienen, a manera de ejemplo: 2

x ∫ e dx,

2

dx

∫ ln( x) , ∫

ex ∫ x dx, 1

2

1 + x 3 dx,

ln( x) ∫ x + 1 dx, 1

π /4

∫ sin( x

2

)dx,

1 + x 4 dx,

π /2

∫ x tan( x)dx, ∫ 0

∫

1

sin( x) dx,

0

1

1

∫0 1 + x 5 .dx

∫ (9 − x

1/ 3 2

)

dx

0

Lo anterior motiva el uso de los métodos de integración numérica que se estudian en lo que sigue; lo primero en considerar se basa en la aproximación de la función f mediante polinomios interpolantes. Para aclarar la contradicción antes señalada, se debe recordar la condición necesaria para que una función sea integrable. Dicha condición se menciona de inmediato, sin demostración:

Proposición 1 (Condición necesaria de Integrabilidad). Si una función f es continua en el intervalo [a, b], entonces la función f es integrable en el intervalo [a, b]. No obstante que las condiciones de la proposición 1 son sumamente generales, no se tiene garantía de que, al aplicar los métodos usualmente conocidos para resolver integrales, se pueda encontrar la antiderivada de una función f(x) cualquiera necesaria para obtener la integral definida. Estos apuntes pretenden ilustrar al lector de forma detallada y lo mas sencillo posible, una de las técnicas básicas que permiten resolver dicha situación, haciendo uso de los métodos o modelos numéricos, a través de la denominada “INTEGRACIÓN APROXIMADA, POR EL MÉTODO DE SIMPSON”.

Universidad Surcolombiana – Neiva – Huila - Colombia

4 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

CÁLCULO DE ÁREAS Uno de los problemas matemáticos más frecuentes es el cálculo del área que se forma entre una función f(x), el eje x y los límites a y b. Por ejemplo, se necesita calcular el área A que aparece en la Fig. 1, reiterando que dicha área esta por debajo de la función f(x) entre los límites a y b:

Fig. 1

Partiendo del hecho que la función f (x) y los valores a y b son conocidos. a se considera como el limite inferior y b se considera como límite superior. En este tipo de problemas se pueden obtener dos tipos de soluciones: 9 Soluciones algebraicas: se obtiene una fórmula precisa y exacta para el área solicitada. 9 Soluciones numéricas: se calcula numéricamente una estimación del área. Desde luego, la soluciones algebraicas son mejores que las numéricas, porque son exactas. Pero a veces, la complejidad de las funciones hace imposible (o difícil) obtener la solución algebraica, por lo que una solución numérica permite ahorrar tiempo.

EL MÉTODO DE SIMPSON Además de aplicar la regla trapezoidal o Rectangular con segmentos o sub áreas cada vez más pequeñas, otra manera de obtener una estimación aún más exacta de una integral, es la de usar polinomios de orden superior para conectar los puntos, en el caso particular del método que usa orden 2, es decir de la forma ax 2 + bx + c . A las fórmulas resultantes de calcular la integral bajo estos polinomios se les conoce como reglas de Simpson. En este procedimiento, se toma el intervalo de anchura 2h, comprendido entre xi y xi + 2 , y se sustituye la función f(x) por la parábola que pasa por tres puntos (xi, yi), (xi+1, yi+1), y (xi+2, yi+2). El valor del área aproximada, sombreada en la

Universidad Surcolombiana – Neiva – Huila - Colombia

5 de 22

Ing Yamil Armando Cerquera Rojas

figura,

se

calcula

con

un

[email protected]

poco

más

de

trabajo

y

el

resultado

es

h [ yi + 4 yi+1 + yi +2 ] , que se demuestra en seguida. 3 DESARROLLO DEL MODELO DE SIMPSON: Deducción del modelo a partir de la Ecuación de la parábola Para efectos de la demostración del método de Simpson, se asume cada sub área 2 como un pequeño arco de parábola de la forma ax + bx + c con límites así: Limite inferior en –h, limite superior en h, por ende la mitad de la pequeña sub área se encontrará en el Punto 0, tal como se ilustra en Fig. 2.

Fig. 2

Se procede a integrar dicho arco de parábola entre los límites descritos se tendrá: h

ax 3 bx 2 ∫−h(ax + bx + c)dx 3 + 2 + cx −h = , h

2

Reemplazando cada uno de los límites, se tiene:

 ah 3 bh 2   ah 3 bh 2   3 + 2 + ch − − 3 + 2 − ch , ahora destruyendo paréntesis se     tendrá:

ah 3 bh 2 ah 3 bh 2 ah 3 + + ch + − + ch = 2 + 2ch , simplificando un poco la 3 2 3 2 3

solución se obtendrá la ecuación 1 que se muestra a continuación.

Universidad Surcolombiana – Neiva – Huila - Colombia

6 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

h

∫ (ax

2

+ bx + c)dx =

−h

h [2ah 2 + 6c] 3

Ec 1

Fig. 3

Observando la Fig 3, en lo que respecta a las notaciones, se puede decir que f ( xi ) = yi = f (−h) , f ( xi +1 ) = yi +1 = f (0) , f ( xi + 2 ) = yi + 2 = f (h) , Entonces se podría obtener el siguiente sistemas de ecuaciones, evaluando la ecuación 2 general de la parábola ax + bx + c en cada uno de los puntos de la pequeña sub área [–h,0-h]:

f (−h) = ah 2 − bh + c , se puede tomar esta altura como y0 = f ( xi ) f (0) = c , se toma esta altura como y1 = f ( xi +1 ) f (h) = ah 2 + bh + c , y esta altura como y 2 = f ( xi + 2 ) De lo anterior se puede decir que:

y 0 + y 2 = 2ah 2 + 2c

Ec 2

y1 = c

Ec 3

Retomando la Ec 1 se puede expresar igualmente de la siguiente manera: h

∫ (ax

2

+ bx + c)dx =

−h

h [2ah 2 + 2c + 4c] 3

Ec 4

Reemplazando las ecuaciones 2 y 3 en la Ec 4 se tiene que: h

∫ (ax

2

+ bx + c)dx =

−h

Universidad Surcolombiana – Neiva – Huila - Colombia

h [ y0 + 4 y1 + y2 ] = A1 3

Ec 5

7 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

Interpretando la ecuación Ec 5 con base en la sub área seleccionada A1 para desarrollar el modelo de Simpson, se diría que el área del segmento es igual a la suma de la altura o función evaluada en el lado izquierdo mas cuatro veces la función evaluada en la parte central de la sub área mas la función evaluada en el lado derecho de la sub área, todo esto multiplicado por el ancho del sub área y dividido por 3. La simple inspección visual de esta figura y la que describe el procedimiento de los trapecios o los rectángulos, confirma que el método de Simpson deberá ser mucho más exacto que los procedimientos mencionados. Si a y b se denominan como x0 y x2 , y f i ( xi ) se representa mediante un polinomio de Lagrange de segundo orden, entonces la integral es:

 ( x − x1 )( x − x2 )  ( ) f x + 0  ( x − x )( x − x )  1 0 2  0  X2  ( x − x0 )( x − x2 )  I= ∫ f ( x1 ) + dx X 0 ( x1 − x0 )( x1 − x 2 )    ( x − x0 )( x − x1 )  f ( x2 )    ( x2 − x0 )( x2 − x1 )  Después de integrar y de reordenar los términos, resulta la siguiente ecuación:

f ( x0 ) + 4 f ( x1 ) + f ( x2 ) 6

I = (b − a )

Ec 5a

Deducción del modelo a partir del polinomio de lagrange Suponga que se tienen los datos:

a

xm

b

f (a)

f ( xm )

f (b)

Donde x m es el punto medio entre a y. b En este caso se tiene que: b

∫ a

b

f ( x)dx ≈ ∫ f 2 ( x)dx a

Donde f 2 ( x ) es el polinomio de interpolación para los datos en la tabla anterior. Se usará el polinomio de Lagrange.

Universidad Surcolombiana – Neiva – Huila - Colombia

8 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

Así, se tiene que:

f 2 ( x) = f (a ) Si se denota h =

( x − a )( x − x m ) ( x − x m )( x − b) ( x − a )( x − b) + f ( xm ) + f (b) (b − a )(b − x m ) (a − x m )(a − b) ( x m − a )( x m − b)

b−a = x m − a = b − x m , entonces: 2

f 2 ( x) = f (a )

( x − x m )( x − b) ( x − a)( x − x m ) ( x − a)( x − b) + f ( xm ) + f (b) (− h)(−2h) (h)(− h) (2h)(h)

Simplificando términos:

f 2 ( x) =

f ( xm ) f (a) f (b) ( x − x m )( x − b) − ( x − a)( x − b) + ( x − a )( x − x m ) 2 2 2h h 2h 2

Se puede notar que cada uno de los términos anteriores, son esencialmente de la misma forma, es decir, una constante por ( x − α )( x − β ) . Así, se calcula la siguiente integral por partes:

∫ ( x − α )( x − β )dx Sea:

u = x −α

du = dx

dv = ( x − β )dx

(x − β )2 v = ∫ ( x − β )dx = 2

Por lo tanto,

∫ ( x − α )( x − β )dx = ( x − α )

(x − β )2 (x − β )2 −∫ dx 2 2

∫ ( x − α )( x − β )dx = ( x − α )

(x − β ) 2 (x − β )3 − 2 6

Se usa esta fórmula para calcular la integral de cada uno de los tres términos de

f 2 ( x) .

f ( xm ) b f (a) b f (b) b ( x − x ) ( x − b ) dx − ( x − a ) ( x − b ) dx + m ∫ ∫ ∫ ∫ ( x − a)( x − xm )dx 2h 2 a h2 a 2h 2 a a Tomando cada término a integrar como I x , se tendría la siguiente expresión. b

f 2 ( x)dx =

b

∫ f 2 ( x)dx = a

f ( xm ) f (a) f (b) I − I2 + I3 2 1 2 2h h 2h 2

Ec. 1

Entonces b

I 1 = ∫ ( x − x m )( x − b)dx = (x − x m ) a

( x − b) 2 ( x − b ) 3 b ]a − 2 6

Universidad Surcolombiana – Neiva – Huila - Colombia

9 de 22

Ing Yamil Armando Cerquera Rojas

I 1 = −(a − x m )

[email protected]

( a − b) 2 ( a − b) 3 ( −2 h ) 2 ( −2 h ) 3 + = −( − h) + 2 6 2 6

4 I 1 = 2h 3 − h 3 3

I1 = b

I 2 = ∫ ( x − a)( x − b)dx = ( x − a ) a

I2 =

2 3 h 3

( x − b ) 2 ( x − b) 3 b − ]a 2 6

(a − b) 3 (−2h) 3 4 = = − h3 b 6 3 4 I 2 = − h3 3 b

I 3 = ∫ ( x − a)( x − x m )dx = ( x − a ) a

I 3 = (b − a )

( x − xm ) 2 ( x − xm ) 3 b − ]a 2 6

(b − x m ) 2 (b − x m ) 3 (a − x m ) 3 − + 2 6 6

h 2 h 3 ( − h) 3 h 3 2h 3 3 I 3 = (2h ) − + =h − = 2 6 6 3 3 I3 =

2h 3 3

Ahora retomando la Ec. 1 se tiene: b

∫ f 2 ( x)dx = a

b

∫ f 2 ( x)dx = a

f ( xm ) f (a) f (b) I − I2 + I3 2 1 2 2h h 2h 2 f (a )  2 3  f ( x m )  4 3  f (b)  2 3   h  − h  +  h − 2h 2  3  h 2  3  2h 2  3 

b

∫ f 2 ( x)dx =

f (a )

a

4 h h + f ( x m ) h + f (b) 3 3 3

h

∫ f 2 ( x)dx ≅ 3 [ f (a) + 4 f ( xm ) + f (b)] Debido al factor

1 h se le conoce como la regla de Simpson de un tercio. 3

Universidad Surcolombiana – Neiva – Huila - Colombia

10 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

b−a para obtener la fórmula final: 2 f (a) + 4 f ( x m ) + f (b) f ( x)dx ≅ (b − a ) 6

En la práctica, se sustituye el valor de h = b

∫ a

Generalización del modelo para n subintervalos Si se toma (b − a ) / 6 ≈ h / 3 , f ( x0 ) = y 0 , f ( x1 ) = y1 , y f ( x2 ) = y 2 , entonces se tiene como solución de la sub área I = mostrado en la ecuación 5.

h ( y0 + 4 y1 + y2 ) , que sería lo mismo 3

Ahora, se sabe que el área que se desea encontrar sería la sumatoria de todas las sub áreas que se calculen. Al igual que los métodos de la regla trapezoidal y de la regla rectangular, entre mas sub áreas tenga la integral a calcular, mas exacto será el valor encontrado. El área aproximada en el intervalo [a, b] es: b

∫ f ( x)dx = A

1

+ A2 + A3 + ... + An

a

Ahora dejando esta ecuación en términos de la ecuación 5 se tendrá:

h

b

∫ f ( x)dx = 3 ( y

0

a

h h + 4 y1 + y 2 ) + ( y 2 + 4 y3 + y 4 ) + ... + ( y 2 n− 2 + 4 y 2 n −1 + y 2 n ) 3 3

Simplificando h / 3 y sumando los términos se tendrá:

h

b

∫ f ( x)dx = 3 ( y

0

+ 4 y1 + 2 y 2 + 4 y3 + 2 y 4 + 4 y5 + 2 y6 + ... + 4 y 2 n −1 + y 2 n ) ,

a

Donde n sería el número de sub áreas en el cual se ha dividido el área que se desea calcular. A manera de ejemplo, si el área a calcular se hubiera dividido en 4 Sub áreas entonces en términos de y la solución seria:

h

b

∫ f ( x)dx = 3 ( y

0

+ 4 y1 + 2 y 2 + 4 y3 + 2 y 4 + 4 y5 + 2 y6 + 4 y7 + y8 )

a

Bien, dependiendo como se agrupen los términos se llegaría a expresar la solución de dos maneras: b

h

∫ f ( x)dx = 3 [ y

0

+ y8 + 4( y1 + y3 + y5 + y7 ) + 2( y 2 + y 4 + y6 )]

Ec 6

a

ó Universidad Surcolombiana – Neiva – Huila - Colombia

11 de 22

Ing Yamil Armando Cerquera Rojas

b

h

∫ f ( x)dx = 3 [ y

0

[email protected]

− y8 + 4( y1 + y3 + y5 + y7 ) + 2( y 2 + y 4 + y 6 + y8 )]

Ec 7

a

Los primeros términos del paréntesis, contienen los valores de la evaluación de la función en los extremos, el segundo, la suma de los términos de índice impar, y el tercero la suma de los términos de índice par. Las dos ecuaciones se pudieran representar en términos de sumatorias de la siguiente manera. La Ec 6: b

∫ a

n n −1 h f ( x)dx =  y 0 + y 2 n + 4∑ y 2 i −1 + 2∑ y 2 i   3  i =1 i =1

Ec 8

n h f ( x)dx =  y 0 − y 2 n + ∑ [4 y 2 i −1 + 2 y 2 i ]  3  i =1

Ec 9

La Ec 7: b

∫ a

Para efectos de programación y en lo que respecta a mi concepto personal, es mejor la solución representada como Ec 9 y con ella se continúa el trabajo. Hay que tener en cuenta que n es el número de sub áreas en la que se divide el área total a calcular y h = dx / 2 . Ahora lo que se conoce en un momento determinado, cuando se desea calcular el valor de la integral definida, son los siguientes términos:

a = Límite inferior b = Límite Superior n = Número de sub áreas f ( x) La función sobre la cual se desea integrar. Con los valores anteriores se pude calcular el valor de dx así: dx = (b − a ) / n y h = dx / 2 . En necesario entonces dejar la ecuación en términos de f ( x ) , a , b y dx ó h así: Los primeros términos: y 0 = f (a ) y y 2 n = f (b) Analizando ahora los términos impares: y1 = f ( a + 1dx / 2) , y3 = f ( a + 3dx / 2) ,

y5 = f (a + 5dx / 2) , por tanto se tendría de manera general:

y 2i −1 = f (a + (2i − 1)dx / 2) ó y 2i −1 = f (a + (2i − 1)h)

Universidad Surcolombiana – Neiva – Huila - Colombia

Ec 10

12 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

Analizando ahora los términos pares: y 2 = f ( a + 1dx) , y 4 = f (a + 2dx) ,

y 6 = f (a + 3dx) , por tanto se tendría de manera general:

y 2i = f (a + idx) ó y 2i = f (a + 2ih)

Ec 11

Dejando la ecuación Ec 9 en términos de lo expresado en las ecuaciones Ec 10 y Ec 11 se tendrá en forma definitiva la solución así: b

∫ a

n h f ( x)dx =  f (a ) − f (b) + ∑ [4 f (a + (2i − 1)dx / 2 + 2 f (a + idx)]  3  i =1

Ec 12

Ejemplos Ejemplo 1: Utilizar la regla de Simpson para aproximar la integral:

1

2

x ∫ e dx . 0

Tenga en cuenta que el valor real es 1.4626…

Fig. 4

Solución: Usando la fórmula directamente con los siguientes datos:

a=0 b =1

f ( x) = e x

2

Si se asume el área a calcular como un solo arco de parábola, se tendría entonces que dx = (b − a ) / 1 = 1 ó h = ( dx / 2) = 0.5 y aplicando la ecuación Ec 5 se tiene que: 1

x ∫ e dx =

0

Universidad Surcolombiana – Neiva – Huila - Colombia

2

h [ y0 + 4 y1 + y2 ] = A 3

13 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

1

∫e

x2

dx =

0

1

∫e

x2

dx =

0

[

2 2 0 .5 0 2 e + 4 e 0. 5 + e 1 3

]

0.5 [1 + 4(1.2840) + 2.7183] 3 1

∫e

x2

dx =1.4757

0

Ahora si compara los resultados obtenidos al aplicar la regla del Trapecio o la regla de los Rectángulos, con respecto al valor real y al valor obtenido por la regla de Simpson tendría que analizar lo siguiente: Integral

Valor Real

Rectangular

Trapezoidal

Simpson

f ( x) = e

1.4626

1.2840

1.8591

1.4757

12.21% 0.1786

27.11% 0.3965

0.90% 0.0131

x2

Er Ea

Vale la pena aclarar que para los tres métodos se trabajó una sola sub área. Desarrollado en MatLab se tendría el siguiente resultado. »syms x »f=exp(x^2); »integral=int(f) integral = -1/2*i*pi^(1/2)*erf(i*x) ERF Error de la función. Y = ERF(X) es el error de la función para cada elemento de X. X debe ser real. El error de la función está definido como: erf(x) = 2/sqrt(pi) * integral desde 0 a x de exp(-t^2) dt. Analice lo anterior. 1

Ejemplo 2: Aplicar la regla de Simpson para aproximar la integral ∫ e x dx si 2

0

se subdivide el área total en 5 intervalos. Solución: En este caso, se identifica n = 5 , y las particiones generadas estarían delimitadas por los puntos P={0.0, 0.2, 0.4, 0.6, 0.8, 1.0} sobre el eje x. Así, aplicando la fórmula:

Universidad Surcolombiana – Neiva – Huila - Colombia

14 de 22

Ing Yamil Armando Cerquera Rojas

b

h

∫ f ( x)dx = 3 ( y

0

[email protected]

+ 4 y1 + 2 y 2 + 4 y3 + 2 y 4 + 4 y5 + 2 y6 + 4 y 7 + 2 y8 + 4 y9 + y10 )

a

Si se asume el área a calcular como cinco pequeños arcos de parábola, se tendría entonces que dx = (b − a ) / 5 = 0.2 ó h = ( dx / 2) = 0.1 y aplicando la ecuación Ec 7 se tiene que: b

h

∫ f ( x)dx = 3 [ y

0

− y10 + 4( y1 + y3 + y5 + y7 + y9 ) + 2( y 2 + y 4 + y6 + y8 + y10 )]

a

y 0 = f (a) = e a = e 0 = 1 , y 2 n = y10 = f (b) = e1 = 2.7183 2

2

2

Los términos impares se encontrarían de acuerdo con la fórmula:

y 2i −1 = f (a + (2i − 1)h) = f ( xi ) Así: y1 = f (h) = f (0.1) = e 0.1 = 1.0101 2

y3 = f (3h) = f (0.3) = e 0.3 = 1.0942 2

y5 = f (5h) = f (0.5) = e 0.5 = 1.2840 2

y7 = f (7h) = f (0.7) = e 0.7 = 1.6323 2

y9 = f (9h) = f (0.9) = e 0.9 = 2.2479 2

Y la sumatoria igual a: 7.2685 Los términos pares se encontrarían de acuerdo a la fórmula y 2 i = f (a + 2ih) así:

y 2 = f (a + 2h) = f (0.2) = e 0.2 = 1.0408 2

y 4 = f (a + 4h) = f (0.4) = e 0.4 = 1.1735 2

y 6 = f (a + 6h) = f (0.6) = e 0.6 = 1.4333 2

y8 = f (a + 8h) = f (0.8) = e 0.8 = 1.8965 2

y10 = f (a + 10h) = f (1) = e1 = 2.7183 2

Y la sumatoria igual a: 8.2624 Por tanto el valor de la integral será igual a: b

∫ f ( x)dx = a

0.1 [1 − 2.7183 + 4(7.2685) + 2(8.2624)] = 1.4627 3

Universidad Surcolombiana – Neiva – Huila - Colombia

15 de 22

Ing Yamil Armando Cerquera Rojas

Integral

f ( x) = e x

2

[email protected]

Valor Real

Rectangular

Trapezoidal

Simpson

1.4626

1.4537

1.48065

1.4627

0.6085% 0.0089

1.2375% 0.0181

0.0068% 0.0001

Er Ea

ex Ejemplo 3: Usar la regla de Simpson para aproximar la integral: ∫ dx . 2 x 4

Solución: Igual que en el ejemplo anterior, se sustituyen los datos de manera directa en la fórmula de Simpson dividiendo el área en cuatro (4) sub áreas. En este caso, se tienen los datos: a=2, b=4, n=4

f ( x) = e x / x dx = (4-2)/4 = 0.5 h = dx/2 = 0.25 Por lo tanto, se tiene que:

h

b

∫ f ( x)dx = 3 [ y

− y10 + 4( y1 + y3 + y5 + y7 + y9 ) + 2( y 2 + y 4 + y6 + y8 + y10 )]

0

a

y 0 = f (a ) = e a / a = e 2 / 2 = 3.6945 , y 2 n = y8 = f (b) = e 4 / 4 = 13.6495 Los

términos

y 2i −1 y1 = y3 = y5 = y7 =

impares

se

encontrarían = f (a + (2i − 1)h) = f ( xi ) así:

de

acuerdo

con

la

fórmula

f (a + h) = f (2 + 0.25) = e 2.25 / 2.25 = 4.2168 f (a + 3h) = f (2.75) = e 2.75 / 2.75 = 5.6882 f (a + 5h) = f (3.25) = e 3.25 / 3.25 = 7.9355 f (a + 7 h) = f (3.75) = e 3.75 / 3.75 = 11.3390

Y la sumatoria igual a: 29.1795 Los términos pares se encontrarían de acuerdo a la fórmula y 2 i = f (a + 2ih) así:

y2 y4 y6 y8

= = = =

f ( a + 2h) = f ( a + 4h) = f ( a + 6h) = f ( a + 8h ) =

f (2.5) = e 2.5 / 2.5 = 4.8730 f (3.0) = e 3.0 / 3.0 = 6.6952 f (3.5) = e 3.5 / 3.5 = 9.4616 f (4.0) = e 4.0 / 4.0 = 13.6495

Y la sumatoria igual a: 34.6793 Por tanto el valor de la integral será igual a:

Universidad Surcolombiana – Neiva – Huila - Colombia

16 de 22

Ing Yamil Armando Cerquera Rojas

b

∫ f ( x)dx = a

[email protected]

0.25 [3.6945 − 13.6495 + 4(29.1795) + 2(34.6793)] = 14.6768 3

Ejemplo 4: Evaluar la función

1.2

∫ f ( x)dx , usando la siguiente tabla: 0

x f(x)

0 0

0.10 6.84

0.30 4.00

0.50 4.20

0.70 5.51

0.95 5.77

1.20 1.0

Fig 5

Solución. Observe en la fig 5 que en el intervalo [0,0.1] se puede aplicar la regla del trapecio, en el intervalo [0.1,0.7] la regla de Simpson de 3/8 y en el intervalo [0.7,1.2] la regla de Simpson de 1/3. Así, se tienen las siguientes integrales:

I1 = I2 = I3 =

0.1

∫ f ( x)dx =

0 0.7

0.1 − 0 [ f (0) + f (0.1)] = 0.842 2

∫

f ( x)dx =

0.7 − 0.1 [ f (0.1) + 3 f (0.3) + 3 f (0.5) + f (0.7)] = 2.7712 8

∫

f ( x)dx =

1.2 − 0.7 [ f (0.7) + 4 f (0.95) + f (1.2)] = 2.4658 6

0.1 1. 2 0. 7

Finalmente, la integral buscada es la suma de las tres integrales anteriores: 1.2

∫ f ( x)dx = 0.842 + 2.7712 + 2.4658 = 6.079 0

Universidad Surcolombiana – Neiva – Huila - Colombia

17 de 22

Ing Yamil Armando Cerquera Rojas

Ejemplo 5: Calcula la integral

[email protected]

3.25

∫ f ( x)dx , usando la siguiente tabla de datos:

−1

x f(x)

-1 2

-0.5 -3

0 1.5

1 -1

1.75 0.5

2.5 0.75

3.25 -2

Si se desarrolla la gráfica con Matlab ajustando los datos de la tabla anterior mediante splines cúbicos, se tendría lo siguiente: » x=[-1 -0.5 0 1 1.75 2.5 3.25]; » y=[2 -3 1.5 -1 0.5 0.75 -2]; » xx=-1:0.01:3.5; » yy=spline(x,y,xx); » plot(x,y,'o',xx,yy)

Fig. 6

Solución combinando varios métodos de integración. Para este caso, se puede aplicar la regla de Simpson de 1/3 en el intervalo [1,0], la regla del trapecio en el intervalo [0,1] y la regla de Simpson de 3/8 en el intervalo [1,3.25]. Así, se tiene las siguientes integrales:

I1 =

0

∫ f ( x)dx ≈

−1 1

I 2 = ∫ f ( x)dx ≈

I3 =

0 3.25

0 − (−1) [ f (−1) + 4 f (−0.5) + f (0)] ≈ −1.41667 6 1− 0 [ f (0) + f (1)] ≈ 0.25 2

∫ f ( x)dx ≈ 1

3.25 − 1 [ f (1) + 3 f (1.75) + 3 f (2.5) + f (3.25)] ≈ 0.210938 8

Por lo tanto, la integral buscada es la suma de las tres integrales anteriores:

Universidad Surcolombiana – Neiva – Huila - Colombia

18 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

3.25

∫ f ( x)dx ≈ −1.41667 + 0.25 + 0.210938 ≈ −0.955729

−1

Vale la pena comentar que no siempre tiene que suceder que se apliquen exactamente las tres reglas. En realidad, esto depende de cómo se encuentran espaciados los intervalos de la tabla de datos y la forma que pueda tener la curva.

Ejemplo 6: Calcular la siguiente Integral:

3

log( x) dx , correspondiente al área 3 1

∫

bajo la curva mostrada en la Fig. 7 entre los límites 1 y 3. Intervalo: 1 ≤ x ≤ 3 Método: Regla de Simpson n = 2 sub intervalos, donde dx = (b-a)/n=1.

Se utiliza la regla de Simpson con dx =1.0 ó h=0.5, n = 2 sub intervalos y la tabla 3

de valores para f ( x) = ∫ 1

ln( x) x3

dx

Fig. 7

x 1.0 1.5 2.0 2.5 3.0

TABLA de valores

f (x)

yi

0.00000000000000 0.22070724963720 0.24506453586714 0.23180525693000 0.21142803353252

Y0 Y1 Y2 Y3 Y4

Y5

Universidad Surcolombiana – Neiva – Huila - Colombia

19 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

El resultado aplicando la siguiente fórmula sería: (Tenga en cuenta que solo se trabajan 2 sub áreas) b

I = ∫ f ( x)dx = a

h [ y0 − y4 + 4( y1 + y3 ) + 2( y2 + y4 )] 3

I= 0.50/3 * (0.000000000000000 - 0.21142803353252 +… 4 * (0.22070724963720 + 0.23180525693000) +… 2 * (0.24506453586714 + 0.21142803353252))

I = 0.41860118858927 Programa en diferentes lenguajes Programación del método de Simpson en lenguaje C. Suponga que la función f evalúa la función f(x). Entonces la siguiente fórmula permite calcular el área de cada una de las secciones: b

∫ a

n h f ( x)dx =  f (a ) − f (b) + ∑ [4 f (a + (2i − 1)dx / 2 + 2 f (a + idx)]  3  i =1

Programa en lenguaje C: Se supone que la función area se debe llamar con los parámetros a,b y n, que serían límite inferior, límite superior y número de sub áreas respectivamente. double area(double a, double b, int n) { double dx= (b-a)/n; double sum= f(a)-f(b); for (int i=1; i<=n; i++) { double ximpar=a+(2*i-1)*dx; double xpar= (a+i*dx); sum+= 4*f(ximpar)+2*f(xpar); } return dx*sum/6; } double f(double x) { double y; // Se describiría la función a la cual se le desea calcular la integral y = ... return y; } Esta solución es análoga a la que se programó para buscar los ceros de una función. El problema de esta solución, es que cuando se desea calcular la integral de varias funciones distintas, hay que programar una función para calcular el área de cada función.

Universidad Surcolombiana – Neiva – Huila - Colombia

20 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

Si se programa para MatLab el programa sería function area=simpson(a,b,n) dx=(b-a)/n; sum= f(a)-f(b) for i=1:n, ximpar=a+(2*i-1)*dx/2 xpar= a+i*dx sum=sum+ 4*f(ximpar)+2*f(xpar); end area= dx*sum/6 Y la función f estaría dada por (Como caso particular se coloca x^2-4, pero solo cambie la función f y le integra lo que desee):

function y=f(x) y =x^2-4; %Puede cambiar esta función

La jerarquía de clases para C. Se crea una clase base abstracta denominada Simpson, que defina la función miembro integral que calcula la integral definida de cualquier función f(x) por el procedimiento de Simpson. public abstract class Simpson { public double integral(double a, double b, int n) { double dx=(b-a)/n; double suma=f(a)-f(b); for(int i=1; i<=n; i++) { suma+=4*f(a+(2*i-1)*dx) + 2*f(a+i*dx); } return (suma*dx/6); } abstract public double f(double x); } En la clase derivada Funcion de Simpson se definirá la función f(x) cuya integral se desea calcular. public class Funcion extends Simpson { public double f(double x) { return Math.cos(x); } } Para hallar la integral definida de esta función entre los límites 0 y p /2, se crea un objeto de la clase Funcion mediante new y se llama desde este objeto a la función integral pasándole en el primer argumento el límite inferior “0”, en el

Universidad Surcolombiana – Neiva – Huila - Colombia

21 de 22

Ing Yamil Armando Cerquera Rojas

[email protected]

segundo el límite superior, “pi/2”, y por último, el número de divisiones del intervalo “10”. double resultado=new Funcion().integral(0.0, Math.PI/2, 10); System.out.println("integral "+resultado); Comparando los resultados obtenidos por este procedimiento y por el procedimiento de los trapecios se puede comprobar la mayor exactitud de éste último.

Universidad Surcolombiana – Neiva – Huila - Colombia

22 de 22