在Python中求解符号线性方程组会得到意外的答案

Question

首先，请原谅即将到来的考试。我需要解下列由坐标变换产生的方程组：

给定v_n，v峎m，x峎p，y峎p，s，d，求解v_yp，v峎z p的系统：

I: v_n = v_yp * (1/(s*d)) + v_xp * y_p * (1/(s*d^2))  
II: v_m = v_zp * (1/(s*d)) + v_xp * z_p * (1/(s*d^2))  
III: v_xp = -1/Jp_11 * (Jp_12 * v_yp + Jp_13 * v_zp)  

术语Jp_ij是矩阵Jp的条目，可以写成：

^{pr2}$

b_0也给出了。现在我要解一，二，三的问题。我在SymPy中尝试了以下命令：

In [156]: EQ1 = v_n - (v_yp * (1/(s*d)) + yp * v_xp * (1/(s*d**2)))

In [157]: EQ2 = v_m - (v_zp * (1/(s*d)) + zp * v_xp * (1/(s*d**2)))

In [158]: EQ3 = v_xp + 1/Jp.row(0).col(0).tolist()[0][0] * (Jp.row(0).col(1).tolist()[0][0] * v_yp + Jp.row(0).col(2).tolist()[0][0] * v_zp)

In [159]: EQ1
Out[159]: v_n - v_yp/(d*s) - v_xp*yp/(d**2*s)

In [160]: EQ2
Out[160]: v_m - v_zp/(d*s) - v_xp*zp/(d**2*s)

In [161]: EQ3
Out[161]: v_xp + (v_yp*sin(phi)*sin(theta) + v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi))

In [162]: sympy.linsolve([EQ1,EQ2,EQ3], (v_xp, v_yp, v_zp))
Out[162]: {(0, d*s*v_n, d*s*v_m)}

In [163]: sympy.linsolve([EQ1,EQ2,EQ3], (v_yp, v_zp))
Out[163]: EmptySet()

In [164]: sympy.linsolve([EQ1,EQ2], (v_yp, v_zp))
Out[164]: {(d*s*v_n - v_xp*yp/d, d*s*v_m - v_xp*zp/d)}

In [165]: sympy.solve([EQ1,EQ2,EQ3], (v_yp, v_zp))
Out[165]: []

因为第三个方程已经设定了vuxp和vyu，vzp之间的关系，所以我只想用vyuyp和vuzp来求解系统。第163行没有给出我想要的答案，第162行的输出是出乎意料的：v_xp是指向太阳观测者的一个方向的速度，只有在观察磁盘中心时才应该为零（假设径向速度等于v_xp（仅在这一点上）被设为v_r=0）。在

我还试着把vxp（III）手动插入I和II。在

In [175]: v_xp = - (v_yp*sympy.sin(phi)*sympy.sin(theta) + v_zp*(-sympy.sin(b0)*sympy.sin(theta)*sympy.cos(phi) + sympy.cos(b0)*sympy.cos(theta)))/(sympy.sin(b0)*sympy.cos(theta) + sympy.sin(theta)*sympy.cos(b0)*sympy.cos(phi))

In [176]: v_xp
Out[176]: (-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi))

In [178]: EQ4 = v_n - v_yp/(s*d) - v_xp*yp/(s*d**2)

In [179]: EQ4
Out[179]: v_n - v_yp/(d*s) - yp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d**2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)))

In [181]: EQ5 = v_m - v_zp/(s*d) - v_xp*zp/(s*d**2)

In [182]: EQ5
Out[182]: v_m - v_zp/(d*s) - zp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d**2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)))

In [183]: sympy.linsolve([EQ4,EQ5], (v_yp, v_zp))
Out[183]: {(d*s*v_n, d*s*v_m)}

同样，我们得到的结果将对应于v\uxp=0。在

然而，当手工计算时，v_yp，v_zp的表达式会变成一些复杂的术语。如果需要的话，我可以稍后再发。我在做符号计算，因为我想检查结果，看看是否还能简化。在

求解这些方程并没有如预期的那样工作。为什么？在

---

PUSH：我在MATLAB中尝试了同样的任务：这是代码和结果。请注意，为了简单起见，我用Jij替换了Jp iu ij。在

代码：

syms xp yp zp v_xp v_yp v_zp s d v_n v_m b0 phi theta J11 J12 J13;

EQ1 = v_xp == (-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi));
EQ2 = v_n == v_yp/(d*s) + v_xp*yp/(s*d^2);
EQ3 = v_m == v_zp/(s*d) + v_xp*zp/(s*d^2);
[sol_vyp,sol_vzp] = solve([EQ1,EQ2,EQ3],[v_yp,v_zp]);
sol_vyp, sol_vzp
[sol_vyp] = solve([EQ2],[v_yp]);
sol_vyp

EQ4 = v_n == v_yp/(d*s) + yp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d^2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)));
EQ5 = v_m == v_zp/(d*s) + zp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d^2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)));
[sol_vyp] = solve([EQ4],[v_yp]);
sol_vyp
[sol_vyp,sol_vzp] = solve([EQ4,EQ5],[v_yp,v_zp]);
sol_vyp, sol_vzp

EQ6 = v_n == v_yp/(d*s) + yp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d^2);
EQ7 = v_m == v_zp/(d*s) + zp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d^2);
[sol_vyp] = solve([EQ6],[v_yp]);
sol_vyp
[sol_vyp,sol_vzp] = solve([EQ6,EQ7],[v_yp,v_zp]);
sol_vyp, sol_vzp

结果：

>> test_transform

sol_vyp =

Empty sym: 0-by-1

sol_vzp =

