2
$\begingroup$

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

enter image description here

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. enter image description here

$\endgroup$
  • 2
    $\begingroup$ 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). $\endgroup$ – Michael E2 Aug 26 '18 at 2:50
  • 1
    $\begingroup$ 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. $\endgroup$ – JimB Aug 26 '18 at 2:52
  • $\begingroup$ Many thank thanks for your valuable suggestions @MichaelE2@JimB $\endgroup$ – user14634 Aug 26 '18 at 4:34
3
$\begingroup$

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[4], Point[{λs, θ} /. root]}],
 ImageSize -> {250, 200}]

enter image description here

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, {{θ, #[[1]]}, {λs, #[[2]]}},
    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[4], Point[{λs, θ} /. roots]}],
 ImageSize -> {250, 200}]

enter image description here

$\endgroup$
  • $\begingroup$ Many thanks to Bob Hanlon for this good answer. $\endgroup$ – user14634 Aug 26 '18 at 4:36
  • $\begingroup$ From this example we can know that: when no warning appears, the root may be wrong. Do you agree with this claim?@Bob Hanlon $\endgroup$ – user14634 Aug 26 '18 at 4:39
  • 1
    $\begingroup$ @user14634 - From this example, I would conclude that even if there is no warning, a root may be returned when one does not exist. $\endgroup$ – Bob Hanlon Aug 26 '18 at 12:35
  • 1
    $\begingroup$ It's not really "wrong". Numerical roots are fuzzy, loosely speaking. Think in terms of the contour plot, where the contours have positive width. $\endgroup$ – Daniel Lichtblau Aug 26 '18 at 15:02
  • $\begingroup$ Thank you! @Bob Hanlon and Daniel Lichtblau $\endgroup$ – user14634 Aug 27 '18 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.