Método del Trapecio en Scilab, Código Integración Numérica.

[EN ESTE BLOG LOS PROGRAMAS SE ENCUENTRAN AL FINAL PARA SU DESCARGA]

El método del trapecio y también conocido como ‘Regla de Simpson‘, aunque esta última en realidad es una ‘optimización’ de la regla de los trapecios. Pues de este método ya hemos hablado aquí por lo que esta entrada no contendrá una explicación muy amplia sobre el funcionamiento del método, por lo que te recomendamos que le des un vistazo a esa entrada y así puede que comprendas mucho mejor el algoritmo del mismo. Pues bueno, como el titulo lo indica, esta vez lo implementaremos en el conocido software de código abierto Scilab.

La esencia del método del trapecio es calcular el área que existe bajo la curva definida por y=f(x) y el eje X, en un intervalo dado [a,b] mediante el uso de Ntrapecios’ o rectángulos, además suponiendo que dicha función es (y debe ser) continua en todo este intervalo, lo anterior lo ilustramos en la figura 1.

Función Y(x), el área bajo la curva, se puede aproximar mediante n trapecios, en este caso 4.

Fig 1. Función Y(x), el área bajo la curva, se puede aproximar mediante n trapecios, en este caso 4.

En este caso nos enfocaremos en explicar el algoritmo que se sigue, así que para una explicación más detallada te invitamos a dar un vistazo aquí.

Algoritmo para el método de los trapecios.

Paso 1.

En primer lugar se divide el intervalo comprendido entre [a,b] en N sub-intervalos, nombrando al ancho de esos sub-intervalos dx (o delta de X), las variables a, b y n pueden ser ingresadas en medio de la ejecución del programa con funciones como ‘input’, etc. Recuerden que b>a.

punto fijoPaso 2.

Se realiza la siguiente serie (la sumatoria) a continuación, parece difícil pero en seguida la desarrollaremos en Scilab, recuerden que el valor de n es el número de sub-intervalos, así f(X1), f(X2),…f(Xn) es la función evaluada en cada sub-intervalo hasta el final.

punto fijo

Como se puede observar, al principio de la serie tenemos el término que denominamos dx pero divido entre 2, es decir que al inicio de la serie cuando vayamos a programar el método solo tendremos que hacer dx/2 multiplicando el resto de la serie, también podemos observar como todos los términos está multiplicados por 2 excepto el primero y el último término, esto es importante tenerlo en cuenta para entender el desarrollo del código a continuación.

Código en Scilab.

Primero que nada, recuerden que el código se encuentra al final, así que no es necesario que copien lo que aparece en las imágenes, ya que lo que buscamos es explicar el código en cada paso.

El desarrollo del programa es muy simple, inicialmente definimos la función, cuyo nombre en nuestro caso ha sido ‘trapecio()’ pero puede ser cualquiera, debe coincidir el nombre de dicha función y del archivo, esta función retornará una variable de nombre ‘area’ la cual contendrá el valor calculado que es aproximado (como todos los métodos numéricos) al valor real de la integral o el área bajo la curva de la función, como argumentos vamos a recibir los límites del intervalo [a,b] en otra palabras recibimos los valores de a y b respectivamente y como último parámetro recibimos el número de sub-intervalos en los que queremos dividir el intervalo [a,b] es decir recibimos el valor de n (cuantos rectángulos queremos usar, a mayor número de rectángulos, mayor es la precisión), lo anterior lo vemos en la siguiente figura.

trap1

A continuación calculamos el primer término de la sumatoria, aquel que como vemos, no está siendo multiplicado por el 2, y posteriormente hacemos uso de un ciclo iterativo y cuya variable de control inicia en 2 para evitar el cálculo del primer y el último termino que no entran en este ciclo, como se muestra a continuación.

trap2

Por último, terminado el ciclo, se procede a calcular el último término de la serie f(Xn) y se multiplica en la variable ‘area’ por el ‘factor común’ (dx/2) que veíamos fuera de la serie, de esta forma ya se da por terminado el método como se ve a continuación.

trap3

¡Dame un Ejemplo!

Como siempre en este sitio, antes de dejarte el código los probamos con ejemplos para mostrarte como se debe usar y que verdaderamente funciona, en este caso vamos a hacer 2 ejemplos:

Ejemplo 1.

Para este primer ejemplo vamos a calcular la integral de la función f(x)=x², en el intervalo 0 a 1, si realizamos el cálculo analítico de esta integral sería F(x)=x³/3 que evaluada en el intervalo [0,1] debe darnos como resultado 1/3 o 0.3333′, entonces a continuación definimos la función antes dicha mediante el comando deff(), como vemos ahora, importante, el nombre de la función debe ser y=f(x):


deff('y=f(x)','y=x^2;')

luego de definida la función le pasamos los argumentos a=0, b=1 y usaremos 500 sub-intervalos, es decir n=500 y damos ‘enter’, esto debe retornar el valor de la integral, como se ve en la imagen:

ejemplo1Como podemos ver, en primer lugar cargamos la función trapecio.sci al workspace de Scilab, recuerden que el archivo debe estar en el folder en el que estemos trabajando en scilab, de lo contrario nos dará un error de archivo no encontrado, en el segundo paso definimos la función y=f(x) en nuestro caso será y=x², finalmente llamamos la función con sus límites de integración y el número de sub-intervalos (n=500).

Ejemplo 2.

En este segundo ejemplo, vamos a calcular la integral entre 0 y PI de y=sin(x) que si la integramos manualmente o le preguntamos a Wolfram Alpha nos daremos cuenta que su valor real será 2:

Integral de sin(x) entre 0 y PI.

Integral de sin(x) entre 0 y PI.

Para este caso realizaremos los mismos pasos que en el anterior ejemplo, le indicaremos a la función que los límites de integración van de 0 a pi (en Scilab la contante pi la obtenemos con %pi).

ejemplo2

Como podemos observar (y era de esperarse) el resultado no es exactamente 2 y esto se debe a que justamente es un método numérico y la esencia de estos es dar una aproximación (lo más cercana posible) al resultado correcto, en este caso es bastante acertado y como ya os había mencionado si aumentamos el valor de N la precisión aumentará ya que se tomarán más puntos.

Como siempre puede descargar el código de este artículo desde aquí.

Si tienes alguna duda o comentario puedes dejarlo abajo, comparte si esto te ha servido y así logramos que esto llegue a más gente, como siempre comentar es agradecer, saludos y éxitos.

Autor: Julio Marulanda.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s