Avoid imaginary result FindRoot

Question

I have a code where I'm using FindRoot for the equation toroottp. I plug in some sample values on zs,zQ,sig,zh and plot toroot to get a feeling of its range.

I have tried using tp[9,9.30685,-15.50169,10] but it produces some errors. I know the error occurs because clearly near zs=1.1 the plot diverges (the plot is symmetric so there are two roots, one in the positive the other in the negative but choosing the positive wil suffice). The range I put for FindRoot say, {tp, 1, 0.8, 1.15} produces an error, but if I reduce the upper limit a little bit, say, {tp, 1, 0.8, 1.12} then it produces no error.

The question is how do I set up FindRoot so that it will not produce an error even if my proposed range goes beyond those that converge. I want to resolve this since I will change my parameters later and use this for other commands, this is a problem since I do not know beforehand where exactly is the root and where it starts to diverge as in the plot. You cannot expect me to check one point at a time just to make sure there is no error, that will be adsurd, but I guess any improvement on my setup will still be beneficial. The major problem here is the upper limit. On the other hand, I think no matter how I change the parameters the plots will have the same form so resolving the upper limit issue wil suffice I guess. Any guidance?

Just a note, I will vary zs say in the range [7,9.9]. The rest are fixed.

d = 3;
ag = 10;
pg = 10;
wp = 40;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
toroottp[tp_?NumericQ, zs_?NumericQ, zQ_?NumericQ, sig_?NumericQ, zh_?NumericQ] := Module[{tpr, zsr, zQr, sigr, zhr}, {tpr, zsr, zQr, sigr, zhr} = Rationalize[{tp, zs, zQ, sig, zh}, 0]; NIntegrate[SetPrecision[(1/z^d) (Sqrt[1/(f[z, zhr] (1 + (((f[zsr, zhr]^2 tpr^2)/f[z, zhr]) - 1) (z/zsr)^(2 d) (1/(1 - f[zsr, zhr] tpr^2))))] - (Sqrt[zQr^(2 d)/(f[z, zhr] (zQr^(2 d) - z^(2 d) (1 + sigr^2 f[zQr, zhr]^-1)^-1))])), wp], {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 300] - NIntegrate[SetPrecision[(1/z^d) (Sqrt[zQr^(2 d)/(f[z, zhr] (zQr^(2 d) - z^(2 d) (1 + sigr^2 f[zQr, zhr]^-1)^-1))]), wp], {z, zsr, zQr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 300] - 1/zQr^(d - 1)]
tp[zs_?NumericQ, zQ_?NumericQ, sig_?NumericQ, zh_?NumericQ] := tp /. FindRoot[SetPrecision[toroottp[tp, zs, zQ, sig, zh] == 0, wp], {tp, 1, 0.8, 1.15}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 200]

Plot[toroottp[tp, 9, 9.30685, -15.50169, 10], {tp, -1.2, 1.2}, PlotStyle -> {Blue, Thickness[0.005]}, PlotRange -> Full, ImageSize -> Large]

tp[9,9.30685,-15.50169,10]
(*During evaluation of In[1315]:= 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[1315]:= 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[1315]:= 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[1315]:= General::stop: Further output of NIntegrate::slwcon will be suppressed during this calculation.
1.111809133761684722509615030076769149700*)

Answer 1

It helps to mark out the singularities or eliminate them. We can get rid of one with a substitution z -> -u^2 + 9; another one appears only for certain values of tp/tpr, so we solve for that one. You wouldn't need to solve for it, if tp stayed in the approximate range {tp, -1.1313, 1.1313}. Omitting Solve saves only 5% of the execution time at most, though. Avoiding the singular point, one can save about as much again by using a higher-order integration rule, Method -> {"GaussKronrodRule", "Points" -> 11} (81.7 sec. vs. 88.7 sec. in the OP's Plot but over {tp, -1.1313, 1.1313}).

d = 3;
ag = 16;
pg = 6; (* N.B. *)
wp = 40;
f[z_, zh_] := 1 - (z/zh)^(d + 1);

toroottp[tp_?NumericQ, zs_?NumericQ, zQ_?NumericQ, sig_?NumericQ, 
   zh_?NumericQ] :=
  Module[{tpr, zsr, zQr, sigr, zhr, integrand, singpts},
   {tpr, zsr, zQr, sigr, zhr} = 
    Rationalize[{tp, zs, zQ, sig, zh}, 0];
   
   integrand = ((1/
          z^d) (Sqrt[
           1/(f[z, 
               zhr] (1 + (((f[zsr, zhr]^2 tpr^2)/f[z, zhr]) - 
                   1) (z/
                    zsr)^(2 d) (1/(1 - f[zsr, zhr] tpr^2))))] - (Sqrt[
            zQr^(2 d)/(f[z, 
                zhr] (zQr^(2 d) - 
                 z^(2 d) (1 + sigr^2 f[zQr, zhr]^-1)^-1))])) Dt[z, 
         u] /. z -> -u^2 + 9 // Simplify[#, u > 0] &);
   singpts = 
    Reverse@
     Union@
      Flatten[
       SolveValues[# && Sqrt[9 - zsr] < u < 3, u, Reals] & /@ 
        Cases[FunctionSingularities[integrand, u], _Equal]];
   
   NIntegrate[SetPrecision[integrand, wp],
     (*{z,0,zsr}*){u, 3, Sequence @@ singpts, Sqrt[9 - zsr]},
     AccuracyGoal -> ag, PrecisionGoal -> pg,
     WorkingPrecision -> wp, MaxRecursion -> 30] - 
    NIntegrate[SetPrecision[
      (1/z^d) (Sqrt[
         zQr^(2 d)/(f[z, 
             zhr] (zQr^(2 d) - 
              z^(2 d) (1 + sigr^2 f[zQr, zhr]^-1)^-1))]), wp], {z, 
      zsr, zQr},
     AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, 
     MaxRecursion -> 30] - 1/zQr^(d - 1)];

tp[zs_?NumericQ, zQ_?NumericQ, sig_?NumericQ, zh_?NumericQ] := 
 tp /. FindRoot[
   SetPrecision[toroottp[tp, zs, zQ, sig, zh] == 0, wp], {tp, 1, 0.8, 
    1.15}, AccuracyGoal -> ag, PrecisionGoal -> pg, 
   WorkingPrecision -> wp, MaxIterations -> 200];