Complex values for a supposed real integral

I have an integral for which I have checked that it is supposed to be real but there some values of zm which gives a complex value.

First I checked the validity range of func1 and func2,

ContourPlot[func1[z, zm, 9.30685, -15.5, 10], {z, 0, 9.30685}, {zm, 2, 14}, Frame -> True, FrameLabel -> {"z", "zm"}, LabelStyle -> Directive[FontFamily -> "Arial", FontSize -> 20], ImageSize -> Large, PlotLegends -> Automatic] // AbsoluteTiming
ContourPlot[func2[z, zm, 9.30685, 10], {z, 9.30685, 14}, {zm, 2, 14}, Frame -> True, FrameLabel -> {"z", "zm"}, LabelStyle -> Directive[FontFamily -> "Arial", FontSize -> 20], ImageSize -> Large, PlotLegends -> Automatic] // AbsoluteTiming

Below you will see that zm~13.12... and zl=9.30685 so based on the contour plot, the function should be real numbers for the range specified.

Now here is my code,

d = 3;
ag = 20;
pg = 20;
wp = 100;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
func1[z_?NumericQ, zm_?NumericQ, zl_?NumericQ, \[Sigma]_?NumericQ, zh_?NumericQ] := Rationalize[FullSimplify[zm^d/(z^d Sqrt[f[z, zh] zm^(2 d) - f[zm, zh] z^(2 d)]) - zl^d/(z^d Sqrt[f[z, zh] (zl^(2 d) - (f[zl, zh] z^(2 d))/(f[zl, zh] + \[Sigma]^2))])], 0]
func2[z_?NumericQ, zm_?NumericQ, zl_?NumericQ, zh_?NumericQ] := Rationalize[FullSimplify[zm^d/(z^d Sqrt[f[z, zh] zm^(2 d) - f[zm, zh] z^(2 d)])], 0]
Sdiff1[zm_?NumericQ, zl_?NumericQ, \[Sigma]_?NumericQ, zh_?NumericQ] := NIntegrate[func1[z, zm, zl, \[Sigma], zh], {z, 10^-6, zl}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp] - 1/zl^(d - 1)
Sdiff2[zm_?NumericQ, zl_?NumericQ, zh_?NumericQ] := NIntegrate[func2[z, zm, zl, zh], {z, zl, zm}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp]
Sdiff[zm_?NumericQ, zl_?NumericQ, \[Sigma]_?NumericQ, zh_?NumericQ] := Sdiff1[zm, zl, \[Sigma], zh] + Sdiff2[zm, zl, zh]

I want to find the root of Sdiff[zm,zl,\[Sigma],zh]=0, so

Plot[Sdiff[zm, 0.930685, -1.55, 1], {zm, 1.31, 1.314}, PlotStyle -> {Blue, Thickness[0.005]}, ImageSize -> Large] // AbsoluteTiming
Plot[Sdiff[zm, 9.30685, -15.5, 10], {zm, 13.1, 13.14}, PlotStyle -> {Blue, Thickness[0.005]}, ImageSize -> Large] // AbsoluteTiming

I have plotted Sdiff[zm,zl,\[Sigma],zh] for different values, i.e. scale everything by factors of 10 (In this case I showed two different cases, zh=1 and zh=10). As you can see there are complex values at some points which is very weird. This will affect how I find the root using FindRoot, for example

In[311]:= Sdiff[1.3122, 0.930685, -1.55, 1]

During evaluation of In[311]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {1.10524401558531152310428738847291229469593603918764874985616942995484548973267025261496827751468857732080682244020164489570106785806413041234318401687*10^-6}. NIntegrate obtained -0.0953680951680931834745123674416812697318528782917380282484631277114279966108197976780600049228239456971990308722276004495808778188013233959933376197893 and 5.85518219686336069826959647527731082173980350941567052572749792113897389621195362705030451018833409494062810446414117340960083869273350500071564734243150.*^-6 for the integral and error estimates.

During evaluation of In[311]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[311]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {1.31219999996483775815022710925413900325446196741877373671252134559465310956494433552531673372086271474655712701563166561726401758962888297501592615831}. NIntegrate obtained 1.20335539716039122886908985026362661137726264579696310786318485276473000853231349656943135755073685036993548318482112178431183944472080815293312408373-2.98666309015455056925534269794783240056722422689357295577995834953594171518743926959095013943264159192680423190751058648637852902332826307958023336608*10^-8 I and 9.852591478099224742584715514655146848865875888877501022195950834394971266735884212697416425973180150042925051691217796179500572518840194547309207848150.*^-8 for the integral and error estimates.

