Integración numérica con MathMore

Independiente del método de integración a utilizar, lo primero por hacer es definir la función a integrar. Siguiendo los delineamentos de la librería matemática MathMore, existen varias opciones para definir esta función. En los ejemplos siguientes, trabajaremos nuestra función como un objeto de C++. Las clases base de MathMore de las que se construyen las funciones definidas por el usuario son:

La función a integrar debe implementar una de estas clases abstractas. Necesitaremos en principio un solo archivo:

  1. Un .h -header file- o archivo de declaración de nuestra clase (opcional: incluiremos las definiciones de sus métodos si no hay un .cc)
  2. (Opcional): Un .cc en donde definiremos el método de evaluación de la función

La documentación para los métodos disponibles de integración en ROOT se encuentra en:

  • MathMore: la libreria MathMore ofrece una interface a la GNU Scientific Library o GSL.
  • GSL implementa algoritmos de integración numéricos estándar que se encuentran ampliamente documentados en libros.

Integración en 1D

El ejemplo siguiente consiste en la integración numérica de una función de 1D. Primero, hay que definir la función a integrar. Crear el archivo:

gedit dfx.h

En este archivo, escribiremos una clase que tiene como clase base IBaseFunctionOneDim

    1#include <Math/IFunction.h>
    2
    3class dfx :  public ROOT::Math::IBaseFunctionOneDim {
    4public: 
    5  
    6  double DoEval(const double xx) const {
    7
    8        double val = sin(xx) / xx;
    9        return val; 
   10  }
   11    
   12  ROOT::Math::IBaseFunctionOneDim* Clone() const
   13  {
   14    return new dfx();
   15  }
   16  
   17protected:
   18  
   19private:
   20  
   21};

  • En el archivo dfx.h hemos declarado y definido la clase dfx, la cual contiene nuestra función a integrar:

  • Como vemos en el código, esta clase tiene como base la clase IBaseFunctionOneDim y debe implementar obligatoriamente los siguientes métodos:

  1. DoEval( ): este método contiene la función a integrar, la evalúa en el punto que recibe como argumento.
  2. Clone( ): método que retorna un apuntador a una función nueva.

Eso seria todo por parte de la definición. Ahora necesitamos construir un objeto para realizar la integral. Ahora abrir un nuevo archivo en el que escribiremos la siguiente función:

gedit oneDintegral.C

  • Copiar en el el siguiente código:

    1#include "dfx.h"
    2#include "Math/GSLIntegrator.h"
    3
    4void oneDintegral() {
    5
    6  // Crear un apuntador a nuestra funcion
    7  dfx * f1 = new dfx();
    8  
    9  // Obtener un apuntador a un metodo de integracion de MC
   10  ROOT::Math::GSLIntegrator * nminteg = new ROOT::Math::GSLIntegrator( 1.E-6, 1.E-8, 1000 );
   11  
   12  // Pasar al integrador, la funcion
   13  nminteg->SetFunction( *(ROOT::Math::IGenFunction*)f1 );
   14 
   15  /// Limites de integracion 
   16  double xmin;
   17  double xmax;
   18  
   19  xmin =   0.0; // minimo variable x
   20  xmax =  5.0; // maximo variable x
   21
   22  // Calcular la integral
   23  double result = nminteg->Integral( xmin , xmax );
   24  
   25  // Imprimir en pantalla el resultado
   26
   27  std::cout << "integral result: " << result  << std::endl;
   28
   29  // Limpiar la memoria   
   30  delete f1;
   31  delete nminteg;
   32
   33}

  • Cada paso esta comentado y la lógica presente bastante clara. ROOT nos ofrece una interface al algoritmo de integración. Los algoritmos de integración pueden fijarse y ajustar según sea el caso. Para ello referirse a la documentación de la clase GSLIntegrator.

  • Para correr esta función, necesitamos iniciar root, cargar esta función y ejecutarla:

root -l
root [0] .L oneDintegral.C++
root [1] oneDintegral()
integral result: 1.54993

  • El resultado puede ser verificado a mano o en Mathematica si usted ha trabajado en un ejemplo mas difícil.

Integración en 2D

Es muy sencillo extender el metodo anterior para hacer integrales de mas de 1 dimension. Por ejemplo, si queremos integrar la funcion siguiente:

Basta con "anidar" una de las integrales mientras se mantiene una de las variables constante.

Por ejemplo, podemos crear el nuevo archivo:

gedit integration2d.h

En este archivo, definimos las dos clases:

  1. La primera contiene la función a integrar. Esta depende explícitamente de las dos variables.
  2. La segunda clase hace la integración con respecto a una de las variables mientras la otra se mantiene fija.

