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

You're using a program to generate that Mathematica input, right? ;) — irrelephant, Jul 31 '12 at 07:28
I do substitutions with Mathematica to get two equations and then copy-paste it 4 times to get 8 equations for 4 points (x1,y1),...,(x4,y4). — mrgloom, Jul 31 '12 at 07:42
@mrgloom Could you please clarify where are `u,v` entering the top equation ? — b.gatessucks, Aug 01 '12 at 13:35
@b.gatessucks u,v just result of operation, for example (u,v)=A[(x,y)] and then we substitut (u,v) into M and get M[(u,v)]. A,B,M are the same operators, just different constants a(ij),b(ij),m(ij). — mrgloom, Aug 01 '12 at 14:03
It appears to be overdetermined: 9 eqns, 8 unknowns. Not surprisingly, a random substitution for a11 et al gives an empty solution set. — Daniel Lichtblau, Aug 02 '12 at 16:36
Okay, I guess m33 was intended to be a variable, so it is not overdetermined. — Daniel Lichtblau, Aug 02 '12 at 16:52
This question appears to be off-topic because it is about Mathematica and there is a specific StackExchange site for that. — Saullo G. P. Castro, Jan 19 '14 at 13:52

score 5 · Accepted Answer · answered Aug 02 '12 at 17:54

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

excuse me, i can't understand it.Please tell me at least how technic that you use called? Also I can give real numbers for a(i,j) b(i,j), so we can test solution. And here http://math.stackexchange.com/questions/177075/solve-non-linear-system-of-equation-in-mathematica people say that this non-linear equation can be solved in 3D space as matrix equation, maybe it will be simpler to solve it in this way? — mrgloom, Aug 03 '12 at 06:48
for example A matrix (0 -1 300 1 0 0 0 0 1) B matrix (-0.4009 -1.0787 446.1463 1.6180 0.8875 -159.2272 0.0003 0.0029 1) but it can be a problem that matrix B was finded with some error. — mrgloom, Aug 03 '12 at 08:06
My suggestion based on the math.stackexchange response would be to try to reframe this as a matrix null space computation. Using a nonzero tolerance you might be able to work past small numerical error in the input parameters. — Daniel Lichtblau, Aug 04 '12 at 02:48
can you be more specific? how matrix null space can help me? — mrgloom, Aug 06 '12 at 06:25
If you have a matrix eqn of the form A.M = M.B with A, B known and M a matrix of unknowns, then A.M-M.B = 0 and this can be reformulated as a matrix-times-vector = 0, hence a null space computation. — Daniel Lichtblau, Aug 06 '12 at 14:33
I find out that there is to method to test if non-trivial solution exist `NullSpace[KroneckerProduct[IdentityMatrix[n], A]-KroneckerProduct[B, IdentityMatrix[n]]]` and `resultant(det(A−λE),det(B−λE),λ)=0` but it seems mathematica give the trivial solution for this test method if I put matrices with parameters. — mrgloom, Aug 07 '12 at 06:12

Solve non-linear system of equation in Mathematica

1 Answers1