第二次All-in啦。比上次积分更过了一点。
这次的题目跟上次的MonteCarlo算法有像似之处。

题目：
1）有100*100二维点矩阵。每一个元素的值是1或者-1（自旋方向。）
2）每一个元素的相互作用能是跟4个相邻元素的相互作用能之和。（存在周期性边界条件。S[0][0]左边的相邻元素是S[99][0]）
Sij跟一个相邻元素Sij+1的互作用能是-Sij*Sij+1.(反映两个自旋方向一样的时候能量最小，显示了元素的自旋方向的取向。)
3）初始条件是100*100个随机元素。
4）计算100万步。
5）每一步，随机选择一个元素。然后
   计算翻转点阵i的值，即Sij  ->  -Sij所带来的能量差deltaE，则点阵的翻转概率为
                1                     deltaE<=0
                e^(-con/temper)       ,deltaE>0     （con=40，temper=温度）
6）观察随时间（步骤）的所有元素的自旋和sigmaS的变化。
老师说，每一个温度，过一段时间后sigmaS会达到平衡值。
要求得出每一个温度跟平衡值的关系。
理想结果是低温的时候，平均值（sigmaS除以元素数）接近1.然后随温度的增加缓慢地下降
，过某一个温度之后，突然下降。（发生相变）
可我做的结果，不管温度再低，还是sigmaS在-200到200之间。平均值在0附近。
求大佬的敏锐的观察和正确的指正。谢谢！

  
/*  S140131022大佬来帮我的话，我很感谢了。*/   

return 50积分;     

下面是我写的代码。计算了温度100的时候的结果。
不用看全部代码，主要看下面的calculate的函数就行了。


#include "stdafx.h"
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#define CON 40 //概率比率系数     
#define E 2.718281828459045
#define STEP 10000000
using namespace std;

int S[100][100];
int randomS()   //返回1或者-1
{
    static bool first=true;
    int result;
    if (first)
    {
        first=false;
        srand((unsigned)time(NULL));
    }
    result=(int)rand()%2;
    if(result)
        return 1;
    else
        return -1;
}

double random(double range)  //返回给定范围内的任意数
{
    return rand()/(1.0+RAND_MAX)*range;
}

int random(int range)  //返回给定范围内的任意整数
{
    return  rand()%range;
}

int energy(int i,int j)    //i是横轴坐标，j是纵轴坐标,返回粒子的相互作用能
{
    if (i == 0)    //考虑周期性边界条件
    {
        if (j == 0)
            return -S[0][0] * S[0][1] - S[0][0] * S[1][0] - S[0][0] * S[0][99] - S[0][0]*S[99][0];
        if (j == 99)
            return -S[0][99] * S[0][98] - S[0][99] * S[1][99] - S[0][99] * S[99][99] - S[0][99] * S[0][0];
        else
            return -S[0][j] * S[0][j - 1] - S[0][j] * S[0][j + 1] - S[0][j] * S[1][j] - S[0][j] * S[99][j];
    }
    if (i == 99)
    {
        if (j == 0)
            return -S[99][0] * S[99][1] - S[99][0] * S[98][0] - S[99][0] * S[0][0] - S[99][0] * S[99][99];
        if (j == 99)
            return -S[99][99] * S[99][98] - S[99][99] * S[98][99] - S[99][99] * S[0][99] - S[99][99] * S[99][0];
        else
            return -S[99][j] * S[99][j - 1] - S[99][j] * S[99][j + 1] - S[99][j] * S[98][j] - S[99][j] * S[0][j];
    }
    if (j == 0)
        return  -S[i][0] * S[i - 1][0] - S[i][0] * S[i + 1][0] - S[i][0] * S[i][1] - S[i][0] * S[i][99];
    if (j == 99)
        return  -S[i][99] * S[i - 1][99] - S[i][99] * S[i + 1][99] - S[i][99] * S[i][98] - S[i][99] * S[i][0];
    else
        return -S[i][j]*S[i-1][j]-S[i][j]*S[i+1][j]-S[i][j]*S[i][j-1]-S[i][j]*S[i][j+1];
}

int reverse(int i)
{
    if (i==1)
        return -1;
    else
        return 1;
}

int energy_new(int i,int j)   //i是横轴坐标，j是纵轴坐标,返回粒子翻转后的相互作用能
{
    if (i == 0)    //考虑周期性边界条件
    {
        if (j == 0)
            return -reverse(-S[0][0]) * S[0][1] - reverse(S[0][0]) * S[1][0] - reverse(S[0][0]) * S[0][99] - reverse(S[0][0]) * S[99][0];
        if (j == 99)
            return -reverse(-S[0][99]) * S[0][98] - reverse(S[0][99]) * S[1][99] - reverse(S[0][99]) * S[99][99] - reverse(S[0][99]) * S[0][0];
        else
            return -reverse(-S[0][j]) * S[0][j - 1] - reverse(S[0][j]) * S[0][j + 1] - reverse(S[0][j]) * S[1][j] - reverse(S[0][j]) * S[99][j];
    }
    if (i == 99)
    {
        if (j == 0)
            return -reverse(S[99][0]) * S[99][1] - reverse(S[99][0]) * S[98][0] - reverse(S[99][0]) * S[0][0] - reverse(S[99][0]) * S[99][99];
        if (j == 99)
            return -reverse(S[99][99]) * S[99][98] - reverse(S[99][99]) * S[98][99] - reverse(S[99][99]) * S[0][99] - reverse(S[99][99]) * S[99][0];
        else
            return -reverse(S[99][j]) * S[99][j - 1] - reverse(S[99][j]) * S[99][j + 1] - reverse(S[99][j]) * S[98][j] - reverse(S[99][j]) * S[0][j];
    }
    if (j == 0)
        return  -reverse(S[i][0]) * S[i - 1][0] - reverse(S[i][0]) * S[i + 1][0] - reverse(S[i][0]) * S[i][1] - reverse(S[i][0]) * S[i][99];
    if (j == 99)
        return  -reverse(S[i][99]) * S[i - 1][99] - reverse(S[i][99]) * S[i + 1][99] - reverse(S[i][99]) * S[i][98] - reverse(S[i][99]) * S[i][0];
    else
        return -reverse(S[i][j])*S[i-1][j]-reverse(S[i][j])*S[i][i+1]-reverse(S[i][j])*S[i][j-1]-reverse(S[i][j])*S[i][j+1];
}