El código propuesto es el siguiente:

    1#include <Math/IFunction.h>
    2#include "Math/GSLIntegrator.h"
    3
    4// 1.0 Esta es la funcion a integrar:
    5
    6class dfxy :  public ROOT::Math::IBaseFunctionOneDim {
    7public: 
    8  
    9  double DoEval(const double yy) const {
   10    
   11    double val = (pow(xx, 2.0) + xx * pow(yy, 3.0) );
   12    return val; 
   13    
   14  }
   15  
   16  ROOT::Math::IBaseFunctionOneDim* Clone() const
   17  {
   18    return new dfxy();
   19  }
   20  
   21  void SetX( double x ) 
   22  {
   23    xx = x;
   24  }
   25  
   26private:
   27  
   28  double xx;
   29  
   30};
   31
   32// 2.0 Esta es la funcion que integra con respecto a una variable y deja la otra constante
   33
   34class dfx :  public ROOT::Math::IBaseFunctionOneDim {
   35public: 
   36  
   37  double DoEval(const double xx) const {
   38  
   39    //La funcion a integrar es:
   40    dfxy * f1 = new dfxy();
   41    
   42    //Integramos con respecto a Y dejando X constante
   43    f1->SetX( xx );
   44    
   45    //Ahora procedemos tal como en el ejemplo para 1D
   46    ROOT::Math::GSLIntegrator * nminteg = new ROOT::Math::GSLIntegrator( 1.E-6, 1.E-8, 1000 );
   47    
   48    // Pasar al integrador, la funcion
   49    nminteg->SetFunction( *(ROOT::Math::IGenFunction*)f1 );
   50    
   51    /// Limites de integracion (los limites no necesariamente deben estar fijos)
   52    double ymin = 1.0;
   53    double ymax = 2.0;
   54    
   55    // Calcular la integral
   56    double result = nminteg->Integral( ymin , ymax );
   57    
   58    return result; 
   59    
   60  }
   61  
   62  ROOT::Math::IBaseFunctionOneDim* Clone() const
   63  {
   64    return new dfx();
   65  }
   66  
   67};

Notar que la clase dfxy contiene un método para fijar la variable sobre la que no se hace la integración SetX( double x). El resto es cuestión de hace el llamado a dfx e integrar, tal como se hizo en el ejemplo para 1D. Fijando los limites en 0 <= x <= 1.0 y 1.0 <= y <= 2.0 obtenemos:

Integración con método de Monte Carlo

Integración númerica multidimensional. Este seria un ejemplo de una funcion de 2-D. Lo primero, es definir la funcion, para ello lo haremos en un archivo:

gedit d2xy.h

    1#include <Math/IFunction.h>
    2
    3class d2xy :  public ROOT::Math::IBaseFunctionMultiDim {
    4public: 
    5  
    6  double DoEval(const double * x) const {
    7
    8        double xx = x[0]; // variable 1
    9        double yy = x[1]; // variable 2
   10        double val = (xx*xx) * (yy*yy);
   11  
   12        return val; 
   13  }
   14  
   15  unsigned int NDim() const
   16  {
   17    return 2;
   18  }
   19  
   20  ROOT::Math::IBaseFunctionMultiDim* Clone() const
   21  {
   22    return new d2xy();
   23  }
   24  
   25protected:
   26  
   27private:
   28  
   29};

  • En el archivo d2xy.h hemos declarado y definido la clase d2xy, la cual contiene nuestra función a integrar:

  • Como vemos en el codigo, esta clase tiene como base la clase IBaseFunctionMultiDim y debe implementar obligatoriamente los siguientes métodos:

  1. DoEval( ): este método contiene la función a integrar, la evalúa según los argumentos que recibe a través de un arreglo de números y devuelve el valor.
  2. Clone( ): método que retorna un apuntador a una función nueva
  3. NDim( ): método que regresa la dimensión de la función, en este caso es 2.

  • Ahora podemos proceder a evaluar la integral. Por simplicidad trabajaremos esto desde una función. Abrir el siguiente archivo:

gedit mcintegral.C

  • Copiar en el el siguiente código:

    1#include "d2xy.h"
    2#include "Math/GSLMCIntegrator.h"
    3
    4void mcintegral() {
    5
    6  // Crear un apuntador a nuestra funcion
    7  d2xy * f1 = new d2xy();
    8  
    9  // Obtener un apuntador a un metodo de integracion de MC
   10  ROOT::Math::GSLMCIntegrator * nminteg = new ROOT::Math::GSLMCIntegrator( ROOT::Math::IntegrationMultiDim::kVEGAS, 1.E-6, 1.E-4, 500000 );
   11  
   12  // Pasar al integrador, la funcion
   13  nminteg->SetFunction( *(ROOT::Math::IMultiGenFunction*)f1 );
   14 
   15  /// Limites de integracion 
   16  double xmin[2];
   17  double xmax[2];
   18  
   19  xmin[0] = -1.0; // minimo variable 1
   20  xmax[0] =  1.0; // maximo variable 1
   21
   22  xmin[1] = -1.0; // minimo variable 2
   23  xmax[1] =  1.0; // maximo variable 2
   24  
   25  // Calcular la integral
   26  double result = nminteg->Integral( xmin , xmax );
   27  
   28  // Imprimir en pantalla el resultado
   29
   30  std::cout << "integral result: " << result << " " << nminteg->ChiSqr() << std::endl;
   31
   32  // Limpiar la memoria   
   33  delete f1;
   34  delete nminteg;
   35
   36}

  • Cada paso esta comentado y la lógica presente bastante clara. ROOT nos ofrece una interface al algoritmo de integracion por metodo de Monte Carlo VEGAS implementado en GSL mediante un objeto de clase GSLMCIntegrator.

  • Para correr esta función, necesitamos iniciar root, cargar esta función y ejecutarla:

root -l
root [0] .L mcintegral.C++
root [1] mcintegral()
integral result: 0.444446 0.230798

  • El resultado puede ser verificado a mano o en Mathematica si usted ha trabajado en un ejemplo mas difícil.

-- AndresOsorio - 27-Jun-2012

Edit | Attach | Watch | Print version | History: r6 < r5 < r4 < r3 < r2 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r6 - 2012-07-24 - AndresOsorio
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright &© 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback