Suma de Riemann

En matemáticas, la suma de Riemann es un método de integración numérica que nos sirve para calcular el valor de una integral definida es decir el área bajo una curva, este método es muy útil cuando no es posible utilizar el Teorema Fundamental del Cálculo. Éstas sumas toman su nombre del matemático alemán Bernhard Riemann.

La suma de Riemann consiste básicamente en trazar un número finito de rectángulos dentro de un área irregular, calcular el área de cada uno de los rectángulos y sumarlos. El problema de este método de integración numérica es que al sumar las áreas se obtiene un margen de error muy grande.

Fuente: Wikipedia

Últimamente estuve trabajando en un programa para resolver cuestiones de integración numéricamente ya que fue una de las encomiendas de mi asesor el Dr. Jesús Manzanares creo que este es un “Baby Step” que nos llevará a resolver algo mucho mas complejo, creí que era relevante incluir la descripción de la integración numérica por lo que incluí la definición de wikipedia que es buena.

Cabe complementar un poco la definición y agregar que entre más rectángulos se le agreguen al problema la aproximación de la solución es mucho más exacta, y de acuerdo a pruebas hechas llegando a 100 o poco más de esa cantidad se logra tener una aproximación de la solución bastante aceptable 🙂

Bueno creo que ya fue demasiada teoría aunque realmente la teoría nunca es suficiente, a continuación se presenta una implementación en web de la Suma de Riemann, para tener información mas amplia de como realizar un programa en web puedes visitar el post ¿Como portar una aplicación existente a una solución web?

La aplicación está diseñada para resolver 3 tipos de integrales, las siguientes:

riemann.tex

(1)   \begin{equation*}     f(x) = \int_0^1{x^2}\, \mathrm{d}x \end{equation*}

(2)   \begin{equation*}     f(x) = \int_0^1{x^3}\, \mathrm{d}x \end{equation*}

(3)   \begin{equation*}     f(x) = \int_0^1 \cos(x)\, \mathrm{d}x \end{equation*}

[/latexpage]

Donde el rango de la integral es totalmente parametrizable en la implementación, y en donde para realizar la integración numérica se deberá hacer una sumatoria, a continuación se enlista el código en C++ de la implementación, seguido de la aplicación en sí corriendo en web.

NumericalIntegral.cpp

/*
 * File:   NumericalIntegral.cpp
 * Author: Yohan Jasdid
 *
 * Created on October 15, 2011, 4:36 PM
 */
#include "cstdlib"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
 
#include "algorithm"         
#include "iostream"      
#include "fstream"       
#include "string"
#include "sstream"       
#include "math.h"
 
#include "getpost.h"
#include "json/reader.h"
#include "json/writer.h"
#include "json/elements.h"
 
using namespace std;
using namespace json;
 
string exec(char* cmd);
void ReturnError(string error);
 
/**
 * Main function
 * @return
 */
