Solving a function inside an integral

Question

I have an integral without a closed form answer and 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}$$

where $a$, $z_s$, and $z_h$ are constants which can be assigned a value ($0 \leq a, z_s \leq 10$, $z_h=10$), $c=c(z_s)$ is an unknown function of $z_s$ which I need to determine from $\eqref{1}$.

After determining $c(z_s)$ in terms of $z_s$ I will plug that into the expression $S$ given by,

\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{2d-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}$.

I have tried the following code

d = 3;
zh = 10;
SeedRandom[2020];
a = RandomReal[{0, 10}];
toroot[c_?NumericQ, z_] := a - c*z^(d + 1)*NIntegrate[x^d*((1 - (z/zh)^(d + 1) x^(d + 1)) (1 - c^2*z^(2 d) x^(2 d)))^(-1/2), {x, 0, 1}]
cz[z_?NumericQ] := c /. FindRoot[toroot[c, z], {c, 0.002, 0.0000001, 10}]
ints[x_?NumericQ, z_] := With[{b = z/zh}, (((2 d - 1)/(d - 1)) cz[z]^2 z^(2 d)) x^d ((1 - (b x)^(d + 1))/(1 - cz[z]^2 (z x)^(2 d)))^(1/2) - ((b^(d + 1) (d + 1))/(2 (d - 1))) x ((1 - cz[z]^2 (z x)^(2 d))/(1 - (b x)^(d + 1)))^(1/2) + (b^(d + 1) x)/((1 - (b x)^(d + 1)) (1 - cz[z]^2 (z x)^(2 d)))^(1/2)];
intS[z_?NumericQ] := NIntegrate[ints[x, z], {x, 0.5, 1}]
functionS[z_?NumericQ] := ((-((1 - cz[z]^2 z^(2 d)) (1 - b^(d + 1)))^(1/2)/(d - 1)) + intS[z] + 1)/(4 z^(d - 1));
function[z_?NumericQ] := Log[10, functionS[z]];
Plot[function[z], {z, 0, 10}, PlotPoints -> 3, AxesLabel -> {"z", "Log S"}, PlotStyle -> {Thick, Blue}, PlotRange -> Full, ImageSize -> Large]

How should I obtain an expression or at least a numerical fit of $c(z_s)$ in terms of $z_s$ given that $\eqref{1}$ does not have a closed form? Maybe Mathematica has a way to give a relation even it is not exact, like determining the integral using approximate methods? In the end, I will plug $c(z_s)$ in $S$ and plot.

***Note: This is an updated version of my original question, I have removed the redundancies in my original post and clarified my problem.

Answer 1

Your problem is: given $a,d,z$, find a $c$ that corresponds to a root of $f(c)=a-cz_s\int \dots$. It happens that here, you want $f$ to involve numerical integration.

d = 3;
SeedRandom[2020];
a = RandomReal[{0, 10}];
toroot[c_?NumericQ, z_] := 
 a - c*z^(d + 1)*NIntegrate[x^d*(1 - c^2*z^(2 d) x^(2 d))^(-1/2), {x, 0, 1}]
cz[z_?NumericQ] := c /. FindRoot[toroot[c, z], {c, 0.002, 0.0000001, 10}]

Plot[cz[z], {z, 0.4, 5}, PlotPoints -> 3, AxesLabel -> {"z", "c[z]"}]

And that the same plot, in red, on your 3D plot:

Note: Depending on the values of the parameters $a,d,z$, you might have issues in the numerical integration.

Note 2: If, when changing the integrands or playing with the parameters, the above approach is not robust enough, you might want to use numerical continuation to retrieve $c(z)$ from one pair $(c(z_0), z_0)$.

---

Another approach based on ContourPlot: extract the slice from your 3D plot with ContourPlot, extract the points from the plot, interpolate. It is not as accurate (see Is Mathematica `ContourPlot` function really so efficient?) but it is easy to understand.

contour = ContourPlot[inta[zs, c] == a, {zs, 0, 2}, {c, 0, 10}, 
   PlotPoints -> 10] // Quiet
data = contour[[1, 1, 1]];
cinter = Interpolation[data, InterpolationOrder -> 1];
   Plot[{cinter[x], cz[x]}, {x, 0.4, 2}, PlotLabels -> {"FindRoot", "ContourPlot"}]

