2014-04-03 5 views
2

부스트 odeint를 사용하여 다차원 적분을 높은 정확도로 계산할 때 권장되는 방법은 무엇입니까? 다음 코드는 -1 내지 2 (F) = X * Y를 통합되지만 분석 용액에 에러 대하여 1 % 이상 (GCC 4.8.2, -std = C++ 0X)이다부스트 odeint를 사용하는 정확한 다차원 적분

#include "array" 
    #include "boost/numeric/odeint.hpp" 
    #include "iostream" 
    using integral_type = std::array<double, 1>; 
    int main() { 
     integral_type outer_integral{0}; 

     double current_x = 0; 
     boost::numeric::odeint::integrate(
     [&](
      const integral_type&, 
      integral_type& dfdx, 
      const double x 
     ) { 
      integral_type inner_integral{0}; 

      boost::numeric::odeint::integrate(
      [&current_x](
       const integral_type&, 
       integral_type& dfdy, 
       const double y 
      ) { 
       dfdy[0] = current_x * y; 
      }, 
      inner_integral, 
      -1.0, 
      2.0, 
      1e-3 
     ); 

      dfdx[0] = inner_integral[0]; 
     }, 
     outer_integral, 
     -1.0, 
     2.0, 
     1e-3, 
     [&current_x](const integral_type&, const double x) { 
      current_x = x; // update x in inner integrator 
     } 
    ); 
     std::cout 
     << "Exact: 2.25, numerical: " 
     << outer_integral[0] 
     << std::endl; 
     return 0; 
    } 

인쇄 :

Exact: 2.25, numerical: 2.19088 

내부 통합에서보다 엄격한 정지 조건을 사용해야합니까, 아니면이를 더 빠르고 정확하게 수행 할 수 있습니까? 감사!

답변

2

우선, ODE 체계 (http://en.wikipedia.org/wiki/Numerical_integration)보다 고차원 적분을 계산하는 더 나은 수치 적 방법이있을 수 있지만, 이것이 odeint의 깔끔한 응용이라고 생각합니다.

그러나 코드의 문제는 내부 통합에 대한 x 값을 업데이트하기 위해 외부 통합에서 관찰자를 사용한다고 가정한다는 것입니다. 그러나 integrate 함수는 조밀 한 출력 스테퍼를 내부적으로 사용합니다. 즉, 실제 시간 간격과 관찰자 호출이 동기화되어 있지 않음을 의미합니다. 따라서 내부 통합을위한 x는 적절한 순간에 업데이트되지 않습니다. 빠른 수정은 runge_kutta4 스테퍼와 함께 integrate_const를 사용하는 것입니다.이 스테퍼는 일정한 단계 크기를 사용하고 외부 루프의 각 단계 후에 옵저버 호출 및 x 업데이트가 호출되도록합니다. 그러나 이것은 통합 루틴의 내부 세부 사항에 의존하는 약간의 해킹입니다. 더 나은 방법은 상태가 실제로 2 차원 적이지만 각 통합이 두 변수 중 하나에서만 작동하는 방식으로 프로그램을 설계하는 것입니다.

+0

감사합니다. current_x가 의도 한대로 업데이트되지 않을 수도 있다고 생각하지 않았습니다. 나는 다차원 통합에 더 적합한 라이브러리를 찾았습니다 : http://ab-initio.mit.edu/wiki/index.php/Cubature (GPL). 이 경우 최소한 9 개의 샘플만으로도 훨씬 작은 오차를 줄 수 있습니다. – user3493721