D3.js 双摆怪行为

Question

我正在使用 D3.js 制作双摆模拟器。我正在使用 4 阶 Runge-Kutta 来积分 ODE，并且我很确定我在使用 mathamatica 检查输出数据时得到了正确的 theta 和 phi 值。然而，第二个钟摆的长度在移动时不断变化，这是不应该发生的，而且钟摆描绘出的路径似乎也不正确。由于我几乎可以肯定我从集成中获得了正确的 phi 和 theta 值，我认为问题可能出在我计算 x 和 y 坐标的方式上，但我不确定。

var width = 710, height = 710;
var margin = {top: -20, right: 30, bottom: 40, left: 40};
var g = 9.81; //so length is in meters 
var l = 0.4; //length of pendulum 1 (top pendulum)
var ell = 0.5; //ibid pendulum 2 (bottom pendulum)
var m = 5; //mass of pendulum 1
var M = 5; // ibid pendulum 2
var theta = 40 * (Math.PI/180); //angle of top pendulum
var phi = 10 * (Math.PI/180); //angle of bottom pendulum
var thetaDot = 0;
var phiDot = 0;

var xBall1 = l*Math.sin(theta);
var yBall1 = -l*Math.cos(theta);
var xBall2 = ell*Math.sin(phi);
var yBall2 = -ell*Math.cos(phi) + yBall1;
var t = 0;
var stepSize = .01;

function uDot(theta, thetaDot, phi, phiDot) { //first ODE
    return(thetaDot);
}

function vDot(theta, thetaDot, phi, phiDot) { //second ODE
    val = -g*(2*m + M)*Math.sin(theta) - M*g*Math.sin(theta - (2*phi)) - 2*Math.sin(theta-phi)*M*((phiDot*phiDot*ell) + thetaDot*thetaDot*l*Math.cos(theta-phi));

    val = val/( l*(2*m + M - M*Math.cos(2*theta - 2*phi)) );
    return(val);
}

function wDot(theta, thetaDot, phi, phiDot) { //third ODE
    return(phiDot);
}

function sDot(theta, thetaDot, phi, phiDot) {
    val = 2*Math.sin(theta-phi)*(thetaDot*thetaDot*l*(m+M) + g*(m+M)*Math.cos(theta) + phiDot*phiDot*ell*M*Math.cos(theta-phi));    
    
    val = val/( ell*(2*m + M - M*Math.cos(2*theta - 2*phi)) );
    return(val);
}

function RK4() { /* 4th order Runge-Kutta solver for system of 4 equations with 4 variables */

        k0 = stepSize * uDot(theta, thetaDot, phi, phiDot);
        l0 = stepSize * vDot(theta, thetaDot, phi, phiDot);
        q0 = stepSize * wDot(theta, thetaDot, phi, phiDot);
        p0 = stepSize * sDot(theta, thetaDot, phi, phiDot);       

        k1 = stepSize * uDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));
        l1 = stepSize * vDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));
        q1 = stepSize * wDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));
        p1 = stepSize * sDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));        

        k2 = stepSize * uDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));
        l2 = stepSize * vDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));
        q2 = stepSize * wDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));
        p2 = stepSize * sDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));        

        k3 = stepSize * uDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);
        l3 = stepSize * vDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);
        q3 = stepSize * wDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);
        p3 = stepSize * sDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);        

        theta = theta + ((1/6) * (k0 + 2*k1 + 2*k2 + k3));
        thetaDot = thetaDot + ((1/6) * (l0 + 2*l1 + 2*l2 + l3));
        phi = phi + ((1/6) * (q0 + 2*q1 + 2*q2 + q3));
        phiDot = phiDot + ((1/6) * (p0 + 2*p1 + 2*p2 + p3));

    }

d3.select("canvas").attr("width", width).attr("height", height);

var canvas = d3.select("canvas");

var context = canvas.node().getContext("2d");

var svg = d3.select("#doublePendulum")
    .attr("width", width)
    .attr("height", height)
    .append("g")
    .attr("transform",
        "translate(" + margin.left + "," + margin.top + ")");

var xAxis = d3.scaleLinear()
    .domain([-1.5, 1.5])
    .range([0, width]);

svg.append("g")
    .attr("transform", "translate(0," + height + ")")
    .call(d3.axisBottom(xAxis));

var yAxis = d3.scaleLinear()
    .domain([-1.5, 1.5])
    .range([height, 0]);

svg.append("g")
    .call(d3.axisLeft(yAxis));


var rod1 = svg.append('line')  /* put before circle code so it doesn't go over the circle */
    .attr('id', 'rod1')
    .attr('x2', xAxis(0))
    .attr('y2', yAxis(0))
    .attr('x1', xAxis(xBall1))
    .attr('y1', yAxis(yBall1))
    .attr('stroke', '#000000')
    .attr('stroke-width', '2px');

var ball1 = svg.append('circle')
    .attr('id', 'ball1')
    .attr('cx', xAxis(xBall1))
    .attr('cy', yAxis(yBall1))
    .attr('r', 10)
    .style('fill', '#000000');

