# Finding input of the function $f(a,b,c) = (d,e)$ such that $|d-e|$ = constant for $d\in[d_0,d_1]$

The actual question contains a few more parameters than mentioned in the title, but the idea will remain very similar.

Lets begin abstractly; I have a certain function $f(a,b,c) = (d,e,g,h)$. I am interested in finding the values of $a,b,c$ such that the difference $|e-g|$ is constant for the values of $e$ over a specific region $d\in[e_0,e_1]$. So basically I want to keep $g$ at a certain distance from $e$, while sweeping $e$ in the desired range. Perhaps the distance cannot be 100% constant, so it should at least fall within a certain range $(1-\epsilon)|e-g|_{D} \leq |e-g| \leq (1+\epsilon)|e-g|_{D}$ where $\epsilon$ is small and $|e-g|_{D}$ is the desired separation. The values of $d$ and $h$ do not interest me; they can be whatever they want to be to first approximation (the more complicated variant is making sure they do not get too close either).

Now, how would one do this? Simply using solve might work for simple equations, but for mine it does not. Given what kind of beast it is, I suppose one has to somehow invert it, perhaps linearlizing it in the region of choice. But perhaps I should first introduce it. Note that I understand that one wants to use a minimal example; I am afraid however that because I came to this problem from a specific example I do not know how to produce it in the absence of the entire set of equations. It is therefore a bit lengthy.

My function FitFun$(\phi,\phi_{fl1},\phi_{fl2},rules)$ is given through the function below, corresponding to the $f(a,b,c)$ from before. Lets stick to that notation in the discussion as it is simpler. (Note that I understand it is desirable to use InputForm, but somehow I cannot apply it to the code below; it gives me {})

FitFun[ϕ_, ϕfl1_, ϕfl2_, rules_] :=
Module[{ω1, ω2, ω3, ωr, J12, J13, J23, J1r, J2r, J3r},
ω1 = ω10*Sqrt[Abs@Cos[b1*(ϕ - a1) + bfl11*ϕfl1 + bfl21*ϕfl2]] /. rules;
ω2 = ω20*Sqrt[Abs@Cos[b2*(ϕ - a2) + bfl12*ϕfl1 + bfl22*ϕfl2]] /. rules;
ω3 = ω30*Sqrt[Abs@Cos[b3*(ϕ - a3) + bfl13*ϕfl1 + bfl23*ϕfl2]] /. rules;
ωr = ωr0 /. rules;
J12 = J120/Sqrt[ω10*ω20]*Sqrt[ω1*ω2] /. rules;
J13 = J130/Sqrt[ω10*ω30]*Sqrt[ω1*ω3] /. rules;
J23 = J230/Sqrt[ω20*ω30]*Sqrt[ω2*ω3] /. rules;
J1r = J1r0/Sqrt[ω10*ωr]*Sqrt[ω1*ωr] /. rules;
J2r = J2r0/Sqrt[ω20*ωr]*Sqrt[ω2*ωr] /. rules;
J3r = J3r0/Sqrt[ω30*ωr]*Sqrt[ω3*ωr] /. rules;
H = ( {
{ω1, J12, J13, J1r},
{J12, ω2, J23, J2r},
{J13, J23, ω3, J3r},
{J1r, J2r, J3r, ωr}
} ) /. rules;
eigs = Eigenvalues[H] /. rules;
Return[eigs]
]


If I now define a bunch of parameters and a set of replacement rules, the above will function

rules = {a1 -> 0.32, a2 -> 0.261, a3 -> 0.512,
b1 -> 2.59, b2 -> 2.96, b3 -> 1.75,
ω10 -> 7.068, ω20 -> 6.811, ω30 -> 7.446, ωr0 -> 12,
J120 -> 0.09, J130 -> 0.006, J230 -> 0.035,
J1r0 -> 0.0002, J2r0 -> 0.0014, J3r0 -> 0.110,
bfl11 -> 0.548, bfl12 -> -0.189, bfl13 -> -0.11,
bfl21 -> 0.06, bfl22 -> -0.63, bfl23 -> 0.05};


