p5.jsで非円形摩擦車を作成する

はじめに

　この記事は去年の７月に制作した非円形摩擦車のコードの説明です。非円形摩擦車とは次のようなものです。

　このように片方の車が円形でなく、もう片方と接触しながら回転します。すべての時間において接触点での速度の大きさと方向が一致しています。これを計算で出そうと思います。

　便利なのでp5.jsで書きました。バージョンは1.11.2です。

コード全文

　strange form wheel remake

/*
  数値の算出元：https://editor.p5js.org/dark_fox/sketches/X4N0iyFCO
  ニュートン法を内側から外側に向かって使う
  ニュートンで書き直しました。お疲れさま。
*/

/*
  非円形摩擦車
  参考：
  https://www.youtube.com/watch?v=uqlpZN6eX2o
  https://www.youtube.com/watch?v=xbONXxrTIcc
*/

const master = {};
const WHEEL_NUM = 4;
const CANVAS_SCALE = Math.floor(Math.min(window.innerWidth, window.innerHeight)*0.99);
const CANVAS_WIDTH = CANVAS_SCALE;
const CANVAS_HEIGHT = CANVAS_SCALE;
const SCALE_FACTOR = CANVAS_SCALE*0.02;

// インテグラル0からtまでのfという関数の1000刻みでのt:0～1の値を格納するための関数
function getIntegralArray(f,dt=1000){
  const result=[0];
  let sum=0;
  for(let i=0;i<dt;i++){
    sum+=(1/dt)*f(i/dt);
    result.push(sum);
  }
  // xに対してresult[x/dt]が疑似的にfの積分になる
  // 間の場合は補間する
  // たとえばfloor(x/dt)やfract(x/dt)を取るなどする
  // 0～dtに入るように適切に制限する
  return result;
}

// newton法関連

// radius関数の0～1におけるdetail刻みでの最大値を出すだけ
// なぜ必要なのかは後述
function maxCurveRadius(radius,dt=1000){
  let r=-Infinity;
  for(let i=0;i<dt;i++){
    r=Math.max(radius(i/dt),r);
  }
  return r;
}


function getIntegralOuter(radius,diffTheta,a,dt=1000){
  let sum=0;
  for(let i=0; i<dt; i++){
    sum+=(1/dt)*diffTheta(i/dt)*radius(i/dt)/(a-radius(i/dt));
  }
  return sum;
}

function getIntegralPower(radius,diffTheta,a,dt=1000){
  let sum=0;
  for(let i=0; i<dt; i++){
    sum+=(1/dt)*diffTheta(i/dt)*radius(i/dt)/Math.pow(a-radius(i/dt),2);
  }
  return sum;
}

// 精度やばいです
function applyNewtonOuter(radius, diffTheta, target = Math.PI*2, dt=1000){
  const r = maxCurveRadius(radius,dt);
  // r+0.1とかr+0.01
  let initialValueDiff = 0.1;
  while(1){
    if(getIntegralOuter(radius, diffTheta, r+initialValueDiff) > target) break;
    initialValueDiff *= 0.1;
  }
  let a0 = r + initialValueDiff;
  for(let k=0; k<20; k++){
    a0 += (getIntegralOuter(radius,diffTheta,a0,dt)-target)/getIntegralPower(radius,diffTheta,a0,dt);
    console.log(a0);
  }
  return a0;
}

