方程未完全求解

Question

我一直在尝试求解二维向量 P 的方程。但是solve之后rhs上还有一些P。这是否意味着 Maxima 无法做到或我做错了什么？

这是它：

load("vect");
declare(".", commutative);

declare(P, nonscalar);
declare([v1,V1,r1], nonscalar);
declare([v2,V2,r2], nonscalar);
declare([w1,W1,m1,I1,w2,W2,m2,I2], scalar);

/* Revolute Constraint */
constraint: v2 + (w2~r2) - (v1 + (w1~r1)) = 0$

/* Velocities after impulse P */
eq1: v1 = V1 - P/m1$
eq2: w1 = W1 - (r1~P) / I1$
eq3: v2 = V2 + P/m2$
eq4: w2 = W2 + (r2~P) / I2$

eq: subst([eq1,eq2,eq3,eq4], constraint)$
solve(eq, P);

（我试图得到一个满足约束条件的冲量方程。我正在关注 Dirk Gregorius 的第二个 post：https://gamedev.net/forums/topic/469531-revolute-joint-usingimpulses/4086845)

Answer 1

按照 Robert Dodier 的建议，我分解了所有向量并分别求解 P[1] 和 P[2]。我有一些东西可以给我答案，但现在我怎样才能把它变成漂亮的矢量形式呢？这是：

load("vect");
declare(".", commutative);

declare(P, nonscalar);
declare([v1,V1,r1], nonscalar);
declare([v2,V2,r2], nonscalar);
declare([w1,W1,m1,I1], scalar);
declare([w2,W2,m2,I2], scalar);

cross_scalar_vector(s,v) := [-s*v[2], s*v[1]]$
/* Revolute Constraint on Linear Velocity */
constraint: v2 + cross_scalar_vector(w2,r2) - (v1 + cross_scalar_vector(w1,r1)) = [0,0]$

/* Sub in velocities after impulse P. */
post_velocities: [
    v1 = V1 - P/m1,
    w1 = W1 - (r1~P) / I1,
    v2 = V2 + P/m2,
    w2 = W2 + (r2~P) / I2
]$
constraint: subst(post_velocities, constraint)$

/* Break up the remaining vectors for solve. */
vectors: [
    P = [P[1], P[2]],
    V1 = [V1[1], V1[2]],
    r1 = [r1[1], r1[2]],
    V2 = [V2[1], V2[2]],
    r2 = [r2[1], r2[2]]
]$
constraint: subst(vectors, constraint)$

/* Break up vector constraint into x and y constraint for solve. */
xconstraint: lhs(constraint)[1] = 0$
yconstraint: lhs(constraint)[2] = 0$

/* Not sure why we need to do this again? */
xconstraint: subst(vectors, xconstraint)$
yconstraint: subst(vectors, yconstraint)$

/* Expand cross products for solve. */
xconstraint: express(xconstraint)$
yconstraint: express(yconstraint)$

solve([xconstraint,yconstraint], [P[1],P[2]]);

Answer 2

我想我已经了解了细节。我不得不手工做一些事情，而 Maxima 主要只是检查我是否做对了。如果目标只是找到解决方案，我想那没关系。

这是我的脚本。您可以通过以下方式执行此操作：maxima --batch=foo.mac 其中 foo.mac 是脚本的名称。

/* adapted from: 
 */

load("vect");

/* I don't want to reorder arguments of cross product. */
remrule ("~", ?\~rule4);

/* I do want to flatten noncommutative multiplication.
 * (It appears that's disabled by vect; ordinarily it happens by default.)
 */
dotassoc: true $

declare(P, nonscalar);
declare([v1,w1,V1,W1,r1], nonscalar);
declare([v2,w2,V2,W2,r2], nonscalar);
declare([m1,I1,m2,I2], scalar);

/* Revolute Constraint */
constraint: v2 - (r2~w2) - (v1 - (r1~w1)) = 0$

/* Velocities after impulse P */
eq1: v1 = V1 - P/m1$
eq2: w1 = W1 - (r1~P) / I1$
eq3: v2 = V2 + P/m2$
eq4: w2 = W2 + (r2~P) / I2$

eq: subst([eq1,eq2,eq3,eq4], constraint);

A(a) := matrix([0, -a[3], a[2]], [a[3], 0, -a[1]], [-a[2], a[1], 0]);

eqa: ev (eq, (u~v) := 'A(u).v);

matchdeclare (pp, lambda ([e], not freeof(P, e)));
matchdeclare (qq, lambda ([e], freeof(P, e)));
defrule (rp, pp + qq = 0, pp = -qq);

eqa1: expand (eqa);
eqa2: apply1 (eqa1, rp);

matchdeclare (aa, lambda ([e], matrixp(e) or listp(e)));
tellsimpafter (I(aa), ident (length (aa)));

matchdeclare ([aa, bb], all);
tellsimpafter (I(aa) . bb, bb);
tellsimpafter (aa . I(bb), aa);

M: -(1/I2)*'A(r2).'A(r2) - (1/I1)*'A(r1).'A(r1) + (1/m2)*I(P) + (1/m1)*I(P);
N: 'A(r2).W2 - 'A(r1).W1 - V2 + V1;

eqa2_factored: M . P = N;

expand (eqa2_factored);
?resimplify (%);
if % # eqa2 then error ("tried to factor eqa2, but something went wrong.");

solution: P = M^^-1 . N;

/* EXAMPLE: */

I1: 20 $
I2: 3 $
m1: 100 $
m2: 12 $
V1: [17, 19, -23] $
V2: [-5, -3, 11] $
W1: [8, 4, 14] $
W2: [-6, -16, 24 ] $
r1: [1/2, 2/3, 3/4] $
r2: [5, 7, 3] $

/* note various subterfuges to ensure evaluation with stated values */

example_M: ev (subst (I(P) = ident(3), M), nouns);
example_N: ev (N, nouns);

example_P: example_M^^-1 . example_N;

subst (P = example_P, ev (eqa2, eval, nouns));

if lhs(%) = rhs(%)
    then print ("TEST PASSED: lhs = rhs")
    else error ("TEST FAILED: lhs # rhs");

如果您需要为不同的参数 r1、r2 等评估 P，我的建议是评估矩阵 M 和向量 N 用你想插入的任何值，然后求解方程 P = M^^-1 . N。明确的解决方案可能会非常混乱。

方程未完全求解

Equation not fully solved for

maxima