4

我正在尝试解决 Mathemtica 中的非线性方程组。我尝试了 Solve 和 NSolve，我还尝试定义 a_{ij} 和 b_{ij} 和 m33=1 数值来简化方程，但是 Mathematica 似乎工作时间过长或者我做错了什么。在 Mathematica 我只是想找到解决方案，但我还需要一些 c/c++ 库来在我的代码中执行此操作。

“运算符”中的主要方程：

M[A[(x,y)]]=B[M[(x,y)]]

其中“运算符”是透视变换：

u= (m13 + m11*x + m12*y)/(m33 + m31*x + m32*y); 

v= (m23 + m21*x +m22*y)/(m33 + m31*x + m32*y);

我在 Mathematica 中的输入：

Solve[(b13 + (b11 (m13 + m11 x1 + m12 y1))/(m33 + m31 x1 + 
         m32 y1) + (b12 (m23 + m21 x1 + m22 y1))/(m33 + m31 x1 + 
         m32 y1))/(b33 + (b31 (m13 + m11 x1 + m12 y1))/(m33 + m31 x1 +
          m32 y1) + (b32 (m23 + m21 x1 + m22 y1))/(m33 + m31 x1 + 
         m32 y1)) == (m13 + (m11 (a13 + a11 x1 + a12 y1))/(a33 + 
         a31 x1 + a32 y1) + (m12 (a23 + a21 x1 + a22 y1))/(a33 + 
         a31 x1 + a32 y1))/(m33 + (m31 (a13 + a11 x1 + a12 y1))/(a33 +
          a31 x1 + a32 y1) + (m32 (a23 + a21 x1 + a22 y1))/(a33 + 
         a31 x1 + 
         a32 y1)) && (b23 + (b21 (m13 + m11 x1 + m12 y1))/(m33 + 
         m31 x1 + m32 y1) + (b22 (m23 + m21 x1 + m22 y1))/(m33 + 
         m31 x1 + m32 y1))/(b33 + (b31 (m13 + m11 x1 + m12 y1))/(m33 +
          m31 x1 + m32 y1) + (b32 (m23 + m21 x1 + m22 y1))/(m33 + 
         m31 x1 + 
         m32 y1)) == (m23 + (m21 (a13 + a11 x1 + a12 y1))/(a33 + 
         a31 x1 + a32 y1) + (m22 (a23 + a21 x1 + a22 y1))/(a33 + 
         a31 x1 + a32 y1))/(m33 + (m31 (a13 + a11 x1 + a12 y1))/(a33 +
          a31 x1 + a32 y1) + (m32 (a23 + a21 x1 + a22 y1))/(a33 + 
         a31 x1 + 
         a32 y1)) && (b13 + (b11 (m13 + m11 x2 + m12 y2))/(m33 + 
         m31 x2 + m32 y2) + (b12 (m23 + m21 x2 + m22 y2))/(m33 + 
         m31 x2 + m32 y2))/(b33 + (b31 (m13 + m11 x2 + m12 y2))/(m33 +
          m31 x2 + m32 y2) + (b32 (m23 + m21 x2 + m22 y2))/(m33 + 
         m31 x2 + 
         m32 y2)) == (m13 + (m11 (a13 + a11 x2 + a12 y2))/(a33 + 
         a31 x2 + a32 y2) + (m12 (a23 + a21 x2 + a22 y2))/(a33 + 
         a31 x2 + a32 y2))/(m33 + (m31 (a13 + a11 x2 + a12 y2))/(a33 +
          a31 x2 + a32 y2) + (m32 (a23 + a21 x2 + a22 y2))/(a33 + 
         a31 x2 + 
         a32 y2)) && (b23 + (b21 (m13 + m11 x2 + m12 y2))/(m33 + 
         m31 x2 + m32 y2) + (b22 (m23 + m21 x2 + m22 y2))/(m33 + 
         m31 x2 + m32 y2))/(b33 + (b31 (m13 + m11 x2 + m12 y2))/(m33 +
          m31 x2 + m32 y2) + (b32 (m23 + m21 x2 + m22 y2))/(m33 + 
         m31 x2 + 
         m32 y2)) == (m23 + (m21 (a13 + a11 x2 + a12 y2))/(a33 + 
         a31 x2 + a32 y2) + (m22 (a23 + a21 x2 + a22 y2))/(a33 + 
         a31 x2 + a32 y2))/(m33 + (m31 (a13 + a11 x2 + a12 y2))/(a33 +
          a31 x2 + a32 y2) + (m32 (a23 + a21 x2 + a22 y2))/(a33 + 
         a31 x2 + 
         a32 y2)) && (b13 + (b11 (m13 + m11 x3 + m12 y3))/(m33 + 
         m31 x3 + m32 y3) + (b12 (m23 + m21 x3 + m22 y3))/(m33 + 
         m31 x3 + m32 y3))/(b33 + (b31 (m13 + m11 x3 + m12 y3))/(m33 +
          m31 x3 + m32 y3) + (b32 (m23 + m21 x3 + m22 y3))/(m33 + 
         m31 x3 + 
         m32 y3)) == (m13 + (m11 (a13 + a11 x3 + a12 y3))/(a33 + 
         a31 x3 + a32 y3) + (m12 (a23 + a21 x3 + a22 y3))/(a33 + 
         a31 x3 + a32 y3))/(m33 + (m31 (a13 + a11 x3 + a12 y3))/(a33 +
          a31 x3 + a32 y3) + (m32 (a23 + a21 x3 + a22 y3))/(a33 + 
         a31 x3 + 
         a32 y3)) && (b23 + (b21 (m13 + m11 x3 + m12 y3))/(m33 + 
         m31 x3 + m32 y3) + (b22 (m23 + m21 x3 + m22 y3))/(m33 + 
         m31 x3 + m32 y3))/(b33 + (b31 (m13 + m11 x3 + m12 y3))/(m33 +
          m31 x3 + m32 y3) + (b32 (m23 + m21 x3 + m22 y3))/(m33 + 
         m31 x3 + 
         m32 y3)) == (m23 + (m21 (a13 + a11 x3 + a12 y3))/(a33 + 
         a31 x3 + a32 y3) + (m22 (a23 + a21 x3 + a22 y3))/(a33 + 
         a31 x3 + a32 y3))/(m33 + (m31 (a13 + a11 x3 + a12 y3))/(a33 +
          a31 x3 + a32 y3) + (m32 (a23 + a21 x3 + a22 y3))/(a33 + 
         a31 x3 + 
         a32 y3)) && (b13 + (b11 (m13 + m11 x4 + m12 y4))/(m33 + 
         m31 x4 + m32 y4) + (b12 (m23 + m21 x4 + m22 y4))/(m33 + 
         m31 x4 + m32 y4))/(b33 + (b31 (m13 + m11 x4 + m12 y4))/(m33 +
          m31 x4 + m32 y4) + (b32 (m23 + m21 x4 + m22 y4))/(m33 + 
         m31 x4 + 
         m32 y4)) == (m13 + (m11 (a13 + a11 x4 + a12 y4))/(a33 + 
         a31 x4 + a32 y4) + (m12 (a23 + a21 x4 + a22 y4))/(a33 + 
         a31 x4 + a32 y4))/(m33 + (m31 (a13 + a11 x4 + a12 y4))/(a33 +
          a31 x4 + a32 y4) + (m32 (a23 + a21 x4 + a22 y4))/(a33 + 
         a31 x4 + 
         a32 y4)) && (b23 + (b21 (m13 + m11 x4 + m12 y4))/(m33 + 
         m31 x4 + m32 y4) + (b22 (m23 + m21 x4 + m22 y4))/(m33 + 
         m31 x4 + m32 y4))/(b33 + (b31 (m13 + m11 x4 + m12 y4))/(m33 +
          m31 x4 + m32 y4) + (b32 (m23 + m21 x4 + m22 y4))/(m33 + 
         m31 x4 + 
         m32 y4)) == (m23 + (m21 (a13 + a11 x4 + a12 y4))/(a33 + 
         a31 x4 + a32 y4) + (m22 (a23 + a21 x4 + a22 y4))/(a33 + 
         a31 x4 + a32 y4))/(m33 + (m31 (a13 + a11 x4 + a12 y4))/(a33 +
          a31 x4 + a32 y4) + (m32 (a23 + a21 x4 + a22 y4))/(a33 + 
         a31 x4 + a32 y4)) && m33 == 1, {m11, m12, m13, m21, m22, m23,
   m31, m32}]