int main() {
    map Post;
    initializePost(Post);   //notice that the variable are passed by reference!
 
    cout << "content-type: application/json" << endl << endl;
 
    string userData;
    string message = "";
    if (Post.find("UserData") != Post.end()) {
        userData = Post["UserData"].c_str();
    } else {
        ReturnError("Some requiered parameters were not provided, try again.");
        return 1;
    }
 
    Object objDocument;
    try {
        std::istringstream jsonData(userData);
        Reader::Read(objDocument, jsonData);
 
        // Get and parse parameters
        String jsonSX = objDocument["SX"];
        String jsonEX = objDocument["EX"];
        String jsonN = objDocument["N"];
        String jsonFnc = objDocument["Fnc"];
 
        double start_x = atof(jsonSX.Value().c_str());
        double end_x = atof(jsonEX.Value().c_str());
        double interval = end_x - (start_x);
        double max_n = atof(jsonN.Value().c_str());
        int fnction = atoi(jsonFnc.Value().c_str());
 
        // We must generate a unique ID for naming our data graph file
        char uuidgen[] = "uuidgen";
        string uuid = exec(uuidgen);
 
        // Define files and folders
        string gpFolder = "../gnuplot/";
        string gpDataFile = gpFolder + uuid + ".dat";          // data file
        string gpGraphFile =  gpFolder + uuid + ".png";        // graph file
        string gpCommandFile =  gpFolder + uuid + ".gnuplot";  // gnuplot command file
 
        // Create data file
        ofstream gpData;
        gpData.open((gpDataFile).c_str());
        if(!gpData) {
            ReturnError("** Failed to create the graph file.");
            return 1;
        }
 
        // search for all n..
        stringstream data;
        for (double i = 1; i <= max_n; i++) {
            double n = i;
            double delta_x = interval / n;
            double x_i = start_x;
            double total = 0.0;
            double function = 0.0;
 
            for (double j = 1; j <= n; j++) {
                x_i = x_i + delta_x;
                switch (fnction) {
                    case 1:
                        function = x_i * x_i;
                        break;
                    case 2:
                        function = x_i * x_i * x_i;
                        break;
                    case 3:
                        function = cos(x_i);
                        break;
                }
                total = total + ((function) * delta_x);
            }
 
            data << data << i << " " << total << endl;
            gpData << i << " " << total << "\n";
        }
        gpData.close();
 
        // We must build a gnuplot commands file
        ofstream gpCommands;
        gpCommands.open((gpCommandFile).c_str());
        if(!gpCommands) {
            ReturnError("** Failed to create the graph file.");
            return 1;
        }
 
        gpCommands << "set title 'Integral Numerica'" << "\n";
        gpCommands << "set term png small" << "\n";
        gpCommands << "set output " << "'" << gpGraphFile << "'" << "\n";
        gpCommands << "set nokey" << "\n";
        gpCommands << "plot '" << gpDataFile << "' using 1:2 with lines" << "\n";
        gpCommands.close();
 
        // Execute gnuplot and generate graph
        string commds = "gnuplot '" + gpCommandFile + "'";
        char *cmds = (char*)commds.c_str();
        exec(cmds);
 
        // Results ***
        // Build a json object to return as response of the request
        message += "Calculations finished without errors.";
        Object objResponse;
        objResponse["Error"] = String("");
        objResponse["Message"] = String(message.c_str());
        objResponse["DataFile"] = String(gpDataFile.c_str());
        objResponse["CommandFile"] = String(gpCommandFile.c_str());
        objResponse["GraphFile"] = String(gpGraphFile.c_str());
        objResponse["PlainData"] = String(data.str());
 
        std::stringstream jsonStream;
        Writer::Write(objResponse, jsonStream);
        cout << jsonStream.str();
    } catch (exception& e) {
        ReturnError("** Something went wrong, please report this problem to the system administrator.");
        return 1;
    }
    return 0;
}
 
/**
 * @brief Returns an error in json format
 * @error The error description to be returned
 *
 * This method executes a bash command using a pipe and a buffer to capture the
 * result of the execution, then return the execution to client code.
 */
void ReturnError(string error) {
  Object objResponse;
  objResponse["Error"] = String(error.c_str());
 
  std::stringstream jsonStream;
  Writer::Write(objResponse, jsonStream);
  cout << jsonStream.str();
}
 
/**
 * @brief Executes a shell command
 * @param cmd The command to execute
 * @return The result of the execution
 *
 * This method executes a bash command using a pipe and a buffer to capture the
 * result of the execution, then return the execution to client code.
 */
string exec(char* cmd) {
    FILE* pipe = popen(cmd, "r");
    if (!pipe)
        return "ERROR";
 
    char buffer[128];
    string result = "";
    while(!feof(pipe)) {
        if(fgets(buffer, 128, pipe) != NULL)
            result += buffer;
    }
    pclose(pipe);
 
    // remove carriage return (new lines)
    result.erase(std::remove(result.begin(), result.end(), '\n'), result.end());
    return result;
}

Los resultados del programa muestran la resolución de la integral numérica para un número de particiones (n) dado, que son los rectángulos a resolver y sumar, entre mayor sea el número de rectángulos mayor será la exactitud del resultado, a continuación se muestran los resultados para n=5, n=50 y n=100.

N=5

N=50

N=100

Como se puede ver en los resultados anteriores el resultado converge rápidamente con el resultado real lo cual es bueno ya que se puede obtener una aproximación muy decente utilizando el poder de procesamiento y la fuerza bruta de la computadora..

Como punto final se anexa a continuación la implementación de dicha integral numérica donde se puede jugar con los parámetros de entrada para obtener un resultado deseado y su gráfica correspondiente.

[pageview url=”http://phoxonics.dyndns.org:5080/NumericalIntegral.htm” width=”100%” height=”100%” scrolling=”auto”]

¡Enjoy! 🙂

Leave a Reply

Your email address will not be published. Required fields are marked *