Question

100分，谁会用Mathematica解这个参数方程组啊!写出代码也行，能解决的话再加一百。

syms a b c x1 y1 z1 x y z;eq1=sym('cosb*cosc*x+sina*sinb*cosc*y-cosa*sinc*y+cosa*sinb*cosc*z+sina*sinc*z=x1') eq2=sym('cosb*sinc*x+sina*sinb*sinc*y+cosa*cosc*y+cosa*sinb*sinc*z-sina*cosc*z=y1')eq3=sym('-sinb*x+sina*cosb*y+cosa*cosb*z=z1')[a,b,c,]=solve(eq1,eq2,eq3,a,b,c)
这是matlab的代码转化成Mathematica的代码。。matlab解不出啊，如果有错误帮忙改一下。
其实原问题是这样的：以球心建立三维直角坐标系，球面上有一点A（x,y,z）由于某种原因球发生了绕球心的旋转，A点坐标变为（x1,y1,z1）。问题是：怎么旋转坐标系使A点在旋转后的坐标系内的坐标也是（x,y,z）. 我是想这样解决的：设绕xyz轴旋转的角度一次为abc，利用旋转矩阵可得式子（x1,y1,z1，1）=（x,y,z，1）Gx*Gy*Gz , 进而求出旋转角度abc。。。。
这样解没问题吧，可怎么就是解不出啊？？？？？？？？？？

Answer 1

不懂Matlab代码，你最好能直接给原方程。不过……也大致猜得到含义：这是关于a，b，c的三个三角方程是吧？那么，如果

1 你希望求得的是符号解而非数值解

2 你还不能对a，b，c乃至参数x, y, z, x1, y1, z1的定义域做出比较好的限制

那么，我可以负责任的告诉你，这个方程不管用什么软件都是不好解的，因为三角函数本身的周期性会使方程通解变得非常复杂，关于这一部分，你可以参考这里：

zhidao.baidu.com/question/481775259.html?oldq=1&from=evaluateTo#reply-box-1208076014

此外，除了限制未知量的定义域之外，另一个可行的方法，是适当换元，规避三角函数。一个可用的换元是万能公式：

(* 注意变形前后的a, b, c含义不同 *)
eqn = {Cos[b]*Cos[c]*x + Sin[a]*Sin[b]*Cos[c]*y - Cos[a]*Sin[c]*y + Cos[a]*Sin[b]*Cos[c]*z + Sin[a]*Sin[c]*z == x1, Cos[b]*Sin[c]*x + Sin[a]*Sin[b]*Sin[c]*y + Cos[a]*Cos[c]*y + Cos[a]*Sin[b]*Sin[c]*z - Sin[a]*Cos[c]*z == y1, -Sin[b]*x + Sin[a]*Cos[b]*y + Cos[a]*Cos[b]*z == z1} /. {Cos[a_] -> (1 - a^2)/(1 + a^2), Sin[a_] -> 2 a/(1 + a^2)}//Simplify
   
   (* {((1 + a^2) (-1 + b^2) (-1 + c^2) x + 2 (-(1 + b^2) c y + b z - b c^2 z + 2 a (b (y - c^2 y) + c z + b^2 c z) + a^2 ((1 + b^2) c y - b z + b c^2 z)))/((1 + a^2) (1 + b^2) (1 + c^2)) == x1, (-(1 + b^2) ((-1 + a^2) y + 2 a z) + (1 + b^2) c^2 ((-1 + a^2) y + 2 a z) + 2 c (-(1 + a^2) (-1 + b^2) x + 2 b (2 a y + z - a^2 z)))/((1 + a^2) (1 + b^2) (1 + c^2)) == y1, (-2 (1 + a^2) b x + 2 a y + z - a^2 z + b^2 (-2 a y - z + a^2 z))/((1 + a^2) (1 + b^2)) == z1} *)

显然，即使化成了多项式方程，这个方程组的阶数还是很高，高阶多项式（尤其是方程组）的求解也是很难的，但是，如果能够对参数的定义域也做出限制，确切地说，是给出参数的具体数值的话，这个方程的求解还是有希望的，比如：

(* 这里随机生成一组参数 *)
{x, y, z, x1, y1, z1} = RandomReal[{1, 2}, 6]

(* 本次生成的参数是：{1.90287, 1.53653, 1.61436, 1.68559, 1.18277, 1.62598} *)

(* 注意上面我已经给方程组取名为 eqn 了 *)