4

1 回答 1

5

为了有任何机会，您将需要实际数字来代替参数变量。否则系统太大而无法处理。

为了说明，我创建了多项式（忽略分母消失的场景），然后对参数进行随机数字替换。我本可以通过替换消除 m33，但选择将 m33-1==0 留在系统中。这样一来，任何一个方程都不需要特殊处理。出于效率考虑，可以考虑对变量中线性的方程子集进行这种消除。

In[40]:= exprs = Apply[Subtract, eqns, {1}];
e2 = Together[exprs];
polys = Numerator[e2];

In[62]:= allvars = Variables[polys];
vars = {m11, m12, m13, m21, m22, m23, m31, m32, m33};
params = Complement[allvars, vars]

Out[64]= {a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, \
b21, b22, b23, b31, b32, b33, x1, x2, x3, x4, y1, y2, y3, y4}

In[69]:= SeedRandom[11111];
substitutions = 
  Thread[params -> RandomInteger[{-1000, 1000}, Length[params]]];

In[71]:= numpolys = polys /. substitutions;

NSolve 确定系统实际上是欠定的，也就是说，您的方程具有冗余（代数相关性，更专业地说）。所以它与一个伪随机超平面相交，然后得到一个有限解集。

In[73]:= Timing[solns = NSolve[numpolys == 0, vars];]

During evaluation of In[73]:= NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with (107814 m11)/118505-(177066 m12)/118505-(164294 m13)/118505+(32943 m21)/23701+(186238 m22)/118505-(126102 m23)/118505-(178233 m31)/118505-(185338 m32)/118505+(141088 m33)/118505 == 1.

