首先，请原谅即将到来的长长的测试墙。我需要解决由坐标变换产生的以下一组方程：
给定 v_n, v_m, x_p, z_p, y_p, s, d，求解以下系统的 v_yp, v_zp：
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 的条目，可以写为：
Jp_11 = sin(theta) * cos(phi) * cos(b_0) + cos(theta) * sin(b_0)
Jp_12 = sin(theta) * sin(phi)
Jp_13 = -sin(theta) * cos(phi) * sin(b_0) + cos(theta) * cos(b_0)
b_0 也给出了。现在我想解决 v_yp 和 v_zp 的 I、II、III。我在 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]: []
由于第三个方程已经设置了 v_xp 和 v_yp,v_zp 的关系，我想仅根据 v_yp 和 v_zp 来求解系统。第 163 行没有给我想要的答案，第 162 行的输出是出乎意料的：v_xp 是朝向太阳观察者方向的速度，只有在查看圆盘中心时才应为零（假设径向速度等于 v_xp （仅在那时），设置为 v_r = 0）。
我还尝试将 v_xp (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_xp = 0 的结果。
然而，当手动进行计算时，v_yp、v_zp 的表达式变成了一些复杂的术语。如果需要，我可以稍后发布。我正在做符号计算，因为我想检查结果，看看它是否仍然可以简化。
求解方程没有按预期工作。为什么？
PUSH：我在 MATLAB 中尝试了相同的任务：这是代码和结果。请注意，为简单起见，我将 Jp_ij 替换为 Jij。
代码：
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)}
再一次，我不知道为什么 SymPy 没有给我期望的答案。
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))}