Solve non-linear system of equation in Mathematica

Question

I'm trying to solve non-linear system of equations in Mathemtica. I tried Solve and NSolve, I also tried to define a_{ij} and b_{ij} and m33=1 numerical to simplify equation, but Mathematica seems to work too long or I doing something wrong.In Mathematica I just trying to find solution, but I also need some c/c++ lib to do this in my code.

Main equation in "operators":

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

where "operator" is perspective transform:

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

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

My input in 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}]

Answer 1

In order to have any chance with this you will need actual numbers in place of the parameter variables. Otherwise the system is too big to handle.

To illustrate, I created polynomials (ignoring denominator vanishing scenarios) and then did random numeric substitutions for the parameters. I could have eliminated m33 via replacement but opted to leave m33-1==0 in the system. That way no special handling was needed for any one equation. One might for efficiency consider doing such elimination on equation subsets that are linear in their variables.

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 decided that the system was actually underdetermined, that is, your equations have a redundancy (algebraic dependence, to put it more technically). So it intersected with a pseudorandom hyperplane and then got a finite solution set.

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}

Here are the solutions in this instance.

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.}}

We check that they satisfy the original expressions with very small numeric residuals.

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.}}