Skip to main content
Mathematica

Return to Answer

added 1533 characters in body
Source Link
Alex Trounev
Alex Trounev
  • 48.8k
  • 3
  • 51
  • 115

If all parameters $a,z_s,z_h,c$ are real, then we have constrainsconstraints 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})$:

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.

If all parameters $a,z_s,z_h,c$ are real, then we have constrains 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})$:

Figure 3

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})$:

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.

Source Link
Alex Trounev
Alex Trounev
  • 48.8k
  • 3
  • 51
  • 115

If all parameters $a,z_s,z_h,c$ are real, then we have constrains 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