| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 2144 人关注过本帖
标题:求大佬! MonteCarlo算法
只看楼主 加入收藏
复旦
Rank: 3Rank: 3
等 级:论坛游侠
威 望:2
帖 子:81
专家分:124
注 册:2018-10-29
结帖率:100%
收藏
已结贴  问题点数:47 回复次数:6 
求大佬! MonteCarlo算法
有一个程序利用MonteCarlo行走方式,模拟势阱里的粒子
二维势阱的半径为40
粒子的行走方式:
    1.每一次只走一步
    2.每一次移动距离为固定值s
    3.每一次行走的方向随机
    4.行走概率:如果由行走所带来的能量差为 deltaE
        则这种改变的发生的概率p为:
        e^(-deltaE*con/T))  ,deltaE>=0
        1                   , deltaE<0

这里E是对距离(从原点到粒子)的函数
E=-0.5cos(0.5pi*d)-cos(1.5*pi*d)      ,d:距离,pi:圆周率
然后设置步骤数(1000),统计这个粒子处于基态的概率(基态是从原点的距离0-14的区域)
最后的目的是得出不同温度下的概率。(对每一个温度做统计计算)

这个问题我大概就结了一个星期。(我太菜了)
写出了代码,但是结果不太理想。 有大佬了解这样的问题吗? 我的代码哪里有问题?
老师说,正确的结果是 概率随温度慢慢下降然后过某一个点,突然下降到0附近(发生相变)。

我的代码是这样的。不到100行。
我拿出我的全部积分奖励大佬。

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;


double dist(double x,double y)
{
    double result;
    result=sqrt(x*x+y*y);
    return result;
}

double energy(double x,double y)
{
    double result;
    result=-0.5*cos(0.5*3.141592*dist(x,y)-cos(1.5*3.141592*dist(x,y)));
    return result;
}

double possibility(double con,double temper)
{
    return pow(2.71828,-con/temper);
}

double calculate (double temper)
{
    int i,step,con=200;          //概率比例系数设为40
    double x[10000],y[10000],tempAngle,tempDist,tempX,tempY,p,count=0;
    step=1000;            //统计一万步的结果
    srand((unsigned)time(NULL));
    tempAngle=rand()%1000;
    tempDist=rand()%40;   //设势阱的半径为步伐的40倍
    x[0]=tempDist*cos((double)tempAngle);
    y[0]=tempDist*sin((double)tempAngle);
    for (i=1;i<step;i++)
    {
        tempAngle=rand()%1000;
        tempX=x[i-1]+cos(tempAngle);
        tempY=y[i-1]+sin(tempAngle);
        if (dist(tempX,tempY)>=40)   //如果超过了边界
        {
            x[i]=x[i-1];
            y[i]=y[i-1];
        }
        else
        {
            if (energy(tempX,tempY)<energy(x[i-1],y[i-1]))   //如果预测位置的能量比原来位置的能量小
            {
                x[i]=tempX;
                y[i]=tempY;
            }
            else
            {
                p=possibility(con,temper);
                p*=1000;
                if ((rand()%1000)<p)   //以概率p
                {
                    x[i]=tempX;
                    y[i]=tempY;
                }
                else
                {
                    x[i]=x[i-1];
                    y[i]=y[i-1];
                }
            }
        }
    }
    for (i=0;i<step;i++)
    {
        if (dist(x[i],y[i])<14)
            count++;
    }

    return count/(double)step;
}

int main()
{
    double temper;
    int i;
    for (temper=100;temper<1000;temper+=10)
    {
        cout << "Temperature:" << temper << ' ' << "Possibility"
        << calculate(temper) << endl;
    }
    scanf("%d",&i);
    return 0;
}


搜索更多相关主题的帖子: 概率 include double result step 
2018-11-07 09:53
rjsp
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:528
帖 子:9007
专家分:53942
注 册:2011-1-18
收藏
得分:47 
连蒙带猜,问题很多啦,只说几个主要的
1. possibility中,你求的是 e^(-con/T),可是要求中写的是 e^(-deltaE*con/T)
2. tempAngle=rand()%1000; 中tempAngle的单位是什么?角度还是弧度?cos函数的参数应当是弧度。%1000反正很奇怪,既不是360的整数倍,又不是PI的整数倍。
顺便问一下,从#include "stdafx.h"看出你用的是VC,但不知道你用的是哪个版本的VC,是VC6,还是VC2017?因为我怕我写的标准C++代码,你用VC等不标准的编译器编译失败。
2018-11-07 10:52
rjsp
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:528
帖 子:9007
专家分:53942
注 册:2011-1-18
收藏
得分:0 
能不能在40半径内出现你老师说的“突然下降到0附近”,与温度及con相关吧,结果还与一开始粒子所处位置相关吧
我是一窍不通,只根据你的要求写
程序代码:
#define _USE_MATH_DEFINES
#include <cmath>
#include <random>
#include <cstdlib>
#include <ctime>
#include <iostream>

