# FindRoot gives the wrong root, but with no warning or error message

I am solving an equation set using FindRoot, as shown in the codes below.

c = 2.99792458*10^8; (*light speed*)

no[λ_] := √(2.240921 + (2.246956 (λ*10^-3)^2)/((λ*10^-3)^2 - 126.9207281) + 0.009676/((λ*10^-3)^2 - 0.01562025)); (*refractive index for o ray*)

ne[λ_] := √(2.126019 + ( 0.784404 (λ*10^-3)^2)/((λ*10^-3)^2 - 123.4034379) +  0.008578/((λ*10^-3)^2 - 0.01199134)); (*refractive index for e ray*)

ne2[θ_, λ_] := Sqrt[1/((Cos[θ])^2/(no[λ])^2 + (Sin[θ])^2/(ne[λ])^2)] ;  (*refractive index at angle θ*)

Vo[λ_] = c/( no[λ] - λ D[no[λ], λ]); (*group velocity for o ray*)

Ve[λ_] = c/( ne[λ] - λ D[ ne[λ], λ]);   (*group velocity for e ray*)

Ve2[θ_, λ_] = c/(ne2[θ, λ] - λ D[ ne2[θ, λ], λ]);    (*group velocity at angle θ*)

λi[λp_, λs_] := 1/(1/λp - 1/λs);(*idler waevelenth*)

λp = 520;
root = FindRoot[{no[λs]/ λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp,  1/Ve2[θ, λp] == 1/Vo[λs]}, {{θ, 1}, {λs, 800}}]


FindRoot gives the following root with no warning.

    {θ -> 0.772903, λs -> 794.375}


When no warning or error message appears, the root should be correct. But when it is verified in the plot, the root is false.

{no[λs]/ λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp, 1/Ve2[θ, λp] == 1/Vo[λs]} /. root

fig1 = ContourPlot[  no[λs]/ λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp, {λs, 350, 2000}, {θ, 0.1, 1.5}, ContourStyle -> {Red}, ImageSize -> {250, 200}];
fig2 = ContourPlot[ 1/Ve2[θ, λp] == 1/Vo[λs], {λs, 350, 2000}, {θ, 0.1, 1.6}, ContourStyle -> {Blue}];
Show[fig1, fig2]


The result is My question is : In Mathematica, when no warning appears, the root should be correct? Is this believe true?

P.S.

If I change λp = 500, we can obtain a true root, also with no warning. • The solution meets the precision and accuracy goals (of MachinePrecision/2). If you raise them a little, e.g.. AccuracyGoal -> 10, you get a warning. The user has to decide how close to equality the solution should give and set the appropriate goals. Equal (==) demands the two sides be closer in value than the default settings for FindRoot (see Internal\$EqualTolerance). Aug 26 '18 at 2:50
• In ContourPlot you should add in PlotPoints -> 200 to get rid of the squiggly red lines. I look forward to the experts giving the right answer but my guess is that 2 things are going on: (1) FindRoot assumes that the root exists and does not have the capability of telling explicitly that the root doesn't exist, and (2) Each side of your equations result in really small numbers and the result satisfies the default precision and accuracy goals. Blowing up the first equation by 1,000 and the second by 10^9 will get you a similar result but also with a warning.
– JimB
Aug 26 '18 at 2:52
• Many thank thanks for your valuable suggestions @MichaelE2@JimB Aug 26 '18 at 4:34

Avoiding machine precision shows that there is no intersection of the two curves for λp = 520. Consequently, FindRoot finds a point near where the two curves are closest.

c = 299792458;(*light speed*)

