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


3 Answers 3


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: enter image description here

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 -> ∞}]
  • $\begingroup$ How's your eqn is so different than OP's? $\endgroup$
    – m0nhawk
    Commented May 3, 2020 at 20:23
  • 1
    $\begingroup$ 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. $\endgroup$ Commented May 3, 2020 at 20:25
  • $\begingroup$ Strange, I can reproduce this only for z1+1 as new variable, and Assumptions -> z>0. Otherwise it's OP's function. $\endgroup$
    – m0nhawk
    Commented May 3, 2020 at 20:33
  • $\begingroup$ What version do you have? $\endgroup$ Commented May 3, 2020 at 20:34
  • $\begingroup$ 12.1. Yours? And did you specify Assumptions for eqn[z_]? $\endgroup$
    – m0nhawk
    Commented May 3, 2020 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<D_c<D_{m}=3.25664$ 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:

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

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[
                     {p, {Re, Im}}]

enter image description here

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

and for $0< z<D_m$ e.g.

z[ 2, OM, OL]//N//Chop
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]

enter image description here

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}]]

enter image description here

Consequently, eqn can never equal the specified value of Dc

Using instead

Dc = 2;

Use NSolve

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

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

Or FindRoot

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]

enter image description here

zvalue // N

(* z == 7.13731 *)

Similarly with Solve

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

enter image description here


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.