function setup() {
  createCanvas(CANVAS_WIDTH, CANVAS_HEIGHT);

  const r1 = (t) => 6+sin(TAU*t)+cos(TAU*t)+sin(2*TAU*t)+sin(4*TAU*t);
  // ここをいじることで不等速回転にできる？
  const theta1 = (t) => TAU*t;
  const diffTheta1 = (t) => TAU;
  master.r1 = r1;
  master.theta1 = theta1;

  const detail = 800;

  //const a = 12.717610663426072;
  const a = applyNewtonOuter(r1, diffTheta1, TAU, detail);
  
  const f = (t) => diffTheta1(t)*r1(t)/(a-r1(t));
  const intF = getIntegralArray(f,detail);
  //console.log(TAU-intF[800]); // ほぼTAU

  const r2 = (t) => a-r1(t);
  // 要するにtheta2はintFなんですが値が離散的なので（配列で与えられるので）
  // 前後の値を算出したうえで補間で出しています
  const theta2 = (t) => {
    const x = fract(t);
    const nx = floor(x*detail);
    const fx = fract(x*detail);
    const leftValue = intF[nx];
    const rightValue = intF[nx+1];
    return leftValue*(1-fx) + rightValue*fx;
  }
  master.a = a;
  master.r2 = r2;
  master.theta2 = theta2;

  const pathDetail = 400;

  const path1 = new Path2D();
  for(let h=0;h<=pathDetail;h++){
    const t = h/pathDetail;
    const x = SCALE_FACTOR*r1(t)*cos(theta1(t));
    const y = SCALE_FACTOR*r1(t)*sin(theta1(t));
    if(h===0){
      path1.moveTo(x,y);
    }else{
      path1.lineTo(x,y);
    }
  }
  master.path1 = path1;

  const path2 = new Path2D();
  for(let h=0;h<=pathDetail;h++){
    const t = h/pathDetail;
    const x = SCALE_FACTOR*r2(t)*cos(theta2(t));
    const y = SCALE_FACTOR*r2(t)*sin(theta2(t));
    if(h===0){
      path2.moveTo(x,y);
    }else{
      path2.lineTo(x,y);
    }
  }
  master.path2 = path2;
}

function draw(){
  background(0);
  const {a,r1,r2,theta1,theta2,path1,path2} = master;
  
  noStroke();
  const elapsedTime = millis()/3000;

  fill(64,128,255);
  push();
  translate(CANVAS_WIDTH*0.5, CANVAS_HEIGHT*0.5);
  // 全体を反時計回りに回しているのでtimeのあとの右側の位置は-timeの時の位置が来る。
  // マイナスになるのはそういうこと。
  // 座標軸をマイナスに回してそれに基づいてpathで描画しているところをイメージすればわかる。
  rotate(-theta1(elapsedTime));
  drawingContext.fill(path1);
  pop();

  fill(64,128,255,128);
  for(let k=0;k<WHEEL_NUM;k++){
    const phi = k*TAU/WHEEL_NUM;
    push();
    // 中心位置に移る計算
    translate(CANVAS_WIDTH*0.5 + a * SCALE_FACTOR*cos(phi), CANVAS_HEIGHT*0.5 + a*SCALE_FACTOR*sin(phi));
    // 鏡映変換で、軸に垂直な直線による折り返しを実行する。
    applyMatrix(-cos(2*phi),-sin(2*phi),-sin(2*phi),cos(2*phi),0,0);
    // 詳細は後述：
    rotate(-theta2(elapsedTime+k/WHEEL_NUM) - theta1(elapsedTime) + theta1(elapsedTime+k/WHEEL_NUM));
    drawingContext.fill(path2);
    pop();
  }
}

実行結果：

中央の車の設計

　極座標で次のように書きます。

  const r1 = (t) => 6+sin(TAU*t)+cos(TAU*t)+sin(2*TAU*t)+sin(4*TAU*t);
  // ここをいじることで不等速回転にできる？
  const theta1 = (t) => TAU*t;
  const diffTheta1 = (t) => TAU;

常に$r_1>0$です。theta1は角度の関数で、$0$から$1$まで変数が動き、$2\pi$まで変化します。これとかみあう式を求めることになります。活躍するのは、微分と積分です。

数学的要請

　二つの摩擦車の式を$r_1,r_2$とするとき、これらは次の関係を満たす必要があります。

r_1(t)+r_2(t)=a,~~~r_1(t)\dot{\theta_1}(t)=r_2(t)\dot{\theta_2}(t).

ただし$a$はある定数とします。このとき、次の曲線：

x_1(t)=r_1(t)\cos\theta_1(t),~~~~y_1(t)=r_1(t)\sin\theta_1(t)~~~(0\leq t\leq 1)

と

x_2(t)=a-r_2(t)\cos\theta_2(t),~~~~y_2(t)=r_2(t)\sin\theta_2(t)~~~(0\leq t \leq 1)

