# Solving a function inside an integral

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;
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.

• Regarding your last block: there errors you have got are due to MMA not substituting the value of z (because of b for instance). – anderstood Sep 19 at 16:36
• @anderstood So the fix there is to add the function With? – mathemania Sep 19 at 17:20
• Not only, I also introduced the dependencies explicitly (this way I can test that each function is evaluated as I expect), used SetDelayed instead of =. – anderstood Sep 19 at 17:23

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;
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.

• This seems to be pertaining to the first expression for $a$ which has a closed form but please look at my UPDATE post. I have tried to use your code corresponding to the update in my post, but it seems to be not working since it shows "The integrand has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}." – mathemania Sep 19 at 16:12
• @mathemania If you read and understand the answer, you'll see that I use nowhere the closed-form expression (btw I use NIntegrate). – anderstood Sep 19 at 16:21
• Yes, but when you substitute it in $S$ and do a plot it brings out an error message as in my first comment. Does this have anything to do with $c$ or the expression of $S$ itself? – mathemania Sep 19 at 16:34
• @mathemania See edit – anderstood Sep 19 at 16:41
• I see, but what if the range of evaluation for $z_s$ is $[0,10]$, is there a way to chop this to avoid those complex values? – mathemania Sep 19 at 17:17
ListLinePlot@
NDSolveValue[{x'[z] == 1, x == 1, toroot[c[z], z] == 0,
c == cz}, c, {z, 0.4, 5}] • This seems to be pertaining to the first expression for $a$ which has a closed form so it is not really an issue on how to evaluate it using the usual NIntegrate and FindRoot but please look at my UPDATE post. – mathemania Sep 19 at 16:09
• @mathemania When I looked at the Q (the first time), there was only one bit of code for the integral. Are you saying that it not the integral you wanted help with? – Michael E2 Sep 19 at 16:21
• Yes there was only one bit of code which corresponds to the first $a$ integral which has a closed form answer, but my question is how to evaluate an integral in general even in the case where there is no closed form answer which I posted in my UPDATE. – mathemania Sep 19 at 16:28
• I think there is already nothing wrong with $cz[z]$, it is when you plug it in $S$ where problem occurs. – mathemania Sep 19 at 17:18
• Interesting combination of ListLinePlot with InterpolatingFunction. – anderstood Sep 19 at 17:28

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.

• Nice, is there a way to do a FindRoot as a list? That is, in your code containing cc[d_?NumericQ, a_?NumericQ, z_?NumericQ] is it possible to input a range of z so that it outputs a range of cc? I'm thinking of Table but seems like it is not working. – mathemania Sep 16 at 9:50
• Try Table[{cc[3, 0.2, z]}, {z, 1, 10, 1/2}] or Plot[cc[3, 0.2, z], {z, 1, 10}]. – Roman Sep 16 at 13:35
• Your suggestion worked, but how to do this if the integral does not have a closed form just like the hypergeometric function? Please see the update on my question. – mathemania Sep 16 at 14:22
• Integrate[x^d/Sqrt[1 - c^2*x^(2*d)], x]==(x^(1 + d)*Hypergeometric2F1[1/2, (1 + d)/(2*d), (1/2)*(3 + 1/d), c^2*x^(2*d)])/(1 + d). Seems important. – Steffen Jaeschke Sep 20 at 9:19