2
$\begingroup$

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

2

For a=1, the plot contains more bumps

1

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$?

$\endgroup$
6
  • 3
    $\begingroup$ Somehow I feel WorkingPrecision -> 5 will always be a problem, except perhaps in academic exercises. What do you think it accomplishes for you? $\endgroup$ – Michael E2 Oct 3 '20 at 20:14
  • $\begingroup$ @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. $\endgroup$ – mathemania Oct 4 '20 at 3:30
  • $\begingroup$ @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. $\endgroup$ – mathemania Oct 4 '20 at 3:35
  • $\begingroup$ @bbgodfrey Please see my update of the question. $\endgroup$ – mathemania Oct 4 '20 at 4:06
  • $\begingroup$ My second comment, now deleted, was in error. Sorry. $\endgroup$ – bbgodfrey Oct 4 '20 at 4:17
4
$\begingroup$

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}]

enter image description here

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}]

enter image description here

$\endgroup$

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.