matarillo.com

The best days are ahead of us.

グラフの自動レイアウトに挑戦 #2 ばねモデル

2022-08-12 20:47:22

作成したアプリケーションとソースコード

物理シミュレーション: 運動方程式とオイラー法

グラフ構造の次は、物理シミュレーションの準備をしよう。シミュレーションでは、座標 (x, y)、速度 (vx, vy)、加速度 (ax, ay)、力 (fx, fy)という、2次元の値を使用する。

JavaScriptには演算子オーバーロードがないので、単純に {x: number, y: number} というオブジェクトで扱う。 前回のPart1では、ノードに座標Rと速度Vのプロパティを定義しておいた。

ところで、このノードの動きをシミュレートするにはどうすればいいのだろうか。そう、物理(力学)の法則、ニュートンの運動方程式だ。中学の理科では「力=重さ×加速度」と習ったと思う。

F = m·a

ところで、加速度 a というのは速度 v の微分(≒時間変化量)だし、速度 v というのは位置 r の微分(≒時間変化量)なので、運動方程式を位置 r を使って書き直すと、位置 r を2回微分したものが式の中に出てくる。

r'' = F/m

物体の位置を正しく求めるためには、力 F を2回積分しないといけない。これは非常に面倒くさいので、積分をしないで近似値を求めることにする。近似値の求め方もいろいろあるのだが、今回はオイラー法と呼ばれる方法を使う。

オイラー法は、時間の流れを短い間隔 Δt で区切ったときの、ある瞬間の値とその次の瞬間の値の差を単純な近似値で求める方法だ。速度の差分は 加速度 × Δt、位置の差分は 速度 × Δt で近似する。精度の高い近似とはいえないのだが、今回は精度の高いシミュレーションをする必要がないので、単純なオイラー法を選択した。

a_k = f_k/m, v_k+1 = v_k + (a_k·Δt), r_k+1 = r_k + (v_k·Δt)

オイラー法をコードに落とすとこんなイメージになる。ちなみに、面倒なのでノードの質量を全部1だと決めた。欲しいのは厳密なシミュレーションではないので、縮尺とか比率に関係するところは適当に決めてしまってよい。

const move = (node, dt, f) => {
    const a = f; // 運動方程式: 質量は1とする
    const r = {
        x: node.r.x + dt * node.v.x,
        y: node.r.y + dt * node.v.y
    };
    const v = {
        x: node.v.x + dt * a.x,
        y: node.v.y + dt * a.y
    };
    node.r = r;
    node.v = v;
};

引数は3つある。最初の引数はノードで、2つ目の引数 dt は定数だ。この値が小さいほど、オイラー法の誤差が小さくなる。今回は dt = 0.1 とする。

3つめの引数 f が、このノードにかかる力だ。では、この力を計算しよう。

ノードにかかる力: ばねの弾性力

ノードとノードをつなぐエッジがばねでできているので、あるノードがばねから受ける力を考える。

ばねの弾性力の図

こんな感じで、ノード間の距離を d 、ばねの自然長をlとしたとき、右上のノードがばねから受ける力 F

F = -k·(d-l)

で表される。ちなみに、定数 k はばねの強さを表す比例定数だ。

コードを書くときは、kl の値は好きに決めていい。後から修正して変えてしまってもかまわない。

const getSpringForce = (n0, n1) => {
    const dx = n0.r.x - n1.r.x;
    const dy = n0.r.y - n1.r.y;
    const dl = Math.sqrt(dx * dx + dy * dy);
    if (dl < Number.EPSILON) {
        const randomForces = [
            { x: 0, y: -1 },
            { x: 1, y: 0 },
            { x: 0, y: 1 },
            { x: -1, y: 0 }
        ];
        return randomForces[Math.floor(Math.random() * randomForces.length)];
    }
    const cos = dx / dl;
    const sin = dy / dl;
    const k = 0.1; // ばね定数
    const l = 60; // ばねの自然長
    const f = -k * (dl - l); // ばねの弾性力
    return { x: f * cos, y: f * sin };
};

途中に不思議なコードが入っているのは、ノード間の距離が0のときの対策だ。現実世界の物体だとありえないが、このシミュレーションでは起こりえる。そうすると角度 θ が計算できなくなってしまう。なので、乱数を使ってランダムな方向に力を向けているのだ。ちょっと強引だが気にしない。

ノードにかかる力: ノードの反発力

ノード同士の反発力は、ばねでつながっているかどうかに関係なく、すべてのノード間で働く。

ノードの反発力の図

