1
$\begingroup$

I have two expressions given by,

$$a = \int_0^1 dx \frac{c z_s^{d+1} 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},$$

$$S = \int_\epsilon^{z_s} dx \frac{1}{x^d} \sqrt{\frac{1}{(1-(x/z_h)^{d+1})(1-c^2 x^{2d})}}\tag{2}\label{2},$$

where $c=c(z_s)$ and $a$ is a function of $z_s$, and $z_h$ is a constant which can be given a value later. Also, even though $a$ is a function of $z_s$ in the expression, in the final expression for $S$ I will just give it a value. $d$ is just a dimension parameter which can be given a value say, $d=3$.

The goal is to find an expression of $S$ such that $c$ is eliminated. Now, one condition which can help achieve that is,

$$\frac{dS}{dz_s} = 0\tag{3}\label{3},$$

actually I need to find an expression for the minimal $S$ (with respect to $z_s$) to which $c$ dictates.

I have asked some related problems before but for eq. $\eqref{2}$ in those cases I made a transformation $x=z_sy$ so that the integral limits there is $[\epsilon/z_s,1]$ and then expand out the finite and divergent terms then I just showed the finite terms there so that I allow the lower limit to be 0.

Solving a function inside an integral

Smoothen the result of FindRoot

In this case, I just start with the original integral $\eqref{2}$ and set a small number for $\epsilon$ (this is just a modification of my previous problem in order to eliminate the divergence contribution but the result should not differ for those regions far from the divergence region).

The problem is when I try to plot those same quantities I get a different result.

The previous code,

d = 3;
zh = 1.5;
torootL[a_, c_?NumericQ, z_] := a - NIntegrate[(c z^(d + 1) xd^(d/(d + 1)))/((1 - xd (z/zh)^(d + 1)) (1 - xd^((2 d)/(d + 1)) c^2 (z)^(2 d)))^(1/2), {xd, 0, 1}, MaxRecursion -> 12, PrecisionGoal -> 6, Method -> {"GlobalAdaptive", "SingularityHandler" -> "DoubleExponential"}]
cz[a_?NumericQ, z_?NumericQ] := c /. FindRoot[torootL[a, c, z], {c, 0.001, 0.0000001, 1000}]
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 -> 12, PrecisionGoal -> 6, Method -> {"GlobalAdaptive", "SingularityHandler" -> "DoubleExponential"}]
functionS[a_, z_] = ((-((1 - cz[a, z]^2 z^(2 d)) (1 - (z/zh)^(d + 1)))^(1/2)/(d - 1)) + intS[a, z] + 1)/(4 z^(d - 1));
Plot[functionS[0.1, z], {z, 1, zh}, PlotPoints -> 8, AxesLabel -> {"z_s", "S"}, AxesStyle -> Directive[Black, 18], PlotStyle -> Thick, PlotRange -> Full, ImageSize -> Large] // Quiet // AbsoluteTiming

produces this plot,

pic1

The present code where I just changed the lines corresponding to the discussion I made above,

d = 3;
zh = 1.5;
torootL[a_, c_?NumericQ, z_] := a - NIntegrate[(c z^(d + 1) xd^(d/(d + 1)))/((1 - xd (z/zh)^(d + 1)) (1 - xd^((2 d)/(d + 1)) c^2 (z)^(2 d)))^(1/2), {xd, 0, 1}, MaxRecursion -> 12, PrecisionGoal -> 6, Method -> {"GlobalAdaptive", "SingularityHandler" -> "DoubleExponential"}]
cz[a_?NumericQ, z_?NumericQ] := c /. FindRoot[torootL[a, c, z], {c, 0.001, 0.0000001, 1000}]

intS[a_?NumericQ, z_?NumericQ] := NIntegrate[(1/x^d) (1/((1 - (x/zh)^(d + 1)) (1 - cz[a, z]^2 x^(2 d))))^(1/2), {x, 10^-10, z}, MaxRecursion -> 20, PrecisionGoal -> 6, Method -> {"GlobalAdaptive"}]
functionS[a_, z_] = intS[a, z] + 1/(z^(d - 1));
Plot[functionS[0.1, z], {z, 1, zh}, PlotPoints -> 8, AxesLabel -> {"z_s", "S"}, AxesStyle -> Directive[Black, 18], PlotStyle -> Thick, PlotRange -> Automatic, ImageSize -> Large] // Quiet // AbsoluteTiming