を考えます。一つ目の曲線を原点中心で負の向きに回して、二つ目の曲線を$(a,0)$中心で正の向きに回します。回す量はそれぞれ

\phi_1 = \theta_1(t_0),~~~\phi_2 = \theta_2(t_0)

です。回したうえで$t=t_0$における位置と速度を見ます。回した結果は、

x_1(t)=r_1(t)\cos(\theta_1(t)-\phi_1),~~~~y_1(t)=r_1(t)\sin(\theta_1(t)-\phi_1)

と

x_2(t)=a-r_2(t)\cos(\theta_2(t)-\phi_2),~~~~y_2(t)=r_2(t)\sin(\theta_2(t)-\phi_2)

です。なので$t=t_0$のときこれらは一致します。微分の結果についてですが、

\begin{align}
\dot{x_1}(t) &= \dot{r_1}(t)\cos(\theta_1(t)-\phi_1)-r_1(t)\dot{\theta_1}(t)\sin(\theta_1(t)-\phi_1),\\
\dot{y_1}(t) &= \dot{r_1}(t)\sin(\theta_1(t)-\phi_1)+r_1(t)\dot{\theta_1}(t)\cos(\theta_1(t)-\phi_1),\\
\dot{x_2}(t) &= -\dot{r_2}(t)\cos(\theta_2(t)-\phi_2)+r_2(t)\dot{\theta_2}(t)\sin(\theta_2(t)-\phi_2),\\
\dot{y_2}(t) &= \dot{r_2}(t)\sin(\theta_2(t)-\phi_2)+r_2(t)\dot{\theta_2}(t)\cos(\theta_2(t)-\phi_2)
\end{align}

なので、

r_1(t)+r_2(t)=a,~~~~~\dot{r_1}(t)=-\dot{r_2}(t)

であることより、微分の方の条件と合わせて、

\dot{x_1}(t_0)=\dot{x_2}(t_0),~~~~~~\dot{y_1}(t_0)=\dot{y_2}(t_0)

が成立します。これらのことより、あの条件を満たす$r_2,\theta_2$を求めて、これに従って描画すればいいことになります。すなわち負の向きと正の向きと互いに逆方向に回すわけです。

関数を求める

　まず

r_2(t)=a-r_1(t)

です。ここで$a$はまだ未知数です。わかるのは定数ということだけです。次に角度が、

\dot{\theta_2}(t)=\frac{\dot{\theta_1}(t)r_1(t)}{a-r_1(t)}

で決まるので、これの積分として、

\theta_2(t)=\int_0^t{\frac{\dot{\theta_1}(t)r_1(t)}{a-r_1(t)}}~dt

とすればいいんですが、ここでかみ合うために、

\theta_2(1)=2\pi

を課す必要があります。そうなると被積分関数が周期$1$なので、

\theta_2(t+1)=2\pi+\theta_2(t)

が成立して$\theta_2$の方の曲線も周期的になります。なので、$a$はこれを満たす必要があります。

\int_0^{1}{\frac{\dot{\theta_1}(t)r_1(t)}{a-r_1(t)}}~dt = 2\pi.

ニュートン法

　これを求めるためにニュートン法を使います。

I(a)=\int_0^{1}{\frac{\dot{\theta_1}(t)r_1(t)}{a-r_1(t)}}~dt - 2\pi

とおくと、これは分子が正なので、$a$が大きくなるほど$-2\pi$に近づいていく下に凸の関数です。その曲線を大きい方から右に進んでいく場合、関数$r_1$の$0$～$1$までのうち一番大きい値の近くで正の無限大に発散することが分かるでしょうか。

というのも$a-r_1(t)$が正の凄く小さい値になるからですね。ここら辺に発散のポイントがあって、そこから下に凸で一定値に近づいていくので、唯一のゼロ点が存在します。これを求めればいいことが分かりました。そこでニュートン法です。

　ニュートン法とは図のように、接線を求めてそれと軸との交点を割り出すことを繰り返してゼロ点に近づいていく手法です。状況に応じた様々な条件が必要ですが、今の場合はゼロ点より小さいことが確実な値を取ってそこから正の方向に近づいていく手法でいけます。それを取っているのが

  const r = maxCurveRadius(radius,dt);
  // r+0.1とかr+0.01
  let initialValueDiff = 0.1;
  while(1){
    if(getIntegralOuter(radius, diffTheta, r+initialValueDiff) > target) break;
    initialValueDiff *= 0.1;
  }
  let a0 = r + initialValueDiff;