Out[311]= -0.0465144 - 2.98666*10^-8 I
In[316]:=
FindRoot[Sdiff[zm, 0.930685, -1.55, 1] == 0, {zm,
1.313}] // AbsoluteTiming

During evaluation of In[316]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {1.43628442733304049176237159421375406043675406787164837826181212170174294734144693316126775482870488514544158898529087874593892511596377848099164854596*10^-6}. NIntegrate obtained -0.0953677430126576890317804757035965852145657587652392085970936658759314126174898801398971900787493739043444088160327727805498519290223977816445364470503 and 2.89397313756015414359427380459096848107793820089910650435934688803685220569080408374145792659908259535152922600074117123194677145733427270103524166053150.*^-6 for the integral and error estimates.

During evaluation of In[316]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {1.31299999996476393815843577120175756096583203717702249833126334557685736161243282560190392580959562320605806579767154678758445722585164366254519505445}. NIntegrate obtained 1.24209549585833258515400729015185658085179856431970135575634837182003541211813883750840953454444322107739595395976933113881528509394769347545102038392-1.11260715736832754271226365950881587154788934169698566310114610778037210620400003096207658071042571105498101884905957622359883099514289681467968087977*10^-7 I and 9.11306124954966805954714643395838397126221015795969449969272935181145135202368019974986545127461924277936248514998091680891708387436445255181327286198150.*^-8 for the integral and error estimates.

During evaluation of In[316]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {1.43628442733304049176237159421375406043675406787164837826181212170174294734144693316126775482870488514544158898529087874593892511596377848099164854596*10^-6}. NIntegrate obtained -0.0953677430126576890317804757035965852145657587652392085970936658759314126174898801398971900787493739043444088160327727805498519290223977816445364470503 and 2.89397313756015414359427380459096848107793820089910650435934688803685220569080408374145792659908259535152922600074117123194677145733427270103524166053150.*^-6 for the integral and error estimates.

During evaluation of In[316]:= General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation.

During evaluation of In[316]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[316]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[316]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

During evaluation of In[316]:= General::stop: Further output of NIntegrate::slwcon will be suppressed during this calculation.

During evaluation of In[316]:= NIntegrate::inumri: The integrand func2[z,1.31311 -8.39439*10^-6 I,0.930685,1] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{1.313106823851375892786563781555742025375366210937500000000000000000000000000000000000000000000000000-8.39438812464856821027488209541189689844031818211078643798828125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^-6 I,1.31310682385137589278656378154532196282846108303459917393652031300261052080212042327483089465064104398127251947970448226866065749184806172138611881869-8.39438812464856821027488209518317028336138899066660076766073388076349356073810504462990043049957127732694718006373758881957959364026202385510465424045*10^-6 I}}.

During evaluation of In[316]:= NIntegrate::inumri: The integrand func2[z,1.31311 -8.43138*10^-6 I,0.930685,1] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{1.313106794078247929746794397942721843719482421875000000000000000000000000000000000000000000000000000-8.43138406127335189024427086756929838884389027953147888183593750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^-6 I,1.31310679407824792974679439793230178198382245897726121707531775353358197619357632018964367708580121668400692114874349854204266096059225964967595336693-8.43138406127335189024427086733956372466951666984520735049820226690615498627879851667470323522444300010977212779156432486485196065538859756075619294786*10^-6 I}}.

During evaluation of In[316]:= NIntegrate::inumri: The integrand func2[z,1.31311 -8.43138*10^-6 I,0.930685,1] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{1.313106794096682072847670497139915823936462402343750000000000000000000000000000000000000000000000000-8.43137899897577064174235966786952189977455418556928634643554687500000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10^-6 I,1.31310679409668207284767049712949576220030015397972888398459406527352356028334368642572578628765041006633065497395603175327088345565398183827797449243-8.43137899897577064174235966763978737353545031654624029810073943077125217704951179703277710408782099359913267973633997017398860967033749757892941962606*10^-6 I}}.

During evaluation of In[316]:= General::stop: Further output of NIntegrate::inumri will be suppressed during this calculation.

During evaluation of In[316]:= FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

Out[316]= {18.7782, {zm -> 1.31311 - 8.43138*10^-6 I}}

It spits out complex values, which I do not know why.

Also, I still do not understand how AccuracyGoal, PrecisionGoal, and WorkingPrecision play any role in these situations after reading through many other questions here in MSE. There some who said to make the numbers exact and use Rationalize & FullSimplify, use higher WorkingPrecision, etc. None of it solved the issue.

