# Find multiple solutions of a non-linear system with more variables that equations

I have a vector P0inc that depends on 7 variables that are {tinc1, tinc2, tinc3, tinc4, tinc5, tinc6, tinc7}. I would like to find let's say 20 solutions of the system P0inc==P0temp, knowing that P0temp is a known vector with numerical values. For this I am using FindInstance as shown below:

FindInstance[P0inc == P0temp, {tinc1, tinc2, tinc3, tinc4, tinc5, tinc6, tinc7}, Reals, 20];


However, it doesn't work: it takes forever to compile and never gives me any result.

I am using the following code, so that maybe this will help you understand the problem better:

1. In the 1st part in my code, I simply define some functions that will help me build the equation that I want to resolve easily:

Ti[i_, j_, b_, d_, a_, t0_, texp_] := ({{Cos[texp[[i, j]] + t0[[j]]], -Sin[texp[[i, j]] + t0[[j]]]*Cos[b[[j]]], Sin[texp[[i, j]] + t0[[j]]]*Sin[b[[j]]], a[[j]]*Cos[texp[[i, j]] + t0[[j]]]}, {Sin[texp[[i, j]] + t0[[j]]], Cos[texp[[i, j]] + t0[[j]]]*Cos[b[[j]]], -Cos[texp[[i, j]] + t0[[j]]]*Sin[b[[j]]], a[[j]]*Sin[texp[[i, j]] + t0[[j]]]}, {0, Sin[b[[j]]], Cos[b[[j]]], d[[j]]}, {0, 0, 0, 1}});
P[i_, b_, d_, a_, t0_, texp_, P7_] := Ti[i, 1, b, d, a, t0, texp].Ti[i, 2, b, d, a, t0, texp].Ti[i, 3, b, d, a, t0, texp].Ti[i, 4, b, d, a, t0, texp].Ti[i, 5, b, d, a, t0, texp].Ti[i, 6, b, d, a, t0, texp].Ti[i, 7, b, d, a, t0, texp].P7;


2. After this, I define the constant parameters of the vector P:

breal = {-90, 90, -90, -90, 90, -90, 0}*Pi/180;
dreal = {455, 0, 770, 0, -531, 0, -131};
areal = {0, -65, 0, -65, 0, -29, 0};
t0real = {0, 0, 0, 0, 0, 0, 0}*Pi/180;
P7real = {0, 0, 50, 1};


3. So now, I generate random values for the variables of the vector P and then calculate P0temp:

tall = RandomReal[360, {1, 7}]*Pi/180;
P0temp = P[1, breal, dreal, areal, t0real, tall, P7real];