var rod2 = svg.append('line') 
    .attr('id', 'rod2')
    .attr('x2', xAxis(xBall1))
    .attr('y2', yAxis(yBall1))
    .attr('x1', xAxis(xBall2))
    .attr('y1', yAxis(yBall2))
    .attr('stroke', '#000000')
    .attr('stroke-width', '2px');

var ball2 = svg.append('circle')
    .attr('id', 'ball2')
    .attr('cx', xAxis(xBall2))
    .attr('cy', yAxis(yBall2))
    .attr('r', 10)
    .style('fill', '#000000');    

function update() {
    context.beginPath();
    context.strokeStyle = '#0026FF';
    context.moveTo(xAxis(xBall2) + margin.left, yAxis(yBall2) + margin.top);    

    t += .01;

    RK4();
    xBall1 = l*Math.sin(theta);
    yBall1 = -l*Math.cos(theta);
    xBall2 = ell*Math.sin(phi);
    yBall2 = -ell*Math.cos(phi) + yBall1;

    context.lineTo(xAxis(xBall2) + margin.left, yAxis(yBall2) + margin.top);
    context.stroke();    


    d3.select('#rod1')
        .attr('y2', yAxis(0))
        .attr('x2', xAxis(0))
        .attr('y1', yAxis(yBall1))
        .attr('x1', xAxis(xBall1));

    d3.select('#ball1')
        .attr('cx', xAxis(xBall1))
        .attr('cy', yAxis(yBall1));

    d3.select('#rod2')
        .attr('y2', yAxis(yBall1))
        .attr('x2', xAxis(xBall1))
        .attr('y1', yAxis(yBall2))
        .attr('x1', xAxis(xBall2));

    d3.select('#ball2')
        .attr('cx', xAxis(xBall2))
        .attr('cy', yAxis(yBall2));

    }

var runApp = setInterval(function () { update(); }, 5);

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/5.7.0/d3.min.js"></script>
<div style="position: relative;">
    <canvas style="position: absolute;"></canvas>
    <svg id="doublePendulum" style="position: absolute; z-index: 1;"></svg>
</div>

Answer 1

你是对的，你没有正确计算 x,y 位置，更具体地说，你没有正确计算 x:

xBall1 = l*Math.sin(theta);
yBall1 = -l*Math.cos(theta);
xBall2 = ell*Math.sin(phi);
yBall2 = -ell*Math.cos(phi) + yBall1;

这与您计算初始位置的方式相同，但在初始化后，xBall2 依赖于 xBall1。

xBall1 = l*Math.sin(theta);
xBall2 = ell*Math.sin(phi) + xBall1;

在不查看其余内部结构的情况下，感觉不太对，但我不确定生成的模式应该是什么样子。但是，这样做可以使两个钟摆始终保持相同的长度，所以我敢说这是正确的：

var width = 710, height = 710;
var margin = {top: -20, right: 30, bottom: 40, left: 40};
var g = 9.81; //so length is in meters 
var l = 0.4; //length of pendulum 1 (top pendulum)
var ell = 0.5; //ibid pendulum 2 (bottom pendulum)
var m = 5; //mass of pendulum 1
var M = 1; // ibid pendulum 2
var theta = 40 * (Math.PI/180); //angle of top pendulum
var phi = 10 * (Math.PI/180); //angle of bottom pendulum
var thetaDot = 0;
var phiDot = 0;

var xBall1 = l*Math.sin(theta);
var yBall1 = -l*Math.cos(theta);
var xBall2 = ell*Math.sin(phi);
var yBall2 = -ell*Math.cos(phi) + yBall1;
var t = 0;
var stepSize = .01;

function uDot(theta, thetaDot, phi, phiDot) { //first ODE
    return(thetaDot);
}

function vDot(theta, thetaDot, phi, phiDot) { //second ODE
    val = -g*(2*m + M)*Math.sin(theta) - M*g*Math.sin(theta - (2*phi)) - 2*Math.sin(theta-phi)*M*((phiDot*phiDot*ell) + thetaDot*thetaDot*l*Math.cos(theta-phi));

    val = val/( l*(2*m + M - M*Math.cos(2*theta - 2*phi)) );
    return(val);
}

function wDot(theta, thetaDot, phi, phiDot) { //third ODE
    return(phiDot);
}

function sDot(theta, thetaDot, phi, phiDot) {
    val = 2*Math.sin(theta-phi)*(thetaDot*thetaDot*l*(m+M) + g*(m+M)*Math.cos(theta) + phiDot*phiDot*ell*M*Math.cos(theta-phi));    
    
    val = val/( ell*(2*m + M - M*Math.cos(2*theta - 2*phi)) );
    return(val);
}

