# Need help getting the intercept of two lines in a ContourPlot

I am trying to find a set of parameters $\{\epsilon,\epsilon_b\}$, so that the two local minima have exactly the difference of 15.1 (see code below) and the left sided local minimum and the local maximum have a distance of 2.5 (see code below). I am able to solve this problem in form of two ContourPlots and Show, but I am not able to find the intercept of the two ContourLines with NSolve. Maybe someone can help me with that and explain why I can solve this problem graphically, but not numerically.

The code:

U = Function[{r, ϵ, ϵb},E^(-((2(r-2σ)^2)/[Sigma]^2)) ϵb + 4 ϵ (-(σ^6/r^6) + σ^12/r^12) + (lc Log[(lc (1 - r^2/lc^2) Sinh[(r (3 - r^2/lc^2))/(lc (1 - r^2/lc^2))])/(r (3 - r^2/lc^2))])/(β σ)];

T = 298.15;
kB = 8.3145*10^-3;
β = 1/(kB T);
lc = 65;
σ = 2.5;

xfCoord[ϵ_?NumberQ, ϵb_?NumberQ] := FindMinimum[{U[r, ϵ, ϵb],0 < r < 2 σ}, {r, σ}]
xf[ϵ_?NumberQ, ϵb_?NumberQ] := r /. FindMinimum[{U[r, ϵ, ϵb],0 < r < 2 σ}, {r, σ}][[2]]
Uf[ϵ_?NumberQ, ϵb_?NumberQ] := FindMinimum[{U[r, ϵ, ϵb],0 < r < 2 σ}, {r, σ}][[1]]
xuCoord[ϵ_?NumberQ, ϵb_?NumberQ] :=FindMinimum[{U[r, ϵ, ϵb],3 σ < r < lc}, {r,3 σ}]
xu[ϵ_?NumberQ, ϵb_?NumberQ] :=r /. FindMinimum[{U[r, ϵ, ϵb],3 σ < r < lc}, {r, 3 σ}][[2]]
Uu[ϵ_?NumberQ, ϵb_?NumberQ] := FindMinimum[{U[r, ϵ, ϵb], 3 σ < r < lc}, {r,3 σ}][[1]]
xbCoord[ϵ_?NumberQ, ϵb_?NumberQ] :=FindMaximum[{U[r,ϵ, ϵb], σ < r <0.9 lc}, {r, 2 σ}]
xb[ϵ_?NumberQ, ϵb_?NumberQ] :=r /. FindMaximum[{U[r,ϵ, ϵb], σ < r <0.9 lc}, {r, 2 σ}][[2]]
Ub[ϵ_?NumberQ, ϵb_?NumberQ] :=FindMaximum[{U[r, ϵ, ϵb], σ < r <0.9 lc}, {r, 2 σ}][[1]]

Δx[ϵ_?NumberQ, ϵb_?NumberQ] := xb[ϵ, ϵb] - xf[ϵ, ϵb]
ΔG[ϵ_?NumberQ, ϵb_?NumberQ] := Uu[ϵ, ϵb] - Uf[ϵ, ϵb]

p1 = ContourPlot[Δx[ϵ, ϵb] == 2.5, {ϵ, 2, 20}, {ϵb, 1, 20},FrameLabel -> {"ϵ","ϵb"}];
p2 = ContourPlot[ΔG[ϵ, ϵb] == 15.1, {ϵ, 2, 20}, {ϵb, 1, 20},FrameLabel -> {"ϵ","ϵb"},ContourStyle -> Green];
Show[p1, p2]


This now does not work:

NSolve[{Δx[ϵ, ϵb] == 2.5, ΔG[ϵ, ϵb] == 15.1}, {ϵ, ϵb}]


It just gives back the same expression, meaning it can not evaluate it, I guess.

• I'm not sure you gave us all your code. I'm seeing [Sigma] in your definition for U. What is [Sigma]? Commented Oct 14, 2015 at 19:06
• $\sigma = 2.5$ (defined in the second block of the code above). Commented Oct 14, 2015 at 19:53
• Then the slash in front of [sigma] should be \ instead of /. So \[sigma]. Commented Oct 14, 2015 at 20:14

