自己试着做个2d物理引擎玩时遇到的一个问题... 两个物体碰撞时计算摩擦冲量似乎很容易受浮点数精度的影响

两个物体碰撞时给予对方的碰撞冲量我按下面这个公式算了出来

其中速度差向量 vr = vA' - vB' , n是两物体分离轴的法线单位向量
另外把角速度也换算成线速度加到速度向量上来计算, 于是vA' = vA + rA×ωA
到这里都没什么问题

现在要求摩擦冲量, 令单位切线向量

$$tangent = normalize(vr - (vr*normal) * normal)$$

于是变成下图这样

按图1的上面的公式求得碰撞冲量大小j之后, 摩擦冲量大小$$jf = μ * j$$, 摩擦力在两个物体相对静止时就会变成0所以还要求一个最大摩擦冲量大小
vr在tangent上的投影的大小 $$vt = vr * tangent$$, 对两个物体施加摩擦冲量之后最多使得这个vt变成0但不能让vt的符号改变

我按照图1下面的4个公式以及△v = vt求得的最大jf是:
$$jfMax = vr * tangent / (1/mA + 1/mB + (rA×tangent)^2/IA + (rB×tangent)^2/IB)$$
于是摩擦冲量 $$tImpulse = tangent * min(jfMax, μ * j)$$

分别对物体A施加 -tImpulse, 对物体B施加 tImpulse, 如此实现摩擦

但是我发现在去除μ*j只适用jfMax的情况下, 最后vt并不等于0而是在0附近上下波动, 多数是0.01~0.5

 bodyA.applyImpulse(tangentImpulse.negate(), cp);
bodyB.applyImpulse(tangentImpulse, cp);
vA = bodyA.velocity.add(rA.scalarCross(bodyA.angularVelocity));
vB = bodyB.velocity.add(rB.scalarCross(bodyB.angularVelocity));
console.log(vA.sub(vB).dot(tangent));

适用摩擦冲量前后, 有时动能反而会增大

 k1 = bodyA.getKinetic() + bodyB.getKinetic();
bodyA.applyImpulse(tangentImpulse.negate(), cp);
bodyB.applyImpulse(tangentImpulse, cp);
k2 = bodyA.getKinetic() + bodyB.getKinetic();
if ((k2 - k1) > 0.1) console.log(((k2 - k1) / k1 * 100 | 0) + '%');

单个小物体会因为这个问题在与其它物体碰撞时明显的"抽搐"

求摩擦冲量的式子我翻了下别人的物理引擎基本都一样, 所以我觉得可能是浮点数的精度误差累积造成的问题
别人的代码里有看到一些误差修正的式子, 但搞不懂是怎么得出来的也不知道实际修正的是个啥...

另外我的一大坨代码如下

 applyImpulse: function () {
            var bodyA = this.bodyA,
                bodyB = this.bodyB,
                cp = this.contactPoint,
                restitution = bodyA.restitution > bodyB.restitution ? bodyA.restitution : bodyB.restitution,
                normal = this.normal,
                rA = cp.sub(bodyA.centroid),
                rB = cp.sub(bodyB.centroid),
                vA = bodyA.velocity.add(rA.scalarCross(bodyA.angularVelocity)),
                vB = bodyB.velocity.add(rB.scalarCross(bodyB.angularVelocity)),
                invIA = bodyA.inverseInertia,
                invIB = bodyB.inverseInertia,
                invMA = bodyA.inverseMass,
                invMB = bodyB.inverseMass,
                vr = vA.sub(vB),
                rsnA = rA.cross(normal),
                rsnB = rB.cross(normal),
                kn = invMA + invMB + rsnA * rsnA * invIA + rsnB * rsnB * invIB,
                j = -(1 + restitution) * vr.dot(normal) / kn,
                impulse = normal.mul(j);

            var tangent = vr.sub(normal.mul(vr.dot(normal))).normalize(),
                rstA = rA.cross(tangent),
                rstB = rB.cross(tangent),
                kt = invMA + invMB + rstA * rstA * invIA + rstB * rstB * invIB,
                jfMax = vr.dot(tangent) / kt,
                sf = Math.sqrt(bodyA.staticFriction * bodyB.staticFriction),
                df = Math.sqrt(bodyA.dynamicFriction * bodyB.dynamicFriction),
                jf = jfMax < Math.abs(j * sf) ? jfMax : Math.abs(j * df),
                tangentImpulse = tangent.mul(jf);

            bodyA.applyImpulse(impulse.add(tangentImpulse.negate()), cp);
            bodyB.applyImpulse(impulse.negate().add(tangentImpulse), cp);

            return this;
        }