Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to integrate the following function:

$$ r_e=\int_0^z\frac{dz}{\sqrt{\Omega_m(1+z)^3+\Omega_\Lambda}} $$

where $\Omega_\Lambda=1-\Omega_M$ and both $\Omega_\Lambda$ and $\Omega_M$ are real and positive.

Doing the definite integral using Integrate doesn't produce a result even after a few hours of running. I found that just doing the indefinite integral and then putting in the limits later produces a result. Why is this?

Additionally, the result is unexpected and discontinuous for different values of $z$ depending on the value of $\Omega_M$. Any suggestions on the best strategy for integrating this?

FunH[omM_,omL_,z_]:=1/((omM(1+z)^3+omL)^(1./2.))
F[omM_,z_]:=Evaluate[Block[{omL=1-omM},Integrate[FunH[omM,omL,z],z]]]
Table[(1+z)^(-1)*(F[omM,z]-F[omM,0]),{omM,0.000001,1,0.1},{z,0.000001,3,0.001}]

I had to copy the code from another machine so I might have missed a bracket somewhere.

EDIT: If anyone is interested, there is a publicly available paper dealing with a transformation of this integral to make the solution easier: http://mnras.oxfordjournals.org/content/412/4/2685.full

Essentially one can transform the above integral using a few simple substitutions into something which is proportionate to the difference between two separate integrals of the form: $$ \int\frac{du}{\sqrt{u^4+u}} $$ This can be easily integrated by Mathematica.

share|improve this question
    
You should be aware that your denominator may have some problems ... –  belisarius Jan 23 '13 at 23:35
add comment

3 Answers

In our case the symbolic definite integral appears to be much more difficult to evaluate. This is because the complex functions underlying the given integral strongly depend on the parameters and may involve branch points and it implies a huge number of possible cases which the system tries to collect in ConditionalExpression. In general one might assume sufficient number of relations between the symbolic parameters, but in our case it seems to fail. As a more convenient way we define this function :

int[ ΩM_, ΩΛ_, x_] := Integrate[ 1/Sqrt[ΩM (1 + z)^3 + ΩΛ], {z, 0, x}]

now we can evaluate definite integrals in various cases :

int[2, 1, 1]
  (1/(3^(1/4)))(-1)^(1/12) 2^(2/3) ( 
     EllipticF[ ArcSin[((-1)^(11/12) Sqrt[1 + 2^(1/3)])/3^(1/4)], (-1)^(1/3)]
     - EllipticF[ ArcSin[((-1)^(11/12) Sqrt[1 + 2 2^(1/3)])/3^(1/4)], (-1)^(1/3)])

Nevertheless built-in rewrite rules are not always capable to give interesting results, e.g. in ver.9 the first and the second integrals are unevaluated unlike in ver.8

{int[1, 2, 3], int[2, Pi/2, 1], int[Sqrt[3], Pi, -1]} // TraditionalForm

enter image description here

Many assumptions underlying the symbolic definite integral appear to be a no-go issue unlike in a case by case basis of the definite integral.

share|improve this answer
    
surprisingly, Integrate[ 1/Sqrt[\[CapitalOmega]M (1 + z^3) + \[CapitalOmega]\[CapitalLambda]], {z, 0, x}, GenerateConditions -> False] takes a long time on 9 as well?? –  chris Jan 24 '13 at 0:21
    
I am currently using Mathematica 8. This function is used in cosmology and I've found a few papers dealing with just this issue. By massaging the function a bit and re-expressing it as the difference of two integrals of the form: Integrate[1/(u^4+u),{u,0,s}] - Integrate[1/(u^4+u),{u,0,s/(1+z)}] where s is a function of omega. Produces satisfactory results. –  Jacobo Blanco Jan 24 '13 at 0:26
    
It takes a bit but not so long (say less than a minute). –  Artes Jan 24 '13 at 0:26
1  
If OmegaM is negative and x at all large then the argument to the square root eventually becomes negative. So there must be other restrictions. –  Daniel Lichtblau Jan 24 '13 at 21:22
1  
I apologize, you're right OmegaM and OmegaLambda are both positive and real –  Jacobo Blanco Jan 25 '13 at 10:09
show 5 more comments

