まず最初に、長い試練の壁が立ちはだかることをお許しください。座標変換から生じる次の一連の方程式を解く必要があります。
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]: []
3 番目の方程式は既に v_xp と v_yp,v_zp の関係を設定しているので、v_yp と v_zp のみでシステムを解きたいと思います。行 163 は私が望む答えを与えず、行 162 の出力は予期しないものです。 (その時点でのみ)、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))}