Solve[eqn, {a, b, c}]
(* 然后我们要换元换回来，注意这里只是简单地使用了反三角函数，要得到全部解，需要使用{Solve[Tan[d/2] == a, d], Solve[Tan[e/2] == b, e], Solve[Tan[f/2] == c, f]}/.% *)
{2 ArcTan[a], 2 ArcTan[b], 2 ArcTan[c]} /. %
(*
{{1.16575 + 18.1583 I, -1.5708 - 1.27296*10^-9 I, -2.93477 - 
   18.2374 I}, {1.16575 - 18.1583 I, -1.5708 + 
   1.27296*10^-9 I, -2.93477 + 18.2374 I}, {-1.97584 - 
   18.1583 I, -1.5708 - 1.27296*10^-9 I, 
  0.206818 + 18.2374 I}, {-1.97584 + 18.1583 I, -1.5708 + 
   1.27296*10^-9 I, 
  0.206818 - 18.2374 I}, {-0.810092 - 0.774135 I, -0.538095 - 
   35.6682 I, -0.958933 - 0.890268 I}, {-0.810092 + 
   0.774135 I, -0.538095 + 35.6682 I, -0.958933 + 
   0.890268 I}, {2.3315 + 0.774135 I, -2.50062 - 35.8224 I, 
  2.18266 + 0.890268 I}, {2.3315 - 0.774135 I, -2.50062 + 35.8224 I, 
  2.18266 - 0.890268 I}, {2.3315 + 0.774135 I, 0.353705 - 37.366 I, 
  2.18266 - 0.890268 I}, {2.3315 - 0.774135 I, 0.353705 + 37.366 I, 
  2.18266 + 0.890268 I}, {-0.810092 - 0.774135 I, 
  1.5708 - 38.4263 I, -0.958933 + 0.890268 I}, {-0.810092 + 
   0.774135 I, 
  1.5708 + 38.4263 I, -0.958933 - 0.890268 I}, {1.16575 + 18.1583 I, 
  1.5708 - 1.62229*10^-8 I, 
  1.01691 + 18.2374 I}, {1.16575 - 18.1583 I, 
  1.5708 + 1.62229*10^-8 I, 
  1.01691 - 18.2374 I}, {-1.97584 - 18.1583 I, 
  1.5708 - 1.62229*10^-8 I, -2.12468 - 18.2374 I}, {-1.97584 + 
   18.1583 I, 1.5708 + 1.62229*10^-8 I, -2.12468 + 18.2374 I}}
*)

如你所见，依旧是非常复杂。

当然了，使用Reduce求解这个方程的通解并非完全不可能。（毕竟在我写这个答案的这段时间程序一直在运行，并没有返回原式。）但是需要多久那就不知道了，你有兴趣可以试试：

(* 这是最基本的代码 *)
Reduce[{Cos[b]*Cos[c]*x + Sin[a]*Sin[b]*Cos[c]*y - Cos[a]*Sin[c]*y + 
    Cos[a]*Sin[b]*Cos[c]*z + Sin[a]*Sin[c]*z == x1,
  Cos[b]*Sin[c]*x + Sin[a]*Sin[b]*Sin[c]*y + Cos[a]*Cos[c]*y + 
    Cos[a]*Sin[b]*Sin[c]*z - Sin[a]*Cos[c]*z == y1,
  -Sin[b]*x + Sin[a]*Cos[b]*y + Cos[a]*Cos[b]*z == z1}, {a, b, c}]

追问

我是想用x1,y1，z1表示a,b,c. x y z 是常量，希望求得的是符号解而非数值解，，如果说条件的话：（x1,y1，z1）（ x ， y ， z ）都是球表面上的点球心（0,0,0）半径是330. a,b,c.的范围0~π。

追答

如果x, y, z是常量那希望能把具体值给我。

追问

sorry ,我都不知道什么是常量了，它代表数值（算是已知的），但是并不是某一个具体的数字。

追答

……那我也无能为力了。真要碰碰运气的话就把这个代码执行个十几小时试试：

Reduce[{Cos[b] Cos[c] x + Sin[a] Sin[b] Cos[c] y - Cos[a] Sin[c] y + Cos[a] Sin[b] Cos[c] z + Sin[a] Sin[c] z == x1, Cos[b] Sin[c] x + Sin[a] Sin[b] Sin[c] y + Cos[a] Cos[c] y + Cos[a] Sin[b] Sin[c] z - Sin[a] Cos[c] z == y1, -Sin[b] x + Sin[a] Cos[b] y + Cos[a] Cos[b] z == z1, x^2 + y^2 + z^2 == 330^2, x1^2 + y1^2 + z1^2 == 330^2, 0 <= a <= Pi, 0 <= b <= Pi, 0 <= c <= Pi}, {a, b, c}, Reals]

Answer 2

lz的方程似乎有问题，x, y, z, x1, y1, z1为常量，数值化后也无解（如果有的话请列出一组）
这种方程即使存在解析解，输出的式子也会非常巨大，可以占几个屏幕，又有什么意义呢？

Answer 3

如果是要解关于x,y,z的方程组,其中a,b,c,x1,y1,z1都是参数的话代码如下:(最好自行检查一下下面方程有没有抄错...我很容易抄错的...)(我只用过mathematica8和9,不确定以前版本是否可用)
s = Solve[{Cos[b] Cos[c] x + Sin[a] Sin[b] Cos[c] y - Cos[a] Sin[c] y + Cos[a] Sin[b] Cos[c] z + Sin[a] Sin[c] z == x1, Cos[b] Sin[c] x + Sin[a] Sin[b] Sin[c] y + Cos[a] Cos[c] y + Cos[a] Sin[b] Sin[c] z - Sin[a] Cos[c] z == y1, -Sin[b] x + Sin[a] Cos[b] y + Cos[a] Cos[b] z == z1}, {x, y, z}];
X = FullSimplify[s[[1, 1, 2]]]
Y = FullSimplify[s[[1, 2, 2]]]
Z = FullSimplify[s[[1, 3, 2]]]

结果:
X是-z1 Sin[b] + Cos[b] (x1 Cos[c] + y1 Sin[c])
Y是Cos[b] Cos[c] (Cos[a] Sec[b] (y1 - x1 Tan[c]) + Sin[a] (z1 Sec[c] + Tan[b] (x1 + y1 Tan[c])))
Z是Cos[a] (z1 Cos[b] + Sin[c] (y1 Sin[b] + x1 Tan[a]) + Cos[c] (x1 Sin[b] - y1 Tan[a]))

如果不是要解这个可以再描述的清楚一些.

追问

我是想用x1,y1，z1表示a,b,c. x y z 是常量，您看？