Start with the transformed integral

$$\begin{align*} \int_0^z\frac{\mathrm du}{\sqrt{a(u+1)^3+b}}&=\int_1^{z+1}\frac{\mathrm du}{\sqrt{au^3+b}}\\ &=\frac{\sqrt[3]{b/a}}{\sqrt{b}}\int_{\frac1{\sqrt[3]{b/a}}}^{\frac{z+1}{\sqrt[3]{b/a}}}\frac{\mathrm du}{\sqrt{u^3+1}} \end{align*}$$

Now, I had already dealt with that elliptic integral in this math.SE answer, so I'm not repeating the work here; suffice it to say that we now have an explicit expression for your integral in terms of EllipticF[]:

With[{a = 4/3, b = 5/2, z = 7},
     {(* integral *)
      NIntegrate[1/Sqrt[a (t + 1)^3 + b], {t, 0, z}, WorkingPrecision -> 20],
      (* closed form solution *)
      N[((b/a)^(1/3)/(3^(1/4) Sqrt[b]))
        (EllipticF[ArcCos[2 Sqrt[3]/(1 + Sqrt[3] + (z + 1)/(b/a)^(1/3)) - 1], (2 + Sqrt[3])/4] -
         EllipticF[ArcCos[2 Sqrt[3]/(1 + Sqrt[3] + 1/(b/a)^(1/3)) - 1], (2 + Sqrt[3])/4]), 20]}]

   {0.97715948427474779292, 0.97715948427474779300}

Alternatively, since EllipticF[ArcCos[u], m] == InverseJacobiCN[u, m], you can also use this function:

With[{a = 4/3, b = 5/2, z = 7},
 N[((b/a)^(1/3)/(3^(1/4) Sqrt[b]))
   (InverseJacobiCN[2 Sqrt[3]/(1 + Sqrt[3] + (z + 1)/(b/a)^(1/3)) - 1, (2 + Sqrt[3])/4] - 
    InverseJacobiCN[2 Sqrt[3]/(1 + Sqrt[3] + 1/(b/a)^(1/3)) - 1, (2 + Sqrt[3])/4]), 20]]
share|improve this answer
add comment

For reasons noted in a comment, I don't think oM should be allowed to go negative. Not sure if it can be arbitrarily large. Let's say we are operating under the assumption that 0 < oM.

This is what I would like to evaluate:

ii = Integrate[1/Sqrt[oM (1 + z)^3 + (1 - oM)], {z, 0, x}, 
  Assumptions -> {x > 0, 0 < oM < 1}]

Unfortunately, this does not succeed. The reason is that the integration code becomes very sulky when it sees elliptics, especially with parameters, due to the difficulty of discerning path singularities.

So here is an old trick. Give it an explicit value for oM, one that you can use to later recover a (possibly correct) parametrized result. I use $1/e$ for this purpose (since I don't see $e$ in the indefinite integral).

Timing[
 ii = Integrate[1/Sqrt[1/E (1 + z)^3 + (1 - 1/E)], {z, 0, x}, 
    Assumptions -> 0 < x];]

(* Out[1]= {21.960000, Null} *)

Numeric check:

Chop[Table[ii, {x, 1., 11., 1}]]

(* Out[2]= {0.742608949781, 1.14591119954, 1.39565792161, 1.56802998702, 
1.69583835424, 1.79538009669, 1.87570826106, 1.94228852684, 
1.9986361674, 2.04712646089, 2.08942907726} *)

Table[NIntegrate[1/Sqrt[1/E (1 + z)^3 + (1 - 1/E)], {z, 0, x}], {x, 1., 11., 1}]

(* Out[3]= {0.742608949781, 1.14591119954, 1.39565792161, 1.56802998702, 
1.69583835424, 1.79538009669, 1.87570826106, 1.94228852684, 
1.9986361674, 2.04712646089, 2.08942907726} *)

Okay, so it's fine at least for one value of the parameter oM, and a handful of values for x. (That proves that it never fails. NEVER.) We can now do

ii2 = ii /. E -> 1/oM

Use at your own risk, though. It could become garbage depending on settings for x and/or oM.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.