function RK4() { /* 4th order Runge-Kutta solver for system of 4 equations with 4 variables */

        k0 = stepSize * uDot(theta, thetaDot, phi, phiDot);
        l0 = stepSize * vDot(theta, thetaDot, phi, phiDot);
        q0 = stepSize * wDot(theta, thetaDot, phi, phiDot);
        p0 = stepSize * sDot(theta, thetaDot, phi, phiDot);       

        k1 = stepSize * uDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));
        l1 = stepSize * vDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));
        q1 = stepSize * wDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));
        p1 = stepSize * sDot(theta+(0.5*k0), thetaDot+(0.5*l0), phi+(0.5*q0), phiDot+(0.5*p0));        

        k2 = stepSize * uDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));
        l2 = stepSize * vDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));
        q2 = stepSize * wDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));
        p2 = stepSize * sDot(theta+(0.5*k1), thetaDot+(0.5*l1), phi+(0.5*q1), phiDot+(0.5*p1));        

        k3 = stepSize * uDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);
        l3 = stepSize * vDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);
        q3 = stepSize * wDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);
        p3 = stepSize * sDot(theta+k2, thetaDot+l2, phi+q2, phiDot+p2);        

        theta = theta + ((1/6) * (k0 + 2*k1 + 2*k2 + k3));
        thetaDot = thetaDot + ((1/6) * (l0 + 2*l1 + 2*l2 + l3));
        phi = phi + ((1/6) * (q0 + 2*q1 + 2*q2 + q3));
        phiDot = phiDot + ((1/6) * (p0 + 2*p1 + 2*p2 + p3));

    }

d3.select("canvas").attr("width", width).attr("height", height);

var canvas = d3.select("canvas");

var context = canvas.node().getContext("2d");

var svg = d3.select("#doublePendulum")
    .attr("width", width)
    .attr("height", height)
    .append("g")
    .attr("transform",
        "translate(" + margin.left + "," + margin.top + ")");

var xAxis = d3.scaleLinear()
    .domain([-1.5, 1.5])
    .range([0, width]);

svg.append("g")
    .attr("transform", "translate(0," + height + ")")
    .call(d3.axisBottom(xAxis));

var yAxis = d3.scaleLinear()
    .domain([-1.5, 1.5])
    .range([height, 0]);

svg.append("g")
    .call(d3.axisLeft(yAxis));


var rod1 = svg.append('line')  /* put before circle code so it doesn't go over the circle */
    .attr('id', 'rod1')
    .attr('x2', xAxis(0))
    .attr('y2', yAxis(0))
    .attr('x1', xAxis(xBall1))
    .attr('y1', yAxis(yBall1))
    .attr('stroke', '#000000')
    .attr('stroke-width', '2px');

var ball1 = svg.append('circle')
    .attr('id', 'ball1')
    .attr('cx', xAxis(xBall1))
    .attr('cy', yAxis(yBall1))
    .attr('r', 10)
    .style('fill', '#000000');

var rod2 = svg.append('line') 
    .attr('id', 'rod2')
    .attr('x2', xAxis(xBall1))
    .attr('y2', yAxis(yBall1))
    .attr('x1', xAxis(xBall2))
    .attr('y1', yAxis(yBall2))
    .attr('stroke', '#000000')
    .attr('stroke-width', '2px');

var ball2 = svg.append('circle')
    .attr('id', 'ball2')
    .attr('cx', xAxis(xBall2))
    .attr('cy', yAxis(yBall2))
    .attr('r', 10)
    .style('fill', '#000000');    

// for debug:
var i = 0;
//


function update() {
    context.beginPath();
    context.strokeStyle = '#0026FF';
    context.moveTo(xAxis(xBall2) + margin.left, yAxis(yBall2) + margin.top);    

    t += .01;

    RK4();
    xBall1 = l*Math.sin(theta);
    yBall1 = -l*Math.cos(theta) ;
    xBall2 = ell*Math.sin(phi) + xBall1;
    yBall2 = -ell*Math.cos(phi) + yBall1;

    context.lineTo(xAxis(xBall2) + margin.left, yAxis(yBall2) + margin.top);
    context.stroke();    


    d3.select('#rod1')
        .attr('y2', yAxis(0))
        .attr('x2', xAxis(0))
        .attr('y1', yAxis(yBall1))
        .attr('x1', xAxis(xBall1));

    d3.select('#ball1')
        .attr('cx', xAxis(xBall1))
        .attr('cy', yAxis(yBall1));
        
   
    d3.select('#rod2')
        .attr('y2', yAxis(yBall1))
        .attr('x2', xAxis(xBall1))
        .attr('y1', yAxis(yBall2))
        .attr('x1', xAxis(xBall2));
        
    d3.select('#ball2')
        .attr('cx', xAxis(xBall2))
        .attr('cy', yAxis(yBall2));
        
        
    // for debug:     

   
    if(i%50 == 0) {
      var dx = yBall1 - yBall2;
      var dy = xBall1 - xBall2;
      var h1 = Math.sqrt(dx*dx+dy*dy);
      var h0 = Math.sqrt(yBall1*yBall1+xBall1*xBall1);
      
      console.log("Pendulum Lengths: " + h0 + ", " + h1 );
    }
    i++;
    
    //        

    }

var runApp = setInterval(function () { update(); }, 5);

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/5.7.0/d3.min.js"></script>
<div style="position: relative;">
    <canvas style="position: absolute;"></canvas>
    <svg id="doublePendulum" style="position: absolute; z-index: 1;"></svg>
</div>

D3.js 双摆怪行为

D3.js double pendulum odd behavior

javascript

d3.js

physics