时间戳

首页 > 解决方案 > 蒙特卡洛积分返回不正确的值，覆盖内存的问题

c - 蒙特卡洛积分返回不正确的值，覆盖内存的问题

问题描述

我在计算圆环质心的积分时遇到问题，它应该返回 (2.4076, 0.16210, 0.0)。

该程序适用于 pi/4 的估计，但是我认为当我尝试使用 setRandomDomain() 函数覆盖现有点时存在问题。

这是我的代码：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define DIM 1000000

double random_double() {
    static const int a = 16807;
    static const int c = 0;
    static const long long m = 2147483647;
    static long long seed = 1;
    seed = (a * seed + c) % m;
    return ((double) seed) / m;
}

typedef struct Domain_t {
    double *x;
    double *y;
    double *z;
} Domain;

void constructDomain(Domain (**p_domain)) {
    *p_domain = malloc(sizeof(Domain));
    if(p_domain == NULL) {
        printf("ERROR: Memory allocation failed\n");
    }
    (*p_domain)->x = malloc(DIM * sizeof(double));
    if ((*p_domain)->x == NULL) {
        printf("ERROR: Memory allocation failed\n");
    }
    (*p_domain)->y = malloc(DIM * sizeof(double));
    if ((*p_domain)->y == NULL) {
        printf("ERROR: Memory allocation failed\n");
    }
    (*p_domain)->z = malloc(DIM * sizeof(double));
    if((*p_domain)->z == NULL) {
        printf("ERROR: Memory allocation failed\n");
    }
}

void delDomain (Domain (**p_domain)) {
    if (p_domain != NULL) {
        free ((*p_domain)->z);
        free ((*p_domain)->y);
        free ((*p_domain)->x);
        free (*p_domain);
    }
}

double radiusFunc(double point_x, double point_y) {
    return sqrt(pow(point_x,2)+pow(point_y,2));
}

double G(double point_x, double point_y, double point_z, int R) {
    return pow(point_z,2)+pow(radiusFunc(point_x,point_y)-(double)R,2);
}

typedef struct Volume_t {
    int R;
    int r;
    int lower_x;
    int upper_x;
    int lower_y;
    int upper_y;
    int lower_z;
    int upper_z;
    int V;
} Volume;

void setVolume(Volume (*p_volume), int R, int r, int x1, int x2, int y1, int y2, int z1, int z2) {
    p_volume->R = R;
    p_volume->r = r;
    p_volume->lower_x = x1;
    p_volume->upper_x = x2;
    p_volume->lower_y = y1;
    p_volume->upper_y = y2;
    p_volume->lower_z = z1;
    p_volume->upper_z = z2;
    if(z1 == 0 && z2 == 0)
        p_volume->V = (x2-x1)*(y2-y1);
    else if(y1 == 0 && y2 == 0)
        p_volume->V = (x2-x1)*(z2-z1);
    else if(x1 == 0 && x2 == 0)
        p_volume->V = (y2-y1)*(z2-z1);
    else
        p_volume->V = (x2-x1)*(y2-y1)*(z2-z1);
}

void setInitialDomain(Domain (**p_domain)) {
    int i;
    for(i=0;i<DIM;i++) {
        (*p_domain)->x[i] = random_double();
        (*p_domain)->y[i] = random_double();
        (*p_domain)->z[i] = random_double();
    }
}

void setRandomDomain(Domain (*p_domain), Domain (**p_new_domain), Volume (*p_volume)) { 
    int i;
    for(i=0;i<DIM;i++) {
        (*p_new_domain)->x[i] = p_domain->x[i]*(double)(p_volume->upper_x - p_volume->lower_x) + (double)p_volume->lower_x;
        (*p_new_domain)->y[i] = p_domain->y[i]*(double)(p_volume->upper_y - p_volume->lower_y) + (double)p_volume->lower_y;
        (*p_new_domain)->z[i] = p_domain->z[i]*(double)(p_volume->upper_z - p_volume->lower_z) + (double)p_volume->lower_z;
    }
}

