# Smoothen the result of FindRoot

I have a set of code for which it involves finding the corresponding c for each a (although I will give a value of a later) and z using the constraint toroot[a,c,z] and then substituting c back into the final expression functionS[a,z]. I also have another function for which there is a change of variable functionSR[l,z] where a->l-0.01.

d = 3;
zh = 1.5;
toroot[a_, c_?NumericQ, z_] := a - NIntegrate[(c z^(d + 1) x^d)/((1 - ((z x)/zh)^(d + 1)) (1 - c^2 (z x)^(2 d)))^(1/2), {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"]
cz[a_?NumericQ, z_?NumericQ] := c /. FindRoot[toroot[a, c, z], {c, 0.0009, 0.0000001, 10000}, WorkingPrecision -> 5]
intS[a_?NumericQ, z_?NumericQ] := NIntegrate[With[{b = z/zh}, (((-1)/(d - 1)) cz[a, z]^2 z^(2 d)) x^d ((1 - (b x)^(d + 1))/(1 - cz[a, z]^2 (z x)^(2 d)))^(1/2) - ((b^(d + 1) (d + 1))/(2 (d - 1))) x ((1 - cz[a, z]^2 (z x)^(2 d))/(1 - (b x)^(d + 1)))^(1/2) + (b^(d + 1) x)/((1 - (b x)^(d + 1)) (1 - cz[a, z]^2 (z x)^(2 d)))^(1/2)], {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"]
functionS[a_, z_] = ((-((1 - cz[a, z]^2 z^(2 d)) (1 - (z/zh)^(d + 1)))^(1/2)/(d - 1)) + intS[a, z] + 1)/(z^(d - 1));
functionSR[l_, z_] = Replace[functionS[a, z], a -> (l - 0.01), Infinity];


My problem is when I try to find the minimum of functionS[a,z] and functionSR[l,z] for some a and l, say a=1 and l=1, it gives me an error. I think it is connected to the behavior of c when a=1 or l=1.

In[23]:= FindMinimum[functionS[1, z], {z, 1.2, 1.5}] //
Quiet // AbsoluteTiming
FindMinimum[functionSR[1, z], {z, 1.2, 1.5}] // Quiet // AbsoluteTiming

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000944548 -0.00149313 I and 0.0006178735732839699 for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000944548 -0.00149313 I and 0.0006178735732839699 for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000949747 -0.00149122 I and 0.000620731102746343 for the integral and error estimates.

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

During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= General::stop: Further output of FindRoot::reged will be suppressed during this calculation.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.436961 -1.38189 I is not a real number at {z} = {1.2}.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.436961 -1.38189 I is not a real number at {z} = {1.2}.

Out[23]= {0.760891, FindMinimum[functionS[1, z], {z, 1.2, 1.5}]}

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286247 -0.0000971587 I and 0.0005426332486649041 for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286247 -0.0000971587 I and 0.0005426332486649041 for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286812 -0.0000961916 I and 0.0005442259497809905 for the integral and error estimates.

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

During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= General::stop: Further output of FindRoot::reged will be suppressed during this calculation.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.439434 -1.36539 I is not a real number at {z} = {1.2}.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.439434 -1.36539 I is not a real number at {z} = {1.2}.

Out[24]= {0.771827, FindMinimum[functionSR[1, z], {z, 1.2, 1.5}]}


For a=0.1, the plot is much more smooth

For a=1, the plot contains more bumps

Is my code badly written to extract c? Are there any changes that can be done? I have read somewhere that Reduce can also be used instead of FindRoot but I am still figuring it out. Also, is using LocalAdaptive as a method for NIntegrate suitable here?

UPDATE: Please note of the typo, I have corrected it. In the plots before, I wrote c=0.1 and c=1 but it should be a=0.1 and a=1.

The expressions of my problem is given by,

$$a = c z_s^{d+1}\int_0^1 dx \frac{x^d}{\sqrt{(1-(z_s/z_h)^{d+1} x^{d+1})(1-c^2 z_s^{2d} x^{2d})}} \tag{1}\label{1}$$

\begin{align} S &= \frac{1}{4 z_s^{d-1}}\Bigg(1 -\frac{\sqrt{(1-c^2 z_s^{2d})(1-b^{d+1})}}{d-1} - \frac{1}{d-1} c^2 z_s^{2d} \int^1_0 dx x^d \sqrt{\frac{(1-(b x)^{d+1})}{(1-c^2(z_s x)^{2d})}}\\ & -\frac{b^{d+1}(d+1)}{2(d-1)} \int^1_0 dx x \sqrt{\frac{(1-c^2(z_s x)^{2d})}{(1-(b x)^{d+1})}}\\ & + b^{d+1}\int^1_0 dx \frac{x}{\sqrt{(1-(b x)^{d+1})(1-c^2(z_s x)^{2d})}}\Bigg) \tag{2}\label{2} \end{align}

where $$b=\frac{z_s}{z_h}$$ and note that $$c=c(z_s)$$ (c=c[z]) although in the code c=c[a,z], $$c$$ should only depend on $$z_s$$ (z) since $$a$$ will be specified in the end.

Also, maybe there is a better way to design finding $$c$$. Actually, I can have another constraint where $$\frac{dS}{dz_s} = 0$$ (that is because in the end I need to minimize $$S$$ with respect to $$z_s$$) and maybe the derivative of $$\eqref{1}$$ with respect to $$z_s$$, so that these can be used to find $$c$$?

• Somehow I feel WorkingPrecision -> 5 will always be a problem, except perhaps in academic exercises. What do you think it accomplishes for you? – Michael E2 Oct 3 '20 at 20:14
• @bbgodfrey How can the only root be zero if at the plot for a=0.1 you can see that if you just look at the range 0.3 < z <1.5 the plot is smooth and not zero although small. ALSO sorry since there is a typo, I wrote c=0.1 and c=1 for the plots but it should be a=0.1 and a=1. – mathemania Oct 4 '20 at 3:30
• @MichaelE2 You are correct, I removed WorkingPrecision and still got the same plot for c, but that is assuming c is the problem. Maybe there is a better way of designing the code for c in accordance to that I need to define a function functionS[a,z] in the end. – mathemania Oct 4 '20 at 3:35
• @bbgodfrey Please see my update of the question. – mathemania Oct 4 '20 at 4:06
• My second comment, now deleted, was in error. Sorry. – bbgodfrey Oct 4 '20 at 4:17

The source of the NIntegrate error messages can be seen from a factor of the integrand, x^d/Sqrt[1-c x^d z^d], of toroot. For c > z^-3, the integrand is singular for some point in the domain, {x, 0, 1}. Moreover, if NIntegrate could integrate through the singularity (and, with help, it can), the result would be a complex number, which (presumably) is undesirable. To proceed, change the variable of integration to xd = x^(d+1) and apply the appropriate Method from here.

toroot[a_, c_?NumericQ, z_] := a - NIntegrate[((1 - xd (z /zh)^(d + 1))
(1 - c^2 xd^(2 d/(d + 1))  z^(2 d)))^(-1/2), {xd, 0, 1}, Method -> {"GlobalAdaptive",
"SingularityHandler" -> "DoubleExponential"}] (c z^(d + 1))/4


In addition, redefine cz to use the secant Method and bound the search for c to between 0 and z^-3.

cz[a_?NumericQ, z_?NumericQ] := c /.
FindRoot[toroot[a, c, z], {c, .5 z^-3, .6 z^-3/2, 0, z^-3}]


(The initial guesses, .5 z^-3 and .6 z^-3, were chosen somewhat arbitrarily.) With this definition, cz returns the correct value of c, if it exists, and z^-3 along with the FindRoot::reged error message otherwise. With these definitions, the two plots in the question can be obtained correctly as follows. For a = 1,

Plot[Check[cz[1, z], Null], {z, 1.42, zh}, AxesLabel -> {z, c},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


Check prevents plotting of a short range around c = 1.42, where no solution exists, although it does not eliminate the corresponding error messages. The second plot, for a = .1, is

LogPlot[Check[cz[.1, z], Null], {z, .2, zh}, AxesLabel -> {z, c},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]