processingの備忘録　-ビュフォンの針- - プログラミングの備忘録

f:id:taq2777:20211007094347p:plain

こんにちは。
今回は「ビュフォンの針」を取り上げます。

ビュフォンの針は、一定間隔$d$で引かれた線の上に長さ$l$の針を$n$本落としたとき、線と針が交わる確率$p$は、 $$ p = \frac{2l}{\pi d} $$ と表されるというものです。
（確率の問題は受験時も避けてきたものなので詳しい話はしません・・・できません・・・）

さて、これをprocessingで動かしたいとなったとき、まずやってみるのは線を引くことですね。
一回描いてくれればそれでいいので、諸々の処理はsetupの中に書くことにします。
一定間隔、例えば100ピクセル間隔で線を引くプログラムを書くと、

float d;

void setup(){
    size(500, 500);
  
    d = 100;
  
    background(0);
  
    for(int i = 0; i < width/d; i ++){ //線を描く
        stroke(255);
        line(0, d*i, width, d*i);
    }
}

続いて、針を描くプログラムを追加します。
ここで、針の長さ$l$は$\displaystyle l = \frac{d}{2}$とすれば$\displaystyle p = \frac{1}{\pi}$と単純化できるため、そうします。
針はline関数で描きますが、始点・終点の座標や針の角度はランダムな値として配列に保存しておいて使っています。

float d, l;
int n = 100;
float[] x1 = new float[n];
float[] y1 = new float[n];
float[] x2 = new float[n];
float[] y2 = new float[n];
float[] rad = new float[n];

void setup(){
    size(500, 500);
  
    d = 100;
    l = d/2;
    for(int i = 0; i < n; i ++){
          x1[i] = random(l, width-l);
          y1[i] = random(l, height-l);
          rad[i] = random(PI);
          x2[i] = x1[i] + l*cos(rad[i]);
          y2[i] = y1[i] + l*sin(rad[i]);
    }
  
    background(0);
  
    for(int i = 0; i < width/d; i ++){ //線を描く
        stroke(255);
        line(0, d*i, width, d*i);
    }
  
    for(int i = 0; i < n; i ++){ //針を描く
        line(x1[i], y1[i], x2[i], y2[i]);
    }
}

次に、線と針が交わっているかどうかを判断する処理を書きます。ここで、針の始点は終点よりもy座標が小さくなるようにしてあるため、始点が線の手前（y座標が小さい）にあり、終点が線の奥（y座標が大きい）にあるという条件をつければ、交わっていると言えます。

float d, l;
int n = 100;
float nc = 0; //交わっている針の本数をncとし、初期値は0
float[] x1 = new float[n];
float[] y1 = new float[n];
float[] x2 = new float[n];
float[] y2 = new float[n];
float[] rad = new float[n];

void setup(){
    size(500, 500);
  
    d = 100;
    l = d/2;
    for(int i = 0; i < n; i ++){
          x1[i] = random(l, width-l);
          y1[i] = random(l, height-l);
          rad[i] = random(PI);
          x2[i] = x1[i] + l*cos(rad[i]);
          y2[i] = y1[i] + l*sin(rad[i]);
    }
  
    background(0);
  
    for(int i = 0; i < width/d; i ++){ //線を描く
        stroke(255);
        line(0, d*i, width, d*i);
    }
  
    for(int i = 0; i < n; i ++){ //針を描く
        for(int j = 0; j < width/d; j ++){
            if(y1[i] <= d*j && y2[i] >= d*j){ //交わっているとき
                line(x1[i], y1[i], x2[i], y2[i]);
                nc ++; //交わっている針の数を+1する
            }                
            if(y1[i] > d*j && y2[i] < d*(j+1)){ //交わっていないとき
                line(x1[i], y1[i], x2[i], y2[i]);
            }
        }
    }
    
    println("crossed lines = " + round(nc)); //交わっている針の本数を表示
}

最後に、確率から$\pi$を計算します。
交わっている針の本数を$n_c$とすると、交わる確率$p$は$\displaystyle p = \frac{n_c}{n}$なので、$\displaystyle \pi = \frac{n}{n_c}$で計算できます。

float d, l;
int n = 100;
float nc = 0; //交わっている針の本数をncとし、初期値は0
float[] x1 = new float[n];
float[] y1 = new float[n];
float[] x2 = new float[n];
float[] y2 = new float[n];
float[] rad = new float[n];

void setup(){
    size(500, 500);
  
    d = 100;
    l = d/2;
    for(int i = 0; i < n; i ++){
          x1[i] = random(l, width-l);
          y1[i] = random(l, height-l);
          rad[i] = random(PI);
          x2[i] = x1[i] + l*cos(rad[i]);
          y2[i] = y1[i] + l*sin(rad[i]);
    }
  
    background(0);
  
    for(int i = 0; i < width/d; i ++){ //線を描く
        stroke(255);
        line(0, d*i, width, d*i);
    }
  
    for(int i = 0; i < n; i ++){ //針を描く
        for(int j = 0; j < width/d; j ++){
            if(y1[i] <= d*j && y2[i] >= d*j){ //交わっているとき
                line(x1[i], y1[i], x2[i], y2[i]);
                nc ++; //交わっている針の数を+1する
            }                
            if(y1[i] > d*j && y2[i] < d*(j+1)){ //交わっていないとき
                line(x1[i], y1[i], x2[i], y2[i]);
            }
        }
    }
    
    float pi = n/nc; //piの計算
    println("crossed lines = " + round(nc)); //交わっている針の本数を表示
    println("calculated pi = " + pi); //piの計算値を表示
}