# Solving equations involving integrals

I need to find the value of $$z$$ for a particular value of $$D_c$$ (eg. $$500$$), but $$z$$ is inside an integral, and I'm not able to use Solve since the integral is giving Hypergeometric2F1 function as the output.

OmegaM = 0.3111;
OmegaLambda = 0.6889;
Dc = 500;

eqn = Integrate[(OmegaM (1 + z1)^3 + OmegaLambda)^(-1/2), {z1, 0, z},
Assumptions -> z > 0]

-1.1473+(1.20482+1.20482z)Hypergeometric2F1[0.333333,0.5,1.33333,-0.451589(1.+z)^3]

zvalue = Solve[eqn == Dc, z]

Solve was unable to solve the system with inexact coefficients or the
system obtained by direct rationalization of inexact numbers present
in the system. Since many of the methods used by Solve require exact
input, providing Solve with an exact version of the system may help.


Is there any other way I can solve this equation?

Also, Integrate is taking some time and I'd like it to be fast since I need to put it in a loop with lots of $$z$$ values to be computed for corresponding $$D_c$$ values.

• There is a factor of c/H0 multiplied with the integral where c=300000, H0=68. I think that will make the equation consistent. Thanks mathematica.stackexchange.com/users/35368/superciocia for suggesting FindRoot, it solves the issue. Now I just want to know if I can speed up this integral somehow. – Divyajyoti May 4 at 5:45

From your integral I get the following equation (same as yours when you plot it):

eqn[z_] :=
3.2566440560469836 - (
3.5857498598223954 Hypergeometric2F1[1/6, 1/2, 7/
6, -(2.2144005143040824/(1 + z)^3)])/Sqrt[1 + z]


Both Solveand NSolve fail.

So I tried FindRoot:

Dc = 3.1;
FindRoot[ eqn[z] - Dc, {z, 0}]

{z -> 523.001}


which agrees with a graphical solution:

I don't think it has a solution for Dc=500 as the eqn flattens out to 3.25664 as $$z\rightarrow \infty$$:

Limit[ eqn[z], {z -> ∞}]

3.25664

• How's your eqn is so different than OP's? – m0nhawk May 3 at 20:23
• It's the same function. I had to change the integrand from $z1+1$ to $z2$ in order for Mathematica to give me an analytical result. Maybe it's because of different versions. – SuperCiocia May 3 at 20:25
• Strange, I can reproduce this only for z1+1 as new variable, and Assumptions -> z>0. Otherwise it's OP's function. – m0nhawk May 3 at 20:33
• What version do you have? – SuperCiocia May 3 at 20:34
• 12.1. Yours? And did you specify Assumptions for eqn[z_]? – m0nhawk May 3 at 20:35

We will demonstrate that the exact formula for $$z$$ reads: $$z=\wp\bigg(\frac{\sqrt{\Omega_M}}{2}D_c+\wp^{-1}\big(1;0,-\frac{4\Omega_\Lambda}{\Omega_M}\big);0,-\frac{4\Omega_\Lambda}{\Omega_M}\bigg)-1$$ where $$\wp(x;g_2,g_3)$$ is the Weierstrass elliptic function, which yields a value $$w$$ in the elliptic integral $$x=\int^{w}_{\infty}\frac{d t}{\sqrt{4t^3-g_2\;t-g_3}}$$ and so generalizing the answer to the original question for any integrand of the form $$\frac{1}{\sqrt{R(t)}}$$, where $$R(t)$$ is a fourth or a third order polynomial in $$t$$. This formula can be implemented in the following way:

z[ Dc_, OM_, OL_]:= WeierstrassP[ Sqrt[OM/4] Dc+ InverseWeierstrassP[ 1, { 0,-4OL/OM}],
{ 0, -4OL/OM}]-1


We rationalize numeric constants to play with the system seamlessly (although this step is not neccesary):

{ OM, OL} = Rationalize[{ OmegaM = 0.3111, OmegaLambda = 0.6889}];


