プログラミングの備忘録

プログラムをつくる過程を残すもの

processingの備忘録 -ランジュバン方程式-

こんにちは。
今回は「ランジュバン方程式」をprocessingで可視化(?)してみます。


本題とは関係無い話になりますが、最近はずっと競プロとProject Eulerをやっていたのでpythonばかり書いて、1ヶ月ほどでしょうか、processingに全く触れていませんでした。
そろそろやらないと忘れそうだということで、競プロ復習編は気になったものだけ個人的にやることにして、Project Eulerも50問目までやって切り良く終わりにしました。
(今後一切やらないという訳ではないのでまたたまに現れるかもしれません。)

ではリハビリがてら何をつくろうかと思っていたところに「ランジュバン方程式」を見つけました。


ランジュバン方程式とは

ランジュバン方程式はブラウン運動を数式で表したもので、以下のような形です。

$$ m \boldsymbol{a} = - \beta \boldsymbol{v} + \eta (t) $$

ここで、$\beta$は抵抗係数、$\eta (t)$は確率的にかかる力を表しています。

ある物体が、速度に比例した抵抗力$- \beta \boldsymbol{v}$と周りからランダムにかかる力$\eta (t)$を受けながら運動する、という感じです。
wikipediaを軽く読んで書けそうだと思って選んだだけなので間違った解釈になっているかもしれませんが…)


おそらくランダムにかかる力$\eta (t)$の項だけにした場合のものが「ランダムウォーク」に相当するものだと思います。
(今のところ、抵抗の部分$- \beta \boldsymbol{v}$はあってもなくても似たものだろうと考えています。)

プログラム化する

では早速書いてみます。
大したことはしておらず、ただ位置、速度、加速度の計算をしてある物体に適用する、というだけです。

PVector pos, vel, acc;
float v_max;
PVector r;
ArrayList<PVector> history = new ArrayList<PVector>();

void setup(){
    size(500, 500);
    noStroke();
    
    pos = new PVector(width/2, height/2);
    vel = PVector.random2D();
    v_max = 2;
    acc = new PVector();
}

void draw(){
    background(0);
    
    //軌跡
    history.add(new PVector(pos.x, pos.y));
    for(PVector v: history){
        fill(255, 0, 0);
        ellipse(v.x, v.y, 2, 2);
    }
    
    //物体
    fill(255);
    ellipse(pos.x, pos.y, 5, 5);
    
    //ランジュバン方程式
    acc = new PVector(-vel.x, -vel.y);
    acc.setMag(0.1);
    r = PVector.random2D();
    acc.add(r);
    
    //物体の移動
    pos.add(vel);
    vel.add(acc);
    vel.limit(v_max);
    acc.mult(0);
    
    //壁での反射
    if(pos.x < 0 || width < pos.x) { vel.x *= -1; }
    if(pos.y < 0 || height < pos.y){ vel.y *= -1; }
}


位置を保存しておく配列historyをつくり、for文でその要素全てに対してellipse()で軌跡を描画しています。


実行すると、白い丸が動いてブラウン運動の軌跡を描いていきます。

具体的にはやりませんが、ランジュバン方程式の左辺をいじるとブラウン運動の挙動も変わります。
しかもこの部分は「力」という感覚的にわかりやすいものなので、例えば上昇気流や風があって下や左右から力がかかるとか、重力があって上から力がかかるとか、というのを加えてみると違った動きになって面白いかと思います。

抵抗係数やランダムさをいじってみるだけでも結構変わります。


過去にもブラウン運動を描画するプログラムを書いていますが、こちらはこちらで違った方法で描いているので良ければ読んでみてください。
(ランダムな力$\eta (t)$の部分を実際にシミュレーションして求めているような形になっています。)

taq.hatenadiary.jp


まとめ

以上、あっさりでしたがprocessingのリハビリがてら「ランジュバン方程式」を実装してみました。