Out[73]= {357.420000, Null}

以下是本例中的解决方案。

In[74]:= solns // N

Out[74]= {{m11 -> -2.22241, m12 -> 0., m13 -> -2.41203, 
  m21 -> -0.539924, m22 -> 2.33146*10^-172, m23 -> -0.585993, 
  m31 -> 0.921382, m32 -> -4.05984*10^-172, 
  m33 -> 1.}, {m11 -> -2.22241, m12 -> 0., m13 -> -2.41203, 
  m21 -> -0.539924, m22 -> 2.33146*10^-172, m23 -> -0.585993, 
  m31 -> 0.921382, m32 -> -4.05984*10^-172, 
  m33 -> 1.}, {m11 -> -2.22241, m12 -> 0., m13 -> -2.41203, 
  m21 -> -0.539924, m22 -> 2.33146*10^-172, m23 -> -0.585993, 
  m31 -> 0.921382, m32 -> -4.05984*10^-172, 
  m33 -> 1.}, {m11 -> -0.029351, m12 -> 0., m13 -> 0.409304, 
  m21 -> 0.0182228, m22 -> -2.19075*10^-169, m23 -> -0.25412, 
  m31 -> -0.0717095, m32 -> 2.05529*10^-169, 
  m33 -> 1.}, {m11 -> -0.029351, m12 -> 0., m13 -> 0.409304, 
  m21 -> 0.0182228, m22 -> -2.19075*10^-169, m23 -> -0.25412, 
  m31 -> -0.0717095, m32 -> 2.05529*10^-169, 
  m33 -> 1.}, {m11 -> -0.029351, m12 -> 0., m13 -> 0.409304, 
  m21 -> 0.0182228, m22 -> -2.19075*10^-169, m23 -> -0.25412, 
  m31 -> -0.0717095, m32 -> 2.05529*10^-169, 
  m33 -> 1.}, {m11 -> 0.541883, m12 -> 0., m13 -> -0.123031, 
  m21 -> -4.58369, m22 -> -5.60174*10^-170, m23 -> 1.0407, 
  m31 -> -4.40445, m32 -> -5.32622*10^-170, 
  m33 -> 1.}, {m11 -> 0.541883, m12 -> 0., m13 -> -0.123031, 
  m21 -> -4.58369, m22 -> -5.60174*10^-170, m23 -> 1.0407, 
  m31 -> -4.40445, m32 -> -5.32622*10^-170, 
  m33 -> 1.}, {m11 -> 0.541883, m12 -> 0., m13 -> -0.123031, 
  m21 -> -4.58369, m22 -> -5.60174*10^-170, m23 -> 1.0407, 
  m31 -> -4.40445, m32 -> -5.32622*10^-170, m33 -> 1.}}

我们检查它们是否满足具有非常小的数值残差的原始表达式。

In[76]:= exprs /. substitutions /. solns

Out[76]= {{-4.44089*10^-16, 1.11022*10^-16, -8.88178*10^-16, 0., 
  0., -1.11022*10^-16, -4.44089*10^-16, 1.11022*10^-16, 
  0.}, {-4.44089*10^-16, 1.11022*10^-16, -8.88178*10^-16, 0., 
  0., -1.11022*10^-16, -4.44089*10^-16, 1.11022*10^-16, 
  0.}, {-4.44089*10^-16, 1.11022*10^-16, -8.88178*10^-16, 0., 
  0., -1.11022*10^-16, -4.44089*10^-16, 1.11022*10^-16, 
  0.}, {-2.22045*10^-16, 1.11022*10^-16, -1.66533*10^-16, 
  1.11022*10^-16, -5.55112*10^-17, 1.11022*10^-16, -1.11022*10^-16, 
  5.55112*10^-17, 0.}, {-2.22045*10^-16, 
  1.11022*10^-16, -1.66533*10^-16, 1.11022*10^-16, -5.55112*10^-17, 
  1.11022*10^-16, -1.11022*10^-16, 5.55112*10^-17, 
  0.}, {-2.22045*10^-16, 1.11022*10^-16, -1.66533*10^-16, 
  1.11022*10^-16, -5.55112*10^-17, 1.11022*10^-16, -1.11022*10^-16, 
  5.55112*10^-17, 0.}, {2.34535*10^-15, -1.11022*10^-15, 
  1.11022*10^-16, 2.22045*10^-16, 1.31839*10^-15, -1.11022*10^-15, 
  1.249*10^-15, -6.66134*10^-16, 
  0.}, {2.34535*10^-15, -1.11022*10^-15, 1.11022*10^-16, 
  2.22045*10^-16, 1.31839*10^-15, -1.11022*10^-15, 
  1.249*10^-15, -6.66134*10^-16, 
  0.}, {2.34535*10^-15, -1.11022*10^-15, 1.11022*10^-16, 
  2.22045*10^-16, 1.31839*10^-15, -1.11022*10^-15, 
  1.249*10^-15, -6.66134*10^-16, 0.}}
于 2012-08-02T17:54:14.003 回答