int deltaEnergy(int i,int j)   //返回能量差
{
    return energy_new(i,j)-energy(i,j);
}

int SUM()
{
    int sum=0;
    for (int i=0;i<100;i++)
    {
        for(int j=0;j<100;j++)
            sum+=S[i][j];
    }
    return sum;
}

double possibility(double con,double temper)
{
    double result;
    result=pow(E,-con/temper);
    return result;
}

void calculate(double temper)
{
    int i, j, step_i, tempS[100][100],count=0;
    int* sum=(int *)malloc(4 * 10000000 * sizeof(int));
    double con, p=1, deltaE;
    for (int i = 0; i < STEP; i++)
        sum[i] = 0;
    for (int i = 0; i < 100; i++)
        for (int j = 0; j < 100; j++)
        {
            S[i][j] = randomS();
            tempS[i][j] = S[i][j];
        }
    for (step_i=0;step_i<STEP;step_i++)
    {
        i = random(100);
        j = random(100);
        //cout << "Time:" << step_i << ", i:" << i << ", j" << j << ", Before change:" << S[i][j];
        deltaE = deltaEnergy(i,j);
        if (deltaE <= 0)
            tempS[i][j] = reverse(S[i][j]);
        else
        {
            con=(double)CON*deltaE;    //能量差包含在con里面
            p=possibility((double)con,temper);
            //cout << "Possibility: " << p;
            if (random(1)<p)   //以概率p
                tempS[i][j]=reverse(S[i][j]);
            else
                tempS[i][j]=S[i][j];
        }

        S[i][j]=tempS[i][j];

        sum[step_i]=SUM();
            //cout << ", After change:"   << S[i][j]  <<", Sum of Spin:" << sum[step_i] << ", Possib:" << p
            //<< ", deltaE:"   <<  deltaE   << ", con"   <<  con   << ", temper"   <<  temper   << endl;
        if (step_i%10000==0)
            cout<< "Time:" << step_i << "Sum:" << sum[step_i] << endl;
    }
}
int main()
{
    double temper;
    calculate(100);
    system("pause");
    return 0;
}
   

温度0.0001，con=4的时候运行，结果是sigmaS逐渐增加，增加到5000，7000.（一共10000个元素）
但是还是有问题。
1）为什么温度这么低才会sigmaS收敛于10000（或者-10000）
2）sigmaS的收敛速度非常慢，而且爬不到10000，中间停止增加或者又降低。

少说这些废话，分数有个屁用

其它我不考虑，只按要求将你的代码写规范

#include <vector>
#include <cmath>
#include <cstdlib>

struct S
{
public:
    S( size_t row, size_t col ) : col(col), matrix(row*col)
    {
        // 随机-1,+1
        for( size_t i=0; i!=matrix.size(); ++i )
            matrix[i] = rand()/(1.0+RAND_MAX) < 0.5 ? -1 : +1;
    }

    int DeltaEnergy_If_Flip( size_t r, size_t c ) const
    {
        const size_t row = matrix.size() / col;
        const int center = matrix[ r*col + c ];
        const int top = matrix[ (r+row-1)%row*col + c ];
        const int bottom = matrix[ (r+1)%row*col + c ];
        const int left = matrix[ r*col + (c+col-1)%col ];
        const int right = matrix[ r*col + (c+1)%col ];
        int Energy_old = 4 -  center * ( top + bottom + left + right );
        int Energy_New = 4 -  (-center) * ( top + bottom + left + right );
        return Energy_New-Energy_old;
    }

    int Sigma() const
    {
        int sum = 0;
        for( size_t i=0; i!=matrix.size(); ++i )
            sum += matrix[i];
        return sum;
    }

    void RandomFlip( double con, double temper )
    {
        const size_t row = matrix.size() / col;
        const size_t r = size_t(rand()/(1.0+RAND_MAX) * row);
        const size_t c = size_t(rand()/(1.0+RAND_MAX) * col);

        if( DeltaEnergy_If_Flip(r,c)<=0 || rand()/(1.0+RAND_MAX)<exp(-con/temper) )
            matrix[ r*col + c ] = -matrix[ r*col + c ];
    }

private:
    size_t col;
    std::vector<int> matrix;
};

#include <iostream>
#include <ctime>
using namespace std;

void calculate( double con, double temper )
{
    S s( 100, 100 );
    for( size_t i=0; i!=10000000; ++i )
    {
        s.RandomFlip( con, temper );

        if( i%10000 == 0 )
            cout << "Time:" << i << "Sum:" << s.Sigma() << '\n';
    }
}

int main( void )
{
    srand( (unsigned)time(NULL) );
    calculate( 40, 100 );
}