Lets look at it in a bit more detail in a region of interest:

    Manipulate[
Plot[Evaluate@(FitFun[ϕ, ϕfl1, ϕfl2, rules]), {ϕ, 0.9,
1.2}, PlotStyle -> {{Blue, Thickness[0.01], Opacity[.3]}, {Red,
Thickness[0.01], Opacity[.3]}, {Green, Thickness[0.01],
Opacity[.3]}, {Cyan, Thickness[0.01], Opacity[.3]}},
PlotRange -> {4, 8}], {{ϕfl1, 0}, -0.4, 0.4,
0.01}, {{ϕfl2, 0}, -0.7, 0.7, 0.01}]


This produces a plot of the form Note that I've taken the approach of using a 2D plot and setting the other two inputs as parameters, as it is a bit easier to visualize things.

Now, let us return and make it all concrete. I basically want to do the following: at $a = 1.02, b = c = 0$, the second output of the function is $e = 5.91$ and the third is $g = 5.37$; a separation of about 0.5. I now want to find the collection $[a,b,c]$ to move $e$ in the domain $[5.91,5.64]$ while keeping $|e-g| \approx 0.5$.

And honestly, with Solve computing in a never ending fashion, I have no idea how to do so. Linearize the system in some way, taylor expand about the point of interest? This is what I wish to find out.

As a final aside, in the picture above you see that the lines change color; in a way output $e$ and $g$ effectively switch when the lines cross, which makes things hard. However, if we keep the distance fixed, they will never cross and thus not switch.

• There's a handy button in the edit box to greekify your code. Makes reading an already lengthy post a bit easier. – wxffles Jul 13 '16 at 20:15

I'd take a more brute force approach than the ones you have suggested. Not because it's better, but it's just how I roll.

For a bit of speed, we can first calculate the Eigenvalues with all the parameters substituted. Then copy and paste the ones we want. They are big ugly things which I won't reproduce here.

efun[ϕ_, ϕfl1_, ϕfl2_] := Root[...];
gfun[ϕ_, ϕfl1_, ϕfl2_] := Root[...];


Now we can just do a brute force table of values (this might take 10 minutes):

spam = Table[{ϕ, ϕfl1, ϕfl2, efun[ϕ, ϕfl1, ϕfl2], gfun[ϕ, ϕfl1, ϕfl2]},
{ϕ, 0.7, 1.4, 0.0125},
{ϕfl1, -0.6, 0.6, 0.0125},
{ϕfl2, -0.9, 0.9, 0.0125}];


Now select the conditions you want:

select = Select[Flatten[spam, 2], 5.64 <= #[] <= 5.91 &&
0.49 <= Abs[#[] - #[]] <= 0.51 &];


And have a look at them:

Graphics3D[Point@select[[All, ;; 3]], Axes -> True] We appear to get a bunch of lines in $abc$-space. The next step could be to cluster the points and fit lines through them. Or find lines directly (some sort of 3D Hough transform maybe). But I couldn't find any built-in functions to do this easily.

• Very nice! I didn't have time to have a look earlier, but I'll go through it now and see if it is what I am looking for, which from the looks of it it is. – user129412 Jul 14 '16 at 9:25
• Yeah, this looks very good. The next challenge will then indeed be clustering these points, which doesn't seem trivial either, but feels like a separate question of its own. Thanks for the simple approach; it is brute indeed but it does the trick! One trivial question though, why do you use root? Does this extract the eigenvalue? I somehow can't get it to work but I'll look into it for a bit, never used the function. – user129412 Jul 14 '16 at 9:46
• @user129412 I just evaluated the Eigenvalues with all the rules, and Mathematica gives the result as a bunch of lengthy Root expressions. Then just copy and paste the two components you are interested in. It's just a crude attempt at speeding things up a bit. – wxffles Jul 14 '16 at 20:46