That's because NSolve[ ] isn't the right tool for the job as it deals primarily with linear and polynomial equations.

Try FindRoot[ ] instead:

FindRoot[{Δx[ϵ, ϵb] == 2.5, ΔG[ϵ, ϵb] == 15.1}, {{ϵ, 5}, {ϵb, 5}}]

(* {ϵ -> 14.6863, ϵb -> 5.28574} *)

U = Function[{r, ϵ, ϵb},
E^(-((2 (r - 2 σ)^2)/σ^2)) ϵb +
4 ϵ (-(σ^6/r^6) + σ^12/
r^12) + (lc Log[(lc (1 -
r^2/lc^2) Sinh[(r (3 - r^2/lc^2))/(lc (1 - r^2/lc^2))])/(r (3 -
r^2/lc^2))])/(β σ)];

T = 298.15;
kB = 8.3145*10^-3;
β = 1/(kB T);
lc = 65;
σ = 2.5;


To restrict a function's arguments to numeric values use NumericQ rather than NumberQ. NumberQ will not accept numeric symbols (e.g., π, E, GoldenRatio).

NumberQ /@ {π, E, GoldenRatio}

(*  {False, False, False}  *)

xfCoord[ϵ_?NumericQ, ϵb_?NumericQ] :=
FindMinimum[
{U[r, ϵ, ϵb], 0 < r < 2 σ},
{r, σ}];
xf[ϵ_?NumericQ, ϵb_?NumericQ] :=
r /. FindMinimum[{U[r, ϵ, ϵb],
0 < r < 2 σ}, {r, σ}][[2]];
Uf[ϵ_?NumericQ, ϵb_?NumericQ] :=
FindMinimum[
{U[r, ϵ, ϵb], 0 < r < 2 σ},
{r, σ}][[1]];
xuCoord[ϵ_?NumericQ, ϵb_?NumericQ] :=
FindMinimum[{U[r, ϵ, ϵb], 3 σ < r < lc},
{r, 3 σ}];
xu[ϵ_?NumericQ, ϵb_?NumericQ] :=
r /. FindMinimum[
{U[r, ϵ, ϵb], 3 σ < r < lc},
{r, 3 σ}][[2]];
Uu[ϵ_?NumericQ, ϵb_?NumericQ] :=
FindMinimum[
{U[r, ϵ, ϵb], 3 σ < r < lc},
{r, 3 σ}][[1]];
xbCoord[ϵ_?NumericQ, ϵb_?NumericQ] :=
FindMaximum[
{U[r, ϵ, ϵb], σ < r < 0.9 lc},
{r, 2 σ}];
xb[ϵ_?NumericQ, ϵb_?NumericQ] :=
r /. FindMaximum[
{U[r, ϵ, ϵb], σ < r < 0.9 lc},
{r, 2 σ}][[2]];
Ub[ϵ_?NumericQ, ϵb_?NumericQ] :=
FindMaximum[
{U[r, ϵ, ϵb], σ < r < 0.9 lc},
{r, 2 σ}][[1]];

Δx[ϵ_?NumericQ, ϵb_?NumericQ] :=
xb[ϵ, ϵb] - xf[ϵ, ϵb];
ΔG[ϵ_?NumericQ, ϵb_?NumericQ] :=
Uu[ϵ, ϵb] - Uf[ϵ, ϵb];

ContourPlot[
{Δx[ϵ, ϵb] ==
2.5, ΔG[ϵ, ϵb] == 15.1},
{ϵ, 2, 20}, {ϵb, 1, 20},
FrameLabel -> {"ϵ", "ϵb"}]


Use FindRoot. The initial estimates can be taken from the plot

FindRoot[{Δx[ϵ, ϵb] ==
2.5, ΔG[ϵ, ϵb] == 15.1},
{{ϵ, 15}, {ϵb, 5}}]

(*  {ϵ -> 14.6863, ϵb -> 5.28574}  *)