Empty sym: 0-by-1

sol_vyp =

d*s*(v_n - (v_xp*yp)/(d^2*s))

sol_vyp =

(v_n + (v_zp*yp*(cos(b0)*cos(theta) - cos(phi)*sin(b0)*sin(theta)))/(d^2*s*(sin(b0)*cos(theta) + cos(b0)*cos(phi)*sin(theta))))/(1/(d*s) - (yp*sin(phi)*sin(theta))/(d^2*s*(sin(b0)*cos(theta) + cos(b0)*cos(phi)*sin(theta))))

sol_vyp =

(d*s*(d*v_n*sin(b0)*cos(theta) + v_m*yp*cos(b0)*cos(theta) - v_n*zp*cos(b0)*cos(theta) + d*v_n*cos(b0)*cos(phi)*sin(theta) - v_m*yp*cos(phi)*sin(b0)*sin(theta) + v_n*zp*cos(phi)*sin(b0)*sin(theta)))/(d*sin(b0)*cos(theta) - zp*cos(b0)*cos(theta) - yp*sin(phi)*sin(theta) + d*cos(b0)*cos(phi)*sin(theta) + zp*cos(phi)*sin(b0)*sin(theta))

sol_vzp =

(d*s*(v_n*zp*sin(phi)*sin(theta) - v_m*yp*sin(phi)*sin(theta) + d*v_m*sin(b0)*cos(theta) + d*v_m*cos(b0)*cos(phi)*sin(theta)))/(d*sin(b0)*cos(theta) - zp*cos(b0)*cos(theta) - yp*sin(phi)*sin(theta) + d*cos(b0)*cos(phi)*sin(theta) + zp*cos(phi)*sin(b0)*sin(theta))

sol_vyp =

(v_n + (J13*v_zp*yp)/(J11*d^2*s))/(1/(d*s) - (J12*yp)/(J11*d^2*s))

sol_vyp =

-(d*s*(J11*d*v_n + J13*v_m*yp - J13*v_n*zp))/(J12*yp - J11*d + J13*zp)

sol_vzp =

-(d*s*(J11*d*v_m - J12*v_m*yp + J12*v_n*zp))/(J12*yp - J11*d + J13*zp)

除第一行外，结果如预期。我需要一些背景信息为什么这在python中不起作用，请。最后一个检查，我在python中也尝试了以下方法，保留矩阵符号：

In [220]: J11,J12,J13 = sympy.symbols('J11 J12 J13')

In [222]: EQ6 = v_n - v_yp/(d*s) - yp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d**2)

In [223]: EQ7 = v_m - v_zp/(d*s) - zp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d**2)

In [228]: sympy.linsolve([EQ6, EQ7], [v_yp, v_zp])
Out[228]: {(d*s*v_n, d*s*v_m)}

再一次，我不知道为什么辛普森没有给我我我期望的答案。在

---

PUSH:sympy.linsolve可以解决相对简单的系统。当sympy.solve继续运行时，它在某种程度上崩溃了。在

In [236]: a, b, c, d, e, f = sympy.symbols('a b c d e f')

In [246]: sympy.linsolve([v_n - a*v_yp - yp*(b*v_yp)], v_yp)
Out[246]: {(v_n/(a + b*yp),)}

In [247]: sympy.linsolve([v_n - a*v_yp - yp*(b*v_yp + c * v_zp)], v_yp)
Out[247]: {(v_n/a,)}

In [245]: sympy.solve([v_n - a*v_yp - yp*(b*v_yp + c * v_zp)], v_yp)
Out[245]: {v_yp: (-c*v_zp*yp + v_n)/(a + b*yp)}

因此使用solve()我可以以正确的形式检索答案。在

In [255]: sympy.linsolve([EQ6], [v_yp])
Out[255]: {(d*s*v_n,)}

In [256]: sympy.solve([EQ6], [v_yp])
Out[256]: {v_yp: (J11*d**2*s*v_n + J13*v_zp*yp)/(J11*d - J12*yp)}

In [257]: sympy.solve([EQ6, EQ7], [v_yp, v_zp])
Out[257]: 
{v_yp: d*s*(J11*d*v_n + J13*v_m*yp - J13*v_n*zp)/(J11*d - J12*yp - J13*zp),
 v_zp: d*s*(J11*d*v_m - J12*v_m*yp + J12*v_n*zp)/(J11*d - J12*yp - J13*zp)}

In [260]: sympy.solve([EQ4, EQ5], [v_yp, v_zp])
Out[260]: 
{v_yp: d*s*(d*v_n*sin(b0)*cos(theta) + d*v_n*sin(theta)*cos(b0)*cos(phi) - v_m*yp*sin(b0)*sin(theta)*cos(phi) + v_m*yp*cos(b0)*cos(theta) + v_n*zp*sin(b0)*sin(theta)*cos(phi) - v_n*zp*cos(b0)*cos(theta))/(d*sin(b0)*cos(theta) + d*sin(theta)*cos(b0)*cos(phi) - yp*sin(phi)*sin(theta) + zp*sin(b0)*sin(theta)*cos(phi) - zp*cos(b0)*cos(theta)),
 v_zp: d*s*(d*v_m*sin(b0)*cos(theta) + d*v_m*sin(theta)*cos(b0)*cos(phi) - v_m*yp*sin(phi)*sin(theta) + v_n*zp*sin(phi)*sin(theta))/(d*sin(b0)*cos(theta) + d*sin(theta)*cos(b0)*cos(phi) - yp*sin(phi)*sin(theta) + zp*sin(b0)*sin(theta)*cos(phi) - zp*cos(b0)*cos(theta))}