# Solving these system of 3 equations

I have a function $f:\mathbb{R}^3\rightarrow \mathbb{R}^3$ and I'd like to find (numerically) roots of $f$.

My function f is the sum of two functions involving Solve, but the timing seems very reasonable, typically:

f[3, 6, 15.] // AbsoluteTiming
(* {0.002128, {-3., -6.83896, 12.6261}} *)


For the takers, f can be imported using

<< "http://pastebin.com/raw/uA9G2uYh" (* takes a few seconds *)


What I first tried but it takes too long (without returning errors)

FindRoot[f[x, y, z], {{x, 3}, {y, 5}, {z, 6}}] (* takes too long *)
NSolve[# == 0 & /. f[x, y, z] && 0 <= x <= 5 && 0 <= y <= 10 && 0 <= z <= 5]


It seems $f$ is linear in the first variable, so to simplify the problem, we can focus on:

g[t_, s_] := f[0, t, s][[2 ;;]]


Using Manipulate, I managed to identify some good initial conditions:

g[6.246, 18.725]
Plot[Evaluate@g[6.246, s], {s, 18.6, 18.8}]
(* {0.00204707, -0.00214908} *)


but still the following command never stops:

FindRoot[g[t, s], {{t, 6.246}, {s, 18.725}}]


Minimizing the norm of g also does not work (= never ends):

FindMinimum[Norm@g[t, s], {{t, 6.246}, {s, 18.725}}]


Question: How can I find the, or some, roots of g?

The equation f[x,y,z] does not lend itself well for finding roots.

I first tried to explore the space of the function, and got some error messages (stating that the equation could not be solved, etc).

n = 3;
A2 = {{0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}, {-2,
1, 0, 0, 0, 0}, {1, -2, 1, 0, 0, 0}, {0, 0, -1, 0, 0, 0}};
lambda = {l1, 1, 1, l4, l5, 0};
ee[p_] := SparseArray[{{p, 1} -> 1, {6, 1} -> 0}] // Normal // Flatten
eqns1 = {(0.9786935821895908 a1 + 0.7848513153478486 a2 +
0.4355596199317594 a3) Cot[
0.4450418679126288 t] + (-0.15544942857434468 a1 +
0.08626792451566725 a2 + 0.19384226684174244 a3) Cot[
1.2469796037174667 t] +
0.30141660916782004 a2 Cot[1.8019377358048383 t] ==
1 + (0.13414301076393745 a1 + 0.24171735309001358 a3) Cot[
1.8019377358048383 t], (1.2204109352796038 a1 +
0.9786935821895927 a2 + 0.5431339622578357 a3) Cot[
0.4450418679126288 t] + (0.28011019135741444 a1 -
0.15544942857434632 a2 - 0.3492916954160916 a3) Cot[
1.2469796037174667 t] + (0.0596992560778067 a1 +
0.10757434232607628 a3) Cot[1.8019377358048383 t] ==
1 + 0.13414301076393761 a2 Cot[1.8019377358048383 t], a1 == 0};
lambda1[tt_] :=
Check[{a1, 1, 1, a2, a3, 0} /. Solve[eqns1 /. t -> tt] // Flatten,
Indeterminate]
eqns2 = # ==
0 & /@ (MatrixExp[s*A2].(lambda + Inverse[A2].ee[2 n]) //
FullSimplify)[[n + 1 ;;]];
lambda2[ll1_, ss_] :=
Check[({l1, 1, 1, l4, l5, 0} /. Solve[eqns2 /. s -> ss, {l4, l5}] //
N // Flatten) /. l1 -> ll1, Indeterminate]
f[l1_, t_, s_] :=
f[l1, t, s] =
Check[(lambda1[t] - lambda2[l1, s])[[{1, 4, 5}]], Indeterminate]

space = Table[{x, y, z} -> f[x, y, z], {x, 0, 5}, {y, 0, 10}, {z, 0,
5}];
Extract[space, Position[space, Rule[a_, Indeterminate]]]

(*{{0, 0, 0} -> Indeterminate, {0, 0, 1} -> Indeterminate, {0, 0, 2} ->
Indeterminate, {0, 0, 3} -> Indeterminate, {0, 0, 4} ->
Indeterminate, {0, 0, 5} -> Indeterminate, {1, 0, 0} ->
Indeterminate, {1, 0, 1} -> Indeterminate, {1, 0, 2} ->
Indeterminate, {1, 0, 3} -> Indeterminate, {1, 0, 4} ->
Indeterminate, {1, 0, 5} -> Indeterminate, {2, 0, 0} ->
Indeterminate, {2, 0, 1} -> Indeterminate, {2, 0, 2} ->
Indeterminate, {2, 0, 3} -> Indeterminate, {2, 0, 4} ->
Indeterminate, {2, 0, 5} -> Indeterminate, {3, 0, 0} ->
Indeterminate, {3, 0, 1} -> Indeterminate, {3, 0, 2} ->
Indeterminate, {3, 0, 3} -> Indeterminate, {3, 0, 4} ->
Indeterminate, {3, 0, 5} -> Indeterminate, {4, 0, 0} ->
Indeterminate, {4, 0, 1} -> Indeterminate, {4, 0, 2} ->
Indeterminate, {4, 0, 3} -> Indeterminate, {4, 0, 4} ->
Indeterminate, {4, 0, 5} -> Indeterminate, {5, 0, 0} ->
Indeterminate, {5, 0, 1} -> Indeterminate, {5, 0, 2} ->
Indeterminate, {5, 0, 3} -> Indeterminate, {5, 0, 4} ->
Indeterminate, {5, 0, 5} -> Indeterminate}*)


f[x,0,z] is indeterminate.

Let's see how the function behaves in the rest of the space.

We can take the norm of the function. Zeros should be where the norm is zero.

DensityPlot3D[Norm[f[x, y, z]], {x, 0, 5}, {y, 3, 5}, {z, 0, 5},
AxesLabel -> {x, y, z}, PlotLegends -> Automatic]


As you can see there are several planes of interest. You should focus on finding the roots within those ranges exclusively.

• Unfortunately, even if I have a good approximation of a root, I don't manage to find converge to a root. Using the norm does not help: FindMinimum[Norm@g[t, s], {{t, 6.246}, {s, 18.725}}] does not end even though it seems to be close from a root. – anderstood Sep 29 '16 at 13:30