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$?
WorkingPrecision -> 5will 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:14a=0.1you can see that if you just look at the range0.3 < z <1.5the plot is smooth and not zero although small. ALSO sorry since there is a typo, I wrotec=0.1andc=1for the plots but it should bea=0.1anda=1. $\endgroup$ – mathemania Oct 4 '20 at 3:30WorkingPrecisionand still got the same plot forc, but that is assumingcis the problem. Maybe there is a better way of designing the code forcin accordance to that I need to define a functionfunctionS[a,z]in the end. $\endgroup$ – mathemania Oct 4 '20 at 3:35