ここですね。getIntegralOuterで出しているのが$I(a)$だ思ってください。最大半径を取るんですが、その値になんか足してゼロ点を上回ってしまっては困るので、小さくなるように修正を施しています。ただ発散ポイントに達しても困るので乗算で減らしています。
　それ以降の処理でニュートン法を使っています... が、$I'(a)$の算出で微分と積分の順序を交換しています。厳密にはいろいろな条件が必要なこの処理ですが、この文脈では安心安全の可積分なので、問題なく実行できます。

パスの作成

　関数ができたらパスを作ります。その前に$\theta_2$ですが、積分の形です。これはどうしているかというと、

  const theta2 = (t) => {
    const x = fract(t);
    const nx = floor(x*detail);
    const fx = fract(x*detail);
    const leftValue = intF[nx];
    const rightValue = intF[nx+1];
    return leftValue*(1-fx) + rightValue*fx;
  }

このように離散値の左右を取って割合で線形補間しています。パスを作ってるのはここ：

  const path1 = new Path2D();
  for(let h=0;h<=pathDetail;h++){
    const t = h/pathDetail;
    const x = SCALE_FACTOR*r1(t)*cos(theta1(t));
    const y = SCALE_FACTOR*r1(t)*sin(theta1(t));
    if(h===0){
      path1.moveTo(x,y);
    }else{
      path1.lineTo(x,y);
    }
  }

特に工夫もなく、moveToとlineToを組み合わせているだけです。path2も同様にして作ります。再利用した方が高速で描画できるのでそうしています。
　なお、masterという形でグローバル変数を一ヶ所にまとめることで直接アクセスしなくて済むので、変数名の節約ができています。

描画

　いよいよ描画ですが、まず中央の摩擦車は先ほど述べた通り逆回転で描画します。

  fill(64,128,255);
  push();
  translate(CANVAS_WIDTH*0.5, CANVAS_HEIGHT*0.5);
  // 全体を反時計回りに回しているのでtimeのあとの右側の位置は-timeの時の位置が来る。
  // マイナスになるのはそういうこと。
  // 座標軸をマイナスに回してそれに基づいてpathで描画しているところをイメージすればわかる。
  rotate(-theta1(elapsedTime));
  drawingContext.fill(path1);
  pop();

その周囲については、回す方向が逆なので、鏡写しにしたうえで同じ方向に回しています。鏡写しはapplyMatrixを使っています。計算式の算出は省略します。4を3にしてみてください。ちゃんとかみ合います。ただ5とかで重なるのは防げないのでご了承ください。これは$r_1$からして4のときにうまく重ならないようなものを用意しているだけです。

  fill(64,128,255,128);
  for(let k=0;k<WHEEL_NUM;k++){
    const phi = k*TAU/WHEEL_NUM;
    push();
    // 中心位置に移る計算
    translate(CANVAS_WIDTH*0.5 + a * SCALE_FACTOR*cos(phi), CANVAS_HEIGHT*0.5 + a*SCALE_FACTOR*sin(phi));
    // 鏡映変換で、軸に垂直な直線による折り返しを実行する。
    applyMatrix(-cos(2*phi),-sin(2*phi),-sin(2*phi),cos(2*phi),0,0);
    // 詳細は後述：
    rotate(-theta2(elapsedTime+k/WHEEL_NUM) - theta1(elapsedTime) + theta1(elapsedTime+k/WHEEL_NUM));
    drawingContext.fill(path2);
    pop();
  }

角度についてはめんどくさいんですが、今回theta1が等速なので、これでOKです。

おわりに

　工夫すると歯車も作れるそうですが、自分はやり方を知りません...
　ちなみにここまで微分と積分を縦横に使うコードもあんま書きませんね。とかく無味乾燥と言われがちなやつらなので、これからも仲良くしてあげてください。
　ここまでお読みいただいてありがとうございました。

予告：内側もやるよ！

　内側についてもやります。