0

我用 C++ 为蒙特卡洛的积分编写了代码，当我使用二维积分时，它与我的其他函数配合得很好。我概括了 N 维积分的代码，在这种特殊情况下，我取 n = 10。我试图整合一个简单的函数 f = x1+x2+x3+x4+x5+....+x10 其中，x1。 ...x10 在限制 [-5.0, 5.0] 内。当我知道绝对结果应该为 0 时，我发现结果存在很大偏差。如果有人好心地查看我的代码并找出我的代码在哪里中断，我将不胜感激。我附上代码如下：

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;


//define multivariate function F(x1, x2, ...xk)            

double f(double x[], int n)
{
    double y;
    int j;
    y = 0.0;
    for (j = 0; j < n; j = j+1)
      {
         y = y + x[j];
      }     
    y = y;

    return y;
}
/*
 *  Function f(x1, x2, ... xk)
*/             

//define function for Monte Carlo Multidimensional integration

double int_mcnd(double(*fn)(double[],int),double a[], double b[], int n, int m)

{
    double r, x[n], p;
    int i, j;

    // initial seed value (use system time) 
    srand(time(NULL));  

    r = 0.0;
    p = 1.0;

    // step 1: calculate the common factor p
    for (j = 0; j < n; j = j+1)
      {
         p = p*(b[j]-a[j]);
      } 

    // step 2: integration
    for (i = 1; i <= m; i=i+1)
    {
        // calculate random x[] points
        for (j = 0; j < n; j = j+1)
        {
            x[j] = a[j] + static_cast <double> (rand()) /( static_cast <double> (RAND_MAX/(b[j]-a[j])));
        }         
        r = r + fn(x,n);
    }
    r = r*p/m;
    return r;
}


double f(double[], int);
double int_mcnd(double(*)(double[],int), double[], double[], int, int); 

int main(int argc, char **argv)
{    



    /* define how many integrals */
    const int n = 10;       

    double b[n] = {5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0,5.0};                    
    double a[n] = {-5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0, -5.0,-5.0};  

    double result;
    int i, m;
    int N = 20;

    cout.precision(6);
    cout.setf(ios::fixed | ios::showpoint); 

    // current time in seconds (begin calculations)
    time_t seconds_i;
    seconds_i = time (NULL);

    m = 2;                // initial number of intervals

    // convert command-line input to N = number of points
    //N = atoi( argv[1] );

    for (i=0; i <=N; i=i+1)
    {
        result = int_mcnd(f, a, b, n, m);
        cout << setw(30)  << m << setw(30) << result <<endl;
        m = m*2;
    }



// current time in seconds (end of calculations)
    time_t seconds_f;
    seconds_f = time (NULL);
    cout << endl << "total elapsed time = " << seconds_f - seconds_i << " seconds" << endl << endl;

    return 0;
}

我得到的输出是这样的：

                            4           -41046426010.339691
                             8           -14557222958.913620
                            16            25601187040.145161
                            32            29498213233.367203
                            64            -2422980618.248888
                           128           -13400105151.286720
                           256           -11237568021.855265
                           512            -5950177645.396674
                          1024            -4726707072.013641
                          2048            -1240029475.829825
                          4096             1890210492.995555
                          8192              573067706.448856
                         16384              356227781.143659
                         32768             -343198855.224271
                         65536              171823353.999405
                        131072             -143383711.461758
                        262144             -197599063.607231
                        524288              -59641584.846697
                       1048576               10130826.266767
                       2097152              100880200.681037

total elapsed time = 1 seconds

这与我的预期输出零相差无几。请帮我修复代码并提前感谢。

4

0 回答 0