double setIntegrand(Domain (*p_domain), char c) {
    double *p_x = p_domain->x;
    double *p_y = p_domain->y;
    double *p_z = p_domain->z;
    if(c=='x')
        return *p_x;
    else if(c=='y')
        return *p_y;
    else if(c=='z')
        return *p_z;
    else
        return 1.;
}

double calculateIntegral(Domain (*p_domain), Volume (*p_volume), char c) {
    int i;
    double F = 0.;
    for(i=0;i<DIM;i++) {
        if(G(p_domain->x[i], p_domain->y[i], p_domain->z[i], p_volume->R)<=(double)p_volume->r) {
            F += setIntegrand(p_domain, c);
        }
    }
    return F*(double)p_volume->V/(double)DIM;
}

int main() {
    Domain *p_initial_domain;
    Domain *p_random_domain;

    constructDomain(&p_initial_domain);
    printf("Point 1: successful\n");
    constructDomain(&p_random_domain);
    printf("Point 2: successful\n");

    setInitialDomain(&p_initial_domain);

    Volume circle, *p_circle;
    p_circle = &circle;
    setVolume(p_circle,0,1,0,1,0,1,0,0);
    setRandomDomain(p_initial_domain, &p_random_domain, p_circle);

    printf("PI/4 is approximately %f\n", calculateIntegral(p_random_domain, p_circle, 'p'));

    Volume torus, *p_torus;
    p_torus = &torus;
    setVolume(p_torus,3,1,1,4,-3,4,-1,1);
    setRandomDomain(p_initial_domain, &p_random_domain, p_torus);

    double M = calculateIntegral(p_random_domain, p_torus, 'p');
    double X = calculateIntegral(p_random_domain, p_torus, 'x');
    double Y = calculateIntegral(p_random_domain, p_torus, 'y');
    double Z = calculateIntegral(p_random_domain, p_torus, 'z');
    printf("rho integral is approximately %f\n", M);
    printf("x integral is approximately %f\n", X);
    printf("y integral is approximately %f\n", Y);
    printf("z integral is approximately %f\n", Z);
    printf("Centre of mass is approximately (%f, %f, %f)\n", X/M, Y/M, Z/M);

    delDomain(&p_initial_domain);
    delDomain(&p_random_domain);

// return pointers??
// array of structs??

    return 0;
}

目前输出：

PI/4 is approximately 0.785436
rho integral is approximately 22.101282
x integral is approximately 22.101801
y integral is approximately -45.953770
z integral is approximately 11.298411
Centre of mass is approximately (1.000023, -2.079235, 0.511211)

任何想法如何解决这个问题？

另外，请有人解释一下我将如何使用返回指针的函数以及为什么创建一个结构数组而不是数组结构可能更好？

标签： cmontecarlo

解决方案

你的问题是你setIntegrand在所有点上循环调用，但你总是拿第一点：

double *p_x = p_domain->x;

// ...
return *p_x;

这将返回double数组中的第一个。请记住，这*x相当于x[0]. 将索引传递给函数：

double setIntegrand(Domain (*p_domain), char c, int i)
{
    if (c == 'x') return p_domain->x[i];
    if (c == 'y') return p_domain->y[i];
    if (c == 'z') return p_domain->z[i];
    return 1.;
}

然后用那个索引调用它。

for (i = 0; i < DIM; i++) {
    if (G(...) <= p_volume->r) {
        F += setIntegrand(p_domain, c, i);
    }
}

至于您的其他问题：使用结构数组可以使在一起的东西（这里是点的三个坐标）保持在附近。您也可以轻松地将一个点传递给具有单个参数的函数。

如果你有一个构造函数，它是一个通过在堆上分配并初始化新内存来创建新事物的函数，返回一个指针是一个有用的习惯用法。（我发现它比将点传递给指针更惯用，但设计fopen_s函数的人并不这么认为。）

让我们把这两个变化放在一起：

typedef struct Point Point;
typedef struct Domain Domain;

struct Point {
    double x;
    double y;
    double z;
};

struct Domain {
    size_t length;
    Point *point;
};

Domain *constructDomain(size_t length)
{
    Domain *dom = malloc(sizeof(*dom));

    if (dom) {
        dom->length = length;
        dom->point = calloc(length, sizeof(*dom->point));

        // check success
    }

    return dom;
}

推荐阅读