4. So now, P0temp is a know vector with numerical values that corresponds to the function P[i_, b_, d_, a_, t0_, texp_, P7_] when texp==tall (the others parameters were fixed in section 2..

5. Finally now, I define the 7 parameters tinc which are the variables of the equation I want to solve and I calculate the expression of the vector P0inc (which depends on tinc) and try to solve the equation:

tinc = Table[0, {i, 1, 1}, {j, 1, 7}];
tinc[[1, 1]] = tinc1;
tinc[[1, 2]] = tinc2;
tinc[[1, 3]] = tinc3;
tinc[[1, 4]] = tinc4;
tinc[[1, 5]] = tinc5;
tinc[[1, 6]] = tinc6;
tinc[[1, 7]] = tinc7;
P0inc = P[1, breal, dreal, areal, t0real, tinc, P7real];
FindInstance[P0inc == P0temp, {tinc1, tinc2, tinc3, tinc4, tinc5, tinc6, tinc7}, Reals, 20];


As I said the FindInstance does not give me any results. It should at least give me the solution tinc==tall, as p0all was calculated from tall.

If I change the vector equation P0inc==P0temp by P0inc[[1]]==P0temp[[1]] or P0inc[[2]]==P0temp[[2]] or P0inc[[3]]==P0temp[[3]], it works! I also tried Norm[P0inc-P0temp]==0, it does not work!

Could you please test this program and see if it is working for you, and if not, any suggestions?

(PS: I am using Mathematica 10.2 and I am a beginner so my program is probably not optimized.)

• Any answers?? any suggestions?? Sep 19, 2018 at 14:50

I rewrote the code a fair bit to simplify working with it:

T[b_, d_, a_, t_] := {{Cos[t], -Cos[b] Sin[t], Sin[b] Sin[t], a Cos[t]},
{Sin[t], Cos[b] Cos[t], -Cos[t] Sin[b], a Sin[t]},
{0 b, Sin[b], Cos[b], d},
{0 b, 0 b, 0 b, ConstantArray[1, Length[b]]}};
P[b_, d_, a_, t_, p7_] := Fold[Dot, Transpose[T[b, d, a, t], {2, 3, 1}]].p7;


I rewrote Ti as T and wrote both T and P in a fashion that is more amenable to Mathematica's automatic vectorization. In many cases, you don't need to pass indices: matrices of the same size can automatically be added and multiplied piecewise, and most mathematical functions from $$\mathbb{R}\rightarrow\mathbb{R}$$ are automatically expanded over a real vector's members. Transpose is used in P to rearrange the result to match the original specification. Fold is used to eliminate a lot of code repetition: Fold[Dot, {a,b,c,d}] == a.b.c.d, which for very long lists prevents a lot of needless typing and decreases the chance of making a typo. Furthermore, it will automatically expand for longer lists as needed.

(* tall=RandomReal[360,7]*Pi/180; *)
tall = {3.802585264448382, 5.800064320152033, 2.600209644180517,
1.9888042059370301, 3.4126484372045818, 2.305359071351031,
2.8964239600958592};
breal = {-90, 90, -90, -90, 90, -90, 0}*Pi/180;
dreal = {455, 0, 770, 0, -531, 0, -131};
areal = {0, -65, 0, -65, 0, -29, 0};
t0real = {0, 0, 0, 0, 0, 0, 0}*Pi/180;
P7real = {0, 0, 50, 1};


I chose a specific tall for replicability, and I replaced it with a single vector rather than a 1x7 matrix, to correspond with the cleanup performed in P and T.

P0temp = P[breal, dreal, areal, t0real + tall, P7real]


{767.031, 284.856, 832.943, 1.}

Note that t0real and tall have been combined into one argument in P. This is intentional, and will also be useful later. They were never used individually in the original definition, regardless.

Let us now investigate P0inc. There's no need to assign the members of tinc individually, just give them spontaneously generated names:

P0inc = P[breal, dreal, areal, t0real + Table[tinc[i], {i, 1, 7}], P7real];


If you investigate the expression, it's a bit messy, but it has been expanded symbolically. The variables to find instances for are then tinc[1] through tinc[7].

As you have noted though, FindInstance takes forever in this case. So does Eliminate and Solve, which are other common methods of handling this.

So instead, we'll consider an objective function instead, of the form Norm[Rationalize[P0temp - P0inc, 0]]. Since this function is defined by an underspecified system, hopefully all of its minima are zeros, and we can then use FindMinimum to find local minima rapidly from random starting points. We use Rationalize here because we need extended precision later.

randomsoln := Quiet@FindMinimum[Norm[Rationalize[P0temp - P0inc, 0]],
Table[{tinc[i], RandomReal[{-10, 10}]}, {i, 1, 7}],
WorkingPrecision -> 40];


This code manages that. Again we use Table to generate the solution variables rather than writing them all manually, along with random initial conditions for them. It turns out that FindMinimum needs added precision for this calculation to get a reasonably small number for the objective function, so I set WorkingPrecision -> 40. If you want closer matches and are willing to spend a bit of extra time, that can easily be increased. You can then extract a random solution by evaluating randomsoln. To extract 20 of them, perhaps try:

Table[randomsoln, {i, 1, 20}]


It's entirely possible some of them will lie on the same solution line, but it's likely that all of them will be different. If you're only interested in the variables, you can use randomsoln[[2]] to discard the objective function's evaluated value and only extract the actual variable values.

Also note that tinc[7] is actually independent of this system, and this method will choose it entirely randomly.

To constrain the solution variables to the interval of $$-\pi$$ to $$\pi$$, the simplest solution is to simply Mod the values to the right place. P is periodic in tinc on all dimensions with a period of $$2\pi$$, due to the trigonometric functions used in its construction. Thus, to transform the solution values appropriately, the following will suffice:

constrainsoln[s_] := {s[[1]],
s[[2]] /. (Rule[x_, y_] :>
Rule[x, Mod[y + π, 2 π] - π])};


This is used on the solutions found by FindMinimum, e.g.:

constrainsoln[randomsoln]


The reason why FindMinimum's constraints don't always work well is because the edges of the solution space will have local minima, which FindMinimum makes no particular attempt to escape.

• Thank you very much for the solution and for taking time to optimize the code and explaining the changes. It wroks!! finally!! I have just one question: I wanted to add constraints to the variables tinc, so I did so by :randomsoln := Quiet@FindMinimum[Norm[Rationalize[P0temp - P0inc, 0]], Table[{tinc[i], RandomReal[{-Pi, Pi}], -Pi, Pi}, {i, 1, 7}], WorkingPrecision -> 40]; It works, however, the objective function's evaluated value is not close to Zero, how can I solve this please? I tried to set AccuracyGoal -> 8, but it did not help. Sep 21, 2018 at 12:28
• Probably a result of local minima in the solution space. I'm not really sure there's a good way to resolve it, other than just discarding all of the values with poor solution values. FindMinimum` isn't really a root-finding method, so it's not actually trying to get to 0. I just picked it because it's usually going to be fast. That said, because it's fairly fast, just re-picking solutions until you have enough might be an option for most value sets. I'll think about it though. Sep 21, 2018 at 15:27
• @Salma The function turns out to be periodic, so it's possible to transform arbitrary solutions to the requested solution space directly. Sep 22, 2018 at 20:37