Clear["Global*"]

d = 3;
ag = 20;
pg = 20;
wp = 100;
f[z_, zh_] := 1 - (z/zh)^(d + 1);

Rationalize function inputs prior to the calculation rather than rationalizing the end result of an inexact calculation. For performance reasons eliminate simplification except that which occurs naturally from the use of exact numbers. Use Module to localize the temporary variables.

func1[z_?NumericQ, zm_?NumericQ, zl_?NumericQ, σ_?NumericQ,
zh_?NumericQ] :=
Module[{zr, zmr, zlr, σr, zhr},
{zr, zmr, zlr, σr, zhr} =
Rationalize[{z, zm, zl, σ, zh}, 0];
zmr^d/(zr^d Sqrt[f[zr, zhr] zmr^(2 d) - f[zmr, zhr] zr^(2 d)]) -
zlr^d/(zr^d Sqrt[
f[zr,
zhr] (zlr^(2 d) - (f[zlr, zhr] zr^(2 d))/(f[zlr,
zhr] + σr^2))])];

func2[z_?NumericQ, zm_?NumericQ, zl_?NumericQ, zh_?NumericQ] :=
Module[{zr, zmr, zlr, zhr},
{zr, zmr, zlr, zhr} = Rationalize[{z, zm, zl, zh}, 0];
zmr^d/(zr^d Sqrt[f[zr, zhr] zmr^(2 d) - f[zmr, zhr] zr^(2 d)])];

Sdiff1[zm_?NumericQ, zl_?NumericQ, σ_?NumericQ, zh_?NumericQ] :=
Module[{zlr = Rationalize[zl, 0]},
NIntegrate[func1[z, zm, zlr, σ, zh], {z, 10^-6, zlr},
AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp] -
1/zlr^(d - 1)]

Sdiff2[zm_?NumericQ, zl_?NumericQ, zh_?NumericQ] :=
Module[{zlr, zmr},
{zlr, zmr} = Rationalize[{zl, zm}, 0];
NIntegrate[func2[z, zmr, zlr, zh], {z, zlr, zmr}, AccuracyGoal -> ag,
PrecisionGoal -> pg, WorkingPrecision -> wp]]

Sdiff[zm_?NumericQ, zl_?NumericQ, σ_?NumericQ, zh_?NumericQ] :=
Sdiff1[zm, zl, σ, zh] + Sdiff2[zm, zl, zh]

Test case,

Sdiff[1.3122, 0.930685, -1.55, 1] // AbsoluteTiming

(* {7.41299, \
-0.046510430706799189085784355939930933139694434574566957872349257327957297540\
0788205351038711701552046} *)

However, if the values of {ag, pg, wp} are only trying to avoid complex-valued results,

{ag = 10, pg = 10, wp = 20};

Sdiff[1.3122, 0.930685, -1.55, 1] // AbsoluteTiming

(* {1.18131, -0.04651043070679918730} *)

Plots

Plot[Sdiff[zm, 0.930685, -1.55, 1], {zm, 1.31, 1.314},
PlotStyle -> {Blue, Thickness[0.005]},
ImageSize -> Medium] // AbsoluteTiming

Plot[Sdiff[zm, 9.30685, -15.5, 10], {zm, 13.1, 13.14},
PlotStyle -> {Blue, Thickness[0.005]},
ImageSize -> Medium] // AbsoluteTiming

• Why do we need to use Module here? Also, can you give a little bit more example for Module? The documentation for it contains only one example, did not quite get it. May 31, 2021 at 7:41
• In each case there are several temporary variables used whose names could potentially conflict with names in the Global namespace. Module localizes the temporary variables to eliminate any potential naming conflicts. For additional examples search this site (Search on Mathematica ...) for Module. May 31, 2021 at 15:54
• It worked, but do you know what exactly caused the complex values to arise? In the plot there are missing lines, i.e. no real solution. May 31, 2021 at 18:22
• With your definitions machine precision inputs are fed into a complicated calculation that causes loss of precision. Sometimes, the loss of precision will cause one or more of the arguments to Sqrt to be negative and give rise to complex numbers. Rationalizing the final imprecise result of a calculation cannot create precision. However, by rationalizing the input prior to the calculation, the complicated calculation is exact and there is no loss of precision. May 31, 2021 at 19:24
• Big thanks for a very informative and straight answer! Jun 1, 2021 at 7:40