math - 在 Mathematica 中求解非线性方程组

标签 math wolfram-mathematica numerical-methods numerical nonlinear-functions

我正在尝试解决 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}]

最佳答案

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

为了说明这一点，我创建了多项式(忽略分母消失的场景)，然后对参数进行了随机数字替换。我本可以通过替换消除 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.}}

关于math - 在 Mathematica 中求解非线性方程组，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/11735078/

上一篇：r - 试图将公式应用于 R 中的每一列，如何将数据提供给公式？

下一篇：ruby-on-rails - Formtastic 复选框提交带有空白条目的列表(除了选择的 ID)

相关文章：

floating-point - 正确舍入两个处理溢出的 float 之和的 sqrt 计算

string - 使 Mathematica 中的字符串操作更加方便

performance - 从给定的数字和运算集创建表达式树，并在 Mathematica 8 或更高版本中找到计算结果为目标数字的表达式树

c++ - 如何创建包含(人工生成的)高斯(正态)分布的 vector ？

c++ - 如何在 C++ 中同时为整数、 float 和 double 据类型重载运算符

wolfram-mathematica - 在 Mathematica 中使用给定颜色设置给定 ListPlot 的所有点

python - python中平流方程四阶龙格-库塔编程

c - 另一个快速三角函数

language-agnostic - 编程中幺半群/半群的例子

java - 从文本文件读取乐透

©2024 IT工具网  联系我们