Let's derive $$z$$: $$D_c=\int^{z}_{0}\frac{d s}{\sqrt{\Omega_M (s+1)^3+\Omega_{\Lambda}}}=\frac{2}{\sqrt{\Omega_M}}\int^{z+1}_{1}\frac{d s}{\sqrt{4 s^3+\frac{4\Omega_{\Lambda}}{\Omega_M}}}=\\=\frac{2}{\sqrt{\Omega_M}}\Bigg(\int^{\infty}_{1}\frac{d s}{\sqrt{4 s^3+\frac{4\Omega_{\Lambda}}{\Omega_M}}}-\int^{\infty}_{z+1}\frac{d s}{\sqrt{4 s^3+\frac{4\Omega_{\Lambda}}{\Omega_M}}}\Bigg)=\\=\frac{2}{\sqrt{\Omega_M}}\Bigg(-\wp^{-1}\big(1;0,-\frac{4\Omega_\Lambda}{\Omega_M}\big)+\wp^{-1}\big(z+1;0,-\frac{4\Omega_\Lambda}{\Omega_M}\big)\Bigg)$$ and this implies our formula for $$z$$.

The formula for $$z$$ is valid in the range $$0 and we can also derive an exact formula for $$D_m$$: $$D_m=\frac{2}{\sqrt{\Omega_M}} \Re\Big( 2\;\omega_{1}(0,g_3)-\wp^{-1}\big(1;0,g_3\big)\Big)$$ where $$\Re$$ is the real part, $$\omega_{1}(0,g_3)$$ is the Weierstrass half period and $$g_3$$ is the Weierstrass invariant, in our case $$g_3=-\frac{4\Omega_{\Lambda}}{\Omega_M}$$, implementing it:

g3=-4OL/OM;
Dm = 2/Sqrt[OM]( 2WeierstrassHalfPeriodW1[{0, g3}]-InverseWeierstrassP[1,{0, g3}])//Re//N

3.25664


Dm been calculated in version 12.1, however in earlier versions one should evaluate simply Dm = -2/Sqrt[OM] InverseWeierstrassP[1,{0, g3}]. This is because InverseWeierstrassP[1,{0, g3}] is computed in a adjacent parallelogram (see e.g this discussion). One should also note better handling of symbolic input in WeierstrssHalfPeriodW1 etc. For the presentation of the structure of $$z$$ being an elliptic function (shifted and rescaled $$\wp$$) we define:

wHP = Through @ { WeierstrassHalfPeriodW1,WeierstrassHalfPeriodW2,
WeierstrassHalfPeriodW3} @{ 0,-4OL/OM}//ReIm
// FullSimplify;
GraphicsRow @ Table[ ContourPlot[ Evaluate @ Table[p[z[x+I y,OM,OL]] ==k,
{k, wHP[[#1,#2]]& @@@ {{2,1},{2,2},{3,1},{3,2}}}],
{x, -8, 8}, {y, -8, 8}, ContourStyle ->Thread[
{Thick,{Red,Darker@Cyan,Darker@Green,Orange}}]],
{p, {Re, Im}}]


There was an assumption that $$z>0$$, however $$D_c=500$$ can be reached for negative $$z$$, e.g.

z[ 500,OM, OL]//N//Chop

-1.73134


and for $$0< z e.g.

z[ 2, OM, OL]//N//Chop

7.13731

OmegaM = 0.3111 // Rationalize;
OmegaLambda = 0.6889 // Rationalize;
Dc = 500;

eqn = Integrate[(OmegaM (1 + z1)^3 + OmegaLambda)^(-1/2), {z1, 0, z},
Assumptions -> z > 0]


For z > 0, eqn is monotonically increasing

Assuming[z > 0, D[eqn, z] > 0 // Simplify]

(* True *)


The maximum value of eqn is

(lim = Limit[eqn, z -> Infinity]) // N

(* 3.25664 *)

LogLinearPlot[{lim, eqn}, {z, 10^-2, 10^4},
PlotLegends -> Placed["Expressions", {.3, .7}]]


Consequently, eqn can never equal the specified value of Dc

Dc = 2;

zvalue = NSolve[{eqn == Dc, z > 0}, z]

(* {{z -> 7.13731}} *)

zvalue = FindRoot[eqn == Dc, {z, 1}]

{z -> 7.13731}


Or Reduce (provides the exact value as a Root expression)

zvalue = Reduce[{eqn == Dc, z > 0}, z]


zvalue // N

(* z == 7.13731 *)


Similarly with Solve

zvalue = Solve[{eqn == Dc, z > 0}, z][[1]]
`