---

Regarding your edit: your syntax does not allow MMA to know the value of z when it integrates, hence the error message.

The following works:

ints[x_?NumericQ, z_] := 
  With[{b = 
     z/zh}, ((((2 d - 1)/(d - 1)) cz[z]^2 z^(2 d)) x^
         d (1 - (b x)^(d + 1))^(1/2))/(1 - 
        cz[z]^2 z^(2 d) x^(2 d)) (((b^(d + 1) (d + 1))/(2 (d - 
              1))) x (1 - cz[z]^2 z^(2 d) x^(2 d))^(1/
           2))/(1 - (b x)^(d + 1))^(1/
         2) + (b^(d + 1) x)/((1 - (b x)^(d + 1))^(1/2) (1 - 
          cz[z]^2 z^(2 d) x^(2 d))^(1/2))];
intS[z_?NumericQ] := NIntegrate[ints[x, z], {x, 0.5, 1}]
functionS = (-((1 - cz[z]^2 z^(2 d)) (1 - b^(d + 1)))^(1/2)/(d - 1) - 
     intS[z] + 1/z^(d - 1))/4;
function = Log[10, functionS];

function /. z -> 0.5
(* -0.0560923 *)

function /. z-> 2.
(* -1.20545 + 1.36438 I *)

Plot[function, {z, 0.4, 1.3}, PlotPoints -> 3, AxesLabel -> {"z", "Log S"}]

Note however that function is complex-valued on some domains, so you cannot always plot it directly.

Answer 2

ListLinePlot@
 NDSolveValue[{x'[z] == 1, x[1] == 1, toroot[c[z], z] == 0, 
   c[1] == cz[1]}, c, {z, 0.4, 5}]

Answer 3

First the obvious part: performing the integral gives the equation

$$ a = \frac{c z^{d+1} \, _2F_1\left(\frac{1}{2},\frac{d+1}{2 d};\frac{3 d+1}{2 d};c^2 z^{2d}\right)}{d+1} $$

Substitute $x=c z^d$: $$ \frac{a}{z} = \frac{x \, _2F_1\left(\frac{1}{2},\frac{d+1}{2 d};\frac{3 d+1}{2 d};x^2\right)}{d+1} $$

For any given value of $\frac{a}{z}$ satisfying $\left|\frac{a}{z}\right|\le \frac{\sqrt{\pi } \Gamma \left(\frac{3 d+1}{2 d}\right)}{(d+1) \Gamma \left(\frac{2d+1}{2 d}\right)}$ you can find the corresponding value of $x$ (and therefore of $c=x/z^d$) by numerical solving,

With[{az = 0.2, d = 3},
  FindRoot[az == (x*Hypergeometric2F1[1/2,(d+1)/(2d),(3d+1)/(2d),x^2])/(d+1),
           {x, az*(d+1)}]]
(*    {x -> 0.705645}    *)

where the starting point of the root search, $x_0\approx\frac{a}{z}(d+1)$, is found from the linear approximation of the hypergeometric term for small $x$.

All packaged up in one function,

cc[d_?NumericQ, a_?NumericQ, z_?NumericQ] :=
  x/z^d /. FindRoot[a/z == (x*Hypergeometric2F1[1/2,(d+1)/(2d),(3d+1)/(2d),x^2])/(d+1),
                    {x, a/z*(d+1)}]

cc[3, 0.2, 1]
(*    0.705645    *)

Alternatively, assuming that $\left|\frac{a}{z}\right|$ is very small, you could use a series-approximation of the hypergeometric series,

Series[(x*Hypergeometric2F1[1/2,(d+1)/(2d),(3d+1)/(2d),x^2])/(d+1), {x, 0, 4}]
(*    x/(d+1) + x^3/(6d+2) + O(x^5)    *)

Use this for cubic (or even higher-order) root-finding:

f3[az_, d_] = Root[Function[x, x/(1 + d) + x^3/(2 (1 + 3 d)) - az], 1];
f3[0.2, 3]
(*    0.724076    *)

f5[az_, d_] = Root[Function[x, x/(1 + d) + x^3/(2 (1 + 3 d)) + (3 x^5)/(8 (1 + 5 d)) - az], 1];
f5[0.2, 3]
(*    0.711057    *)

In all these Root objects I've picked the first root because it is the only real-valued one.