no[λ_] := √(2.240921 + (2.246956 (λ*10^-3)^2)/((λ*10^-3)^2 - 126.9207281) +
0.009676/((λ*10^-3)^2 - 0.01562025)) //
Rationalize[#, 0] & // Simplify // Evaluate;(*refractive index for o ray*)

ne[λ_] := √(2.126019 + (0.784404 (λ*10^-3)^2)/((λ*10^-3)^2 - 123.4034379) +
0.008578/((λ*10^-3)^2 - 0.01199134)) //
Rationalize[#, 0] & // Simplify // Evaluate;(*refractive index for e ray*)

ne2[θ_, λ_] := Sqrt[1/(Cos[θ]^2/no[λ]^2 +
Sin[θ]^2/ne[λ]^2)];(*refractive index at angle θ*)
Vo[λ_] = c/(no[λ] - λ D[no[λ], λ]);(*group velocity for o ray*)
Ve[λ_] = c/(ne[λ] - λ D[ne[λ], λ]);(*group velocity for e ray*)
Ve2[θ_, λ_] = c/(ne2[θ, λ] - λ D[ne2[θ, λ], λ]);(*group velocity at angle θ*)
λi[λp_, λs_] := 1/(1/λp - 1/λs);(*idler wavelength*)

λp = 520;

eqns = {
no[λs]/λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp,
1/Ve2[θ, λp] == 1/Vo[λs]} //
Simplify;

root = FindRoot[eqns, {{θ, 1}, {λs, 800}},
WorkingPrecision -> 16,
MaxIterations -> 1000]

(* {θ -> 0.7729018777020580, λs -> 794.3739989651066} *)

eqns /. root

(* {False, False} *)

fig1 = ContourPlot[
no[λs]/λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp,
{λs, 350, 1500}, {θ, 1/10, 6/5},
ContourStyle -> Red,
WorkingPrecision -> 16,
MaxRecursion -> 5];
fig2 = ContourPlot[
1/Ve2[θ, λp] == 1/Vo[λs],
{λs, 350, 2000}, {θ, 1/10, 8/5},
ContourStyle -> Blue,
WorkingPrecision -> 16];
Show[fig1, fig2,
Graphics[{AbsolutePointSize, Point[{λs, θ} /. root]}],
ImageSize -> {250, 200}] For λp=500, there are two roots.

λp = 500;

eqns2 = {
no[λs]/λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp,
1/Ve2[θ, λp] == 1/Vo[λs]} // Simplify;

roots = FindRoot[eqns2, {{θ, #[]}, {λs, #[]}},
WorkingPrecision -> 16,
MaxIterations -> 1000] & /@
{{7/10, 700}, {9/10, 900}}

(* {{θ -> 0.7287476137541260, λs ->
693.9615475238889}, {θ -> 0.8693491960341734, λs ->
868.7277923184990}} *)

eqns2 /. roots

(* {{True, True}, {True, True}} *)

fig3 = ContourPlot[
no[λs]/λs + ne2[θ, λi[λp, λs]]/λi[λp, λs] == ne2[θ, λp]/λp,
{λs, 350, 1500}, {θ, 1/10, 6/5},
ContourStyle -> Red,
WorkingPrecision -> 16,
MaxRecursion -> 5];
fig4 = ContourPlot[1/Ve2[θ, λp] == 1/Vo[λs],
{λs, 350, 2000}, {θ, 1/10, 8/5},
ContourStyle -> Blue,
WorkingPrecision -> 16];
Show[fig3, fig4,
Graphics[{AbsolutePointSize, Point[{λs, θ} /. roots]}],
ImageSize -> {250, 200}]
` • Many thanks to Bob Hanlon for this good answer. Aug 26 '18 at 4:36
• From this example we can know that: when no warning appears, the root may be wrong. Do you agree with this claim?@Bob Hanlon Aug 26 '18 at 4:39
• @user14634 - From this example, I would conclude that even if there is no warning, a root may be returned when one does not exist. Aug 26 '18 at 12:35
• It's not really "wrong". Numerical roots are fuzzy, loosely speaking. Think in terms of the contour plot, where the contours have positive width. Aug 26 '18 at 15:02
• Thank you! @Bob Hanlon and Daniel Lichtblau Aug 27 '18 at 15:51