c言語でのサンプル数分のｘ自乗の平均とｔの関係のグラフを求めたい。

Question

プログラムのフローチャートの通りにプログラムを組んだのですが、どこが間違っているのかいまいちわからなくて困っています。ran1.cとgasdev.cは正常にうごきます。langevin3.cのforの中がいまいちわかりません。
エラーは出ませんが、望んでいる結果と違う結果となっています。望んでいる結果は通常の目盛り(片対数や両対数ではない)で線形の関係になります。具体的には、[0, 1]の範囲でy=xを描きます。
追記
修正しましたが、実行の際に動作が停止してしまいました。

コード
//langevin3.cのプログラム。
#include <stdio.h>
#include <math.h>
#include "gasdev.c"

double gasdev(long *idum);

int main()
{
    FILE *fp;
    fp = fopen("langevin3.dat", "w");
    long idum;
    float F, ehta, gamma, kT, h, ave;
    int i, im, j, nsmp;
    float t, g, a, x2;
    float x[100];
    nsmp = 100;    
    idum = -123456789;
    F = 0;
    im = 1000;
    h = 1.0e-3;
    kT = 298 * 1.381e-23;
    ehta = 0.89e-3;
    gamma = 6 * acos(-1.0) * ehta * 0.5e-6;
    a = sqrt(2*kT*h/gamma);
    t = 0.0;
    for (i=1;i<=im;i++) {
        x[j] = 0.0;
        for (j=1;j<=nsmp;j++){
            x[j] += a*gasdev(&idum);
            if (i % 50 == 0) {
                x2 = 0.0;
                x2 += x[j]*x[j];                
                ave = x2/nsmp;
                fprintf(fp, "%f\t%e\n", t, ave * 1.0e12);
            } else {
                ;
            }
            t += h;            
        }         
    }
    fclose(fp);        
}

コード
//ran1.cのプログラム
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran1(long *idum)
{
    long j;
    long k;
    static long iy=0;
    static long iv[NTAB];
    double temp;

    if (*idum <= 0 || !iy) {
        if (-(*idum) < 1) *idum=1;
        else *idum = -(*idum);
        for (j=NTAB+7;j>=0;j--) {
            k=(*idum)/IQ;
            *idum=IA*(*idum-k*IQ)-IR*k;
            if (*idum < 0) *idum += IM;
            if (j < NTAB) iv[j] = *idum;
        }
        iy=iv[0];
    }
    k=(*idum)/IQ;
    *idum=IA*(*idum-k*IQ)-IR*k;
    if (*idum < 0) *idum += IM;
    j=iy/NDIV;
    iy=iv[j];
    iv[j] = *idum;
    if ((temp=AM*iy) > RNMX) return RNMX;
    else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX

コード
//gasdev.cのプログラム
#include <math.h>
#include "ran1.c"

double gasdev(long *idum)
{
    double ran1(long *idum);
    static int iset = 0;
    static double gset;
    double fac, rsq, v1, v2;

    if (iset == 0) {
        do {
            v1 = 2.0*ran1(idum) - 1.0;
            v2 = 2.0*ran1(idum) - 1.0;
            rsq = v1*v1 + v2*v2;
        } while (rsq >= 1.0 || rsq == 0.0);
        fac = sqrt(-2.0*log(rsq)/rsq);
        gset = v1*fac;
        iset = 1;
        return v2*fac;
    } else {
        iset = 0;
        return gset;
    }
}

Answer 1

float x[100];

初期値を設定していないため、x[1]以外は不定値です。

---

xをすべて0初期化したいなら以下のようにするだけで良いです。

float x[100] = { 0 };

Comments

[Okuhira007]

x[j] = 0.0;をjのforループの直後に追加しましたが、結果はうまくいきませんでした。
お忙しい中、ありがとうございます。

[e-watt]

そこで0を代入すると、次の行で積算していた部分が台無しになりませんか？
x[j] += a*gasdev(&idum);　// 0.0 にされた要素に加算（≒実質的に代入）

素直にループに入る前に初期化しましょう（独立にループを書く）

[Okuhira007]

了解しました。ありがとうございます。

[SHOMI]

for (i=1;i<=im;i++) {
x[j] = 0.0;
for (j=1;j<=nsmp;j++){
x[j]に代入後にjが変わるのでこれでは意味がないです。
jは1周目は不定値、2周目以降はnsmp+1になっています。
はじめにxの全要素を0初期化したいだけなら回答に追記した方法で良いです。

[Okuhira007]

解決しました！ありがとうございます。