反発力は、電磁力や重力を参考に、逆2乗の法則に従うことにする。つまり、力 F の大きさが距離 d の2乗に反比例するものとする。

F = g/d^2

ちなみに、定数 g は反発力の強さを表す比例定数だが、ばねのときと同じように、好きな値に決めていい。

const getReplusiveForce = (n0, n1) => {
    const dx = n0.r.x - n1.r.x;
    const dy = n0.r.y - n1.r.y;
    const dldl = dx * dx + dy * dy;
    const dl = Math.sqrt(dldl);
    if (dl < Number.EPSILON) {
        const randomForces = [
            { x: 0, y: -1 },
            { x: 1, y: 0 },
            { x: 0, y: 1 },
            { x: -1, y: 0 }
        ];
        return randomForces[Math.floor(Math.random() * randomForces.length)];
    }
    const cos = dx / dl;
    const sin = dy / dl;
    const g = 500; // 比例定数
    const f = g / dldl; // 2次元での反発は距離に反比例
    return { x: f * cos, y: f * sin };
};

さっきと同じで、距離0のときはランダムに力を向けている。今回は、距離0のときに力の向きが決まらないだけでなく、力の大きさが無限大になってしまう。

ノードにかかる力: 摩擦力

さて、Eadesのばねモデルはこれで完璧のような気がするが、ここで終わりにすると想定外の結果になる。いつまでたってもノードが動き続けていて止まらないのだ。やりたかったことは自動レイアウトなのだから、いつまでも動き続けているのは具合が悪い。

ではどうすればいいかというと、摩擦力を加えるとよい。要するに、速度がだんだん小さくなっていくような仕掛けが必要なのだ。

今回のモデルでは空気抵抗を参考にしよう。力が速度と逆向きに加わることにする。

空気抵抗(空気の摩擦力)の図

力Fの大きさは速度 v の大きさに比例する。

F = -μ·v

摩擦力でも、比例定数 μ は好きな値にしてよい。なお、ソースコード上は μ という記号の変わりにアルファベット m を用いた。

const getFrictionalForce = (n) => {
    const m = 0.3; // 比例定数 < 1
    return { x: -m * n.v.x, y: -m * n.v.y }; // 摩擦力は速度の逆向きに比例
};

第2回のまとめ

では、ここまでをまとめよう。

まず、すべてのノードについて、各ノードにかかる力を足し合わせてオイラー法を適用するのだ。そのような関数(厳密にはラムダ式)は moveAll という名前にする。

const moveAll = () => {
    const dt = 0.1; // 時間の差分
    for (let i = 0; i < nodes.length; i++) {
        const n0 = nodes[i];
        const sf = { x: 0, y: 0 };
        for (let j = 0; j < n0.links.length; j++) {
            const l = n0.links[j];
            const n1 = nodes[l];
            const f = getSpringForce(n0, n1);
            sf.x += f.x;
            sf.y += f.y;
        }
        const rf = { x: 0, y: 0 };
        for (let j = 0; j < nodes.length; j++) {
            if (i === j) continue;
            const n1 = nodes[j];
            const f = getReplusiveForce(n0, n1);
            rf.x += f.x;
            rf.y += f.y;
        }
        const ff = getFrictionalForce(n0);
        const f = {
            x: sf.x + rf.x + ff.x,
            y: sf.y + rf.y + ff.y
        };
        move(n0, dt, f);
    }
};

そして、更新されたノードの位置をSVGに反映する。

const updateElements = () => {
    for (const edge of edges) {
        const n0 = nodes[edge.n0];
        const n1 = nodes[edge.n1];
        const line = edge.element;
        line.setAttribute("x1", n0.r.x);
        line.setAttribute("y1", n0.r.y);
        line.setAttribute("x2", n1.r.x);
        line.setAttribute("y2", n1.r.y);
    }
    for (const node of nodes) {
        const g = node.element;
        g.style.transform = `translate(${node.r.x}px, ${node.r.y}px)`;
    }
};

最後に、Window.requestAnimationFrame() 関数を使ってアニメーションを描画する。

window.onload = () => {
    init();
    arrangeInCircle(200);
    createElements();

    // アニメーションを更新するコールバック関数
    function loop(ts) {
        moveAll();
        updateElements();
        window.requestAnimationFrame(loop);
    }
    window.requestAnimationFrame(loop);
};

Part2はここまで。実行すると、Eadesのばねモデルに従ってノードがアニメーションする。


グラフの自動レイアウトに挑戦 インデックスへ戻る