double dist( double x, double y )
{
    return hypot( x, y ); // 如果你的VC不认识 hypot 就用 _hypot
}

double random( double range )
{
    static bool bfirst = true;
    if( bfirst )
    {
        bfirst = false;
        srand( unsigned(time(NULL)) );
    }
    return rand()/(1.0+RAND_MAX) * range;
}

double energy( double x, double y )
{
    double d = dist( x, y );
    return -0.5*cos(0.5*M_PI*d)-cos(1.5*M_PI*d); // 假如你的公式中用的是角度,那么这里得先乘以M_PI/180
}

double possibility( double x_from, double y_from, double x_to, double y_to, double con, double temper )
{
    double e_from = energy( x_from, y_from );
    double e_to = energy( x_to, y_to );
    return e_from>e_to ? 1 : exp((e_from-e_to)*con/temper);
}

double calculate( double temper )
{
    const int con = 200;
    const double radius = 40; // 势阱半径40
    const unsigned step = 10000; // 模拟一万次
    unsigned count = 0; // 基态数量

    double x, y; // 当前位置
    {
        double a = random( 2*M_PI ); // 随机角度
        double r = random( radius ); // 随机半径
        x = r * cos( a );
        y = r * sin( a );
    }
    for( unsigned i=0; i!=step; ++i )
    {
        double a = random( 2*M_PI ); // 随机角度
        double x_new = x + cos( a ); // 新位置
        double y_new = y + sin( a );

        if( dist(x_new,y_new) < radius ) // 不越界
        {
            double p = possibility( x,y, x_new,y_new, con, temper ); // 从(x,y)到(x_new,y_new)的跳跃几率
            if( random(1) < p ) //
            {
                x = x_new;
                y = y_new;
            }
        }

        if( dist(x,y) < 14 )
            ++count;
    }

    return count*1.0/step;
}

using namespace std;

int main( void )
{
    for( int temper=100; temper<1000; temper+=10 )
    {
        cout << "Temperature: " << temper
             << " Possibility: " << calculate(temper) << '\n';
    }
}

2018-11-07 12:11
复旦
Rank: 3Rank: 3
等 级:论坛游侠
威 望:2
帖 子:81
专家分:124
注 册:2018-10-29
收藏
得分:0 
原来这里有错! possobility里没写了deltaE。
非常感谢!
然后tempAngle的单位是radian。啊!原来这里也不对。受到的是随机整数。
我用的是vs2010里的c++控制台程序。保存一次后用vs2017打开的。 因为vs2017里没有c++。
我晚上看代码了。
很头疼的问题找到了问题点了。
感谢大佬!
2018-11-07 15:49
rjsp
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:528
帖 子:9007
专家分:53942
注 册:2011-1-18
收藏
得分:0 
以下是引用复旦在2018-11-7 15:49:39的发言:

原来这里有错! possobility里没写了deltaE。
非常感谢!
然后tempAngle的单位是radian。啊!原来这里也不对。受到的是随机整数。
我用的是vs2010里的c++控制台程序。保存一次后用vs2017打开的。 因为vs2017里没有c++。
我晚上看代码了。
很头疼的问题找到了问题点了。
感谢大佬!
如果要我说的话,
E=-0.5cos(0.5pi*d)-cos(1.5*pi*d)
就不对,这是一个周期震荡函数,震荡周期大约在 4 左右。而理论上它应该是d越大则E越大的函数。
因为这个E公式不对,所以接下来的几率也就不对,也就没有粒子趋向于圆心的行为,导致最终结果完全是随机的。

“vs2010里的c++控制台程序”
------ VS 包含 VC++、、VC# 等等,所以不要说你用的是VS,M$说过“VS不是编译器、VC才是”

“vs2017里没有c++”
------ 是没有C++,但有VC++2017
2018-11-07 16:56
复旦
Rank: 3Rank: 3
等 级:论坛游侠
威 望:2
帖 子:81
专家分:124
注 册:2018-10-29
收藏
得分:0 
E函数长得这样。
图片附件: 游客没有浏览图片的权限,请 登录注册

温度上升时,跃迁到更高的能量的地方的概率接近于1,粒子的行为越来越接近随机行走。
所以处于基态的统计概率随温度上升,逐渐下降(过某一个点突然下降),下降到某个值。
(基态的面积除以势阱的面积,比这个值稍微高一丢丢。)
大概是这么个思路。
2018-11-07 19:40
复旦
Rank: 3Rank: 3
等 级:论坛游侠
威 望:2
帖 子:81
专家分:124
注 册:2018-10-29
收藏
得分:0 
学到了您的random函数用法。
代码主要我自己修改了。结果虽然不能说完美,但是已经很不错了。
谢谢!我要结贴。 (以后的日子。。。)
2018-11-07 23:28
快速回复:求大佬! MonteCarlo算法
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.020316 second(s), 10 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved