Monte Carlo methods in C++

Asked

Viewed 470 times

3

This is the code I created to calculate in C++ using the Monte Carlo methods, the average distance between the points evenly distributed in a cube. I don’t understand why in the output file "dista_media_quando_se_varia_o_numero_de_configuracoes.txt", the value of the second column is always the same. I wanted the value to vary until it stabilized in a value. I put the formulas used in an image attached. The goal is that with the txt file you can create a graph like the one in attachment where in the x-axis is the logarithm of the values on the left of the txt file and the y-axis is the value of the right column of the txt file. How can I change the file to get the value the mean distance over time (right column of txt file)?

inserir a descrição da imagem aqui inserir a descrição da imagem aqui

#include<iostream>
#include<cstdlib>
#include<math.h>
#include<fstream>
#include<iomanip>

using namespace std;

void distanciamediatotal(int &numero_configuracoes, int &numero_de_pontos, const char* name_of_file, int &semente, int caso){
        double x[numero_de_pontos];
        double y[numero_de_pontos];
        double z[numero_de_pontos];
        double dmedia_auxiliar;
        double dij = 0.0;
        double dijsomatorio_quando_i_menor_j = 0.0;
        double distancia_media_total = 0.0;
        double somatorio_d_media = 0.0;
        long long int i;
        long long int j;
        long long int k;
        long long int m;

        ofstream file;
        file.open(name_of_file);

        switch(caso){
                case 0:
                        for(k=0;k<numero_configuracoes;k++){
                                dijsomatorio_quando_i_menor_j=0.0;
                                distancia_media_total = 0.0;
                                srand48(semente);

                                for(i=0;i<numero_de_pontos;i++){
                                        x[i]=drand48();
                                        y[i]=drand48();
                                        z[i]=drand48();
                                }

                        for(j=1;j<numero_de_pontos;j++){
                                for (i=0; i<j;i++){

                                        dij=sqrt(pow((x[j]-x[i]),2)+pow((y[j]-y[i]),2)+pow((z[j]-z[i]),2));

                                        dijsomatorio_quando_i_menor_j+=dij;
                                }
                        }


           dmedia_auxiliar=dijsomatorio_quando_i_menor_j*2/(double)(numero_de_pontos*(numero_de_pontos-1));
                                somatorio_d_media+= dmedia_auxiliar;
                                distancia_media_total = somatorio_d_media/(k+1);

                                file << k+1 << "\t" << setprecision(10) << distancia_media_total << endl;
                }
                file.close();
                break;
                    case 1:
                        for(m=2;m<=numero_de_pontos;m++){
                                somatorio_d_media = 0.0;

                                for (k=0;k<numero_configuracoes;k++){
                                        srand48(semente);
                                        dijsomatorio_quando_i_menor_j=0.0;
                                        distancia_media_total = 0.0;

                                        for(i=0;i<numero_de_pontos;i++){
                                                x[i]=drand48();
                                                y[i]=drand48();
                                                z[i]=drand48();
                                        }
                                        for(j=1;j<m;j++){
                                                for (i=0; i<j;i++){

dij=sqrt(pow((x[j]-x[i]),2)+pow((y[j]-y[i]),2)+pow((z[j]-z[i]),2));
                                                    dijsomatorio_quando_i_menor_j+=dij;
                                            }
                                    }
                                    dmedia_auxiliar=dijsomatorio_quando_i_menor_j*2/(double)(m*(m-1));
                                    somatorio_d_media+= dmedia_auxiliar;
                            }
                            distancia_media_total = somatorio_d_media/numero_configuracoes;
                            file << m << "\t" << setprecision(10) << distancia_media_total << endl;
                    }
                    file.close();
            break;
                            default:
                            cout << "Teste nao identificado \n" << endl;
    }
}

int main (){

    long long int numero=1;
    long long int numerototal=10;
    long long int Nd;
    long long int i;
    int semente=17;
    srand48(semente);
    double x=0.0;
    double y=0.0;
    double distancia;
    double pi=0.0;
    double errocomlogaritmo=0.0;

    // Exercicio 1 - Calculo do valor de pi

    //Gerar pontos distribuidos uniformemente num quadrado unitario;
    //Calculo da distancia de cada ponto a origem;
    //Determinacao dos pontos que se encontram dentro do semicirculo
    //Estimar o valor de pi;
    //Fazer dois graficos: o primeiro do valor da estimativa de pi
    //em funcao do numero de pontos gerados; o segundo da dependencia
    //funcional em funcao do numero de pontos gerados.


    /*ofstream file;
    file.open("Dados_para_o_numero_pi.txt");

    for(int k=0;k<numerototal;k++){
            numero*=10;
            Nd=0;
            for(i=0;i<numero;i++){
                    x = drand48();
                    y = drand48();
                    distancia = pow(x,2)+pow(y,2);
                    //Se os pontos que estão no quadrado estao tambem dentro do
                    //semicirculo, soma mais uma unidade a Nd
                    if(distancia<=1.0)
                       Nd++;
            }

            pi = 4*Nd/(double)numero;
            errocomlogaritmo = log10(fabs(M_PI-pi));


            file << numero << "\t" << setprecision(10) << pi << "\t" << k+1 << "\t" << setprecision(8) <<
            errocomlogaritmo << "\t" << endl;
    }
    file.close();*/

    //Exercicio 2 - Pontos numa caixa

    //Resolucao do integral usando o Metodo de Monte Carlo

    //varia-se o numero de configuracoes, mantendo fixo o numero de pontos
    int numero_configuracoes_maximo = 100000;
    int numero_de_pontos=1000;
    distanciamediatotal(numero_configuracoes_maximo, numero_de_pontos, "distancia_media_quando_se_varia_o_numero_de_configuracoes.txt", semente, 0);

    //Varia-se o numero de pontos, mantendo fixo o numero de configuracoes
    int numero_configuracoes = 1000;
    //distanciamediatotal(numero_configuracoes, numero_de_pontos, "distancia_media_quando_se_varia_o_numero_de_pontos.txt", semente, 1);

    return 0;
}
No answers

Browser other questions tagged

You are not signed in. Login or sign up in order to post.