produces the plot,

pic2

where the Method I used in the present is only GlobalAdaptive since if I use the Method in the old code it produces a flat line which I don't know why. Anyways, the old and present codes should produce the same plot especially in those region close to $z_s = 1.5$ to which in the first image you can see a minimum occurring around $z_s \approx 1.4$.

Why is this so? I'm still trying to understand the NIntegrate Integration Strategies documentation although it is hard to understand for me since I'm not a heavy user maybe until now.

I was thinking of maybe there is a different strategy to eliminate $c$ in the expression for $S$ using $\eqref{1}$ $\eqref{2}$ $\eqref{3}$ at the same time. Any comments?

UPDATE:

I noticed that when I changed the lower limit of the integral in intS I can get the same curve but with different magnitude so I think that cannot be correct.

intS[a_?NumericQ, z_?NumericQ] := NIntegrate[(1/x^d) (1/((1 - (x/zh)^(d + 1)) (1 - cz[a, z]^2 x^(2 d))))^(1/2), {x, 10^-2, z}, MaxRecursion -> 20, PrecisionGoal -> 6, Method -> {"GlobalAdaptive"}]

pic1

intS[a_?NumericQ, z_?NumericQ] := NIntegrate[(1/x^d) (1/((1 - (x/zh)^(d + 1)) (1 - cz[a, z]^2 x^(2 d))))^(1/2), {x, 1, z}, MaxRecursion -> 20, PrecisionGoal -> 6, Method -> {"GlobalAdaptive"}]

pic2

You can see that when I changed the lower limit to 10^-2 the curve already appeared but the magnitude is wrong, and when I set it to 1 the curve is there and the magnitude is closer to the original. However, I think I cannot just tune the lower limit, it must be the case that $\epsilon \leq z_s \leq z_h$ since $\epsilon$ is a regulator. Note that since that is the case, $z_s$ should approach $\epsilon$ first then take the limit where $\epsilon \rightarrow 0$. So in the case of Mathematica I can just set $\epsilon$ to the machine-wise zero, I think $10^{-10}$ is considered to be 0 but not really setting 0.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ It is true that $\frac{d a}{d z_s} = 0$ ? $\endgroup$
    – Cesareo
    Commented Oct 13, 2020 at 9:52
  • $\begingroup$ @Cesareo Good question, I'm thinking that $a$ should depend on $z_s$ as seen from $\eqref{1}$, although in the end I can just give it a value. I think I should change it to $a=a(z_s)$. $\endgroup$
    – mathemania
    Commented Oct 13, 2020 at 9:58
  • 1
    $\begingroup$ It looks like this problem has no unique solution. May be we need to add some restrictions? $\endgroup$
    – Alex Trounev
    Commented Oct 16, 2020 at 11:49
  • $\begingroup$ @AlexTrounev The only restrictions in the problem are that $S$ must be a minimum for some $z_s$ and $c$ is constrained by this fact. $\endgroup$
    – mathemania
    Commented Oct 18, 2020 at 14:42
  • $\begingroup$ @mathemania Are $a, c, z_s$ real? If $a, c, z_s$ are real then we have some restrictions for $c, z_s$ due to convergence of integral. $\endgroup$
    – Alex Trounev
    Commented Oct 18, 2020 at 14:54

1 Answer 1

Reset to default
1
$\begingroup$

If all parameters $a,z_s,z_h,c$ are real, then we have constraints dictated by existence of integral $a(z_s,z_h,c)$. As an example let put $d=3, z_h=3/2, a=1/10$, from this data we can compute $c(z_s)$ as follows (we made substitution $y=x^{d+1})$:

d = 3; zh = 3/2;
a[zs_?NumericQ, c_?NumericQ] := 
 c zs^(d + 1)/(d + 1) NIntegrate[
   Sqrt[1/(1 - (zs/zh)^(d + 1) y)/(1 - 
       c^2 zs^(2 d) y^(2 d/(d + 1)))], {y, 0, 1}]

Since a[] is real, we can plot this function with using constrains for existence of integral

ContourPlot[(a[zs, c] - 1/10) Boole[c^2 zs^(2 d) <= 1] == 0, {zs, 0, 
  zh}, {c, 0, 1}, PlotPoints -> 50]

Figure 1 The line on this picture is function $c(z_s)$ we got from the existence of integral. Now we can calculate this function directly. Let check that

 cc[zs_] := 
 c /. FindRoot[(a[zs, c] - 1/10) Boole[c^2 zs^(2 d) <= 1] == 0, {c, 0,
     1}, Method -> "Brent"]

 cc[.781598]

Out[]= 0.99999

This is starting point we use in a loop

zs = .79; 
cr[0] = cc[zs]; Do[zs = zs + .01; 
 cr[i] = c /. 
   FindRoot[(a[zs, c] - 1/10) Boole[c^2 zs^(2 d) <= 1] == 0, {c, 
     cr[i - 1]}];, {i, 1, (zh - 79/100) 100}]

lst1 = Table[{(79 + i) .01, cr[i]}, {i, 0, 150 - 79}];

Plot this list and compare with Figure1 above

ListPlot[Join[{{.781598, 1}}, lst1], 
 AxesLabel -> {"\!\(\*SubscriptBox[\(z\), \(s\)]\)", "c"}]

Figure 2 Now we can interpolate data and use it to compute integral $S(z_s)$

fc = Interpolation[Join[{{.781598, 1}}, lst1]];

S[zs_?NumericQ] := 
 NIntegrate[
  x^-d/Sqrt[(1 - (x/zh)^(d + 1)) (1 - fc[zs]^2 x^(2 d))], {x, .782, 
   zs}]

Finally we plot $S(z_s)$ and check that there is no local minimum with $S'(z_s)=0$. Also pay attention that due to $c(z_s)$ definition we define $\epsilon $ as $\epsilon =0.782$

Plot[S[zs], {zs, 0.782, zh}, 
 AxesLabel -> {"\!\(\*SubscriptBox[\(z\), \(s\)]\)", "S"}]

Figure 3

Now we can extend this research for general case when a follows a list, for example

d = 3; zh = 3/2; a0 = {1/20, 1/10, 1/5}; 
a[zs_?NumericQ, c_?NumericQ] := 
 c zs^(d + 1)/(d + 1) NIntegrate[
   Sqrt[1/(1 - (zs/zh)^(d + 1) y)/(1 - 
       c^2 zs^(2 d) y^(2 d/(d + 1)))], {y, 0, 1}];
Table[plot[i] = 
  ContourPlot[(a[zs, c] - a0[[i]]) Boole[c^2 zs^(2 d) <= 1] == 0, {zs,
     0, zh}, {c, 0, 1}, PlotPoints -> 50, 
   FrameLabel -> {"\!\(\*SubscriptBox[\(r\), \(s\)]\)", "c"}, 
   PlotLabel -> a0[[i]]], {i, 3}]

Figure 4

We can extract data directly from plot[i] as follows

point = Table[
   plot[i] // 
    Cases[#, GraphicsComplex[points_, ___] :> points, Infinity] &, {i,
     Length[a0]}];

Using point we can define c=f[i][rs] and S[rs,i], and plot it as

Table[f[i] = Interpolation[point[[i]] // First], {i, Length[a0]}];

x0[i_] := Last[point[[i]] // First] // First

S[zs_?NumericQ, i_] := 
 NIntegrate[
  x^-d/Sqrt[(1 - (x/zh)^(d + 1)) (1 - f[i][zs]^2 x^(2 d))], {x, x0[i],
    zs}]

Table[Plot[S[zs, i], {zs, x0[i], zh}, 
  AxesLabel -> {"\!\(\*SubscriptBox[\(z\), \(s\)]\)", "S"}, 
  PlotLabel -> a0[[i]]], {i, Length[a0]}]

Figure 5

Note, this solution is available due to constraints follow from the suggestion that a, rh,rs,c are real only.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.