Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica : Numerical Integration

  1. Aug 13, 2012 #1

    Hepth

    User Avatar
    Gold Member

    So I have to numerically integrate a function that is quite difficult. I could post it, but its long.
    Regardless, there are 5 variables I am integrating over. The first is independent. The next 4 are subsequently dependent on the integration limits of the previous. Not a big deal.

    something like

    Integrate[F[x,y,z,t],{x,XCUT, 1},{y,0,Sqrt[1-x^2]},{z,0,1-x-y},{t,-1,1}]

    but more.

    Now, my problem comes in that XCUT is very close to 1. So say 0.977. That means that y only varies from 0 to something small, and z is also small. So its a very small slice of phase space that I'm integrating over.

    What happens is my basic numerical integration outputs something like :

    NIntegrate obtained 3.02466*10^-6 and 3.08509*10^-6 for the integral and error estimates.

    So the error is appx the integral.

    So what else can I do? Monte carlo seems to give me an answer without error, but seems to be heavily dependent on the maximum number of points I allow it to use.

    Is there a way to do a HIGHLY ACCURATE numerical integration, where some of the limits might be close to poles. Time is not an issue, I can throw it on the cluster.

    Pumping up WorkingPrecision doesn't seem to help much. Same with max recursion.
     
  2. jcsd
  3. Aug 13, 2012 #2
    What happens with

    Integrate[F[x,y,z,t],{x,977/1000, 1},{y,0,Sqrt[1-x^2]},{z,0,1-x-y},{t,-1,1}]

    and there are absolutely no decimal points anywhere in the definition of F[]?
     
  4. Aug 13, 2012 #3

    Dale

    Staff: Mentor

    Beyond specifying everything that you see as exact numbers, one other thing you can try is to make a transformation [itex]x= 1 + \chi[/itex] everywhere.
     
  5. Aug 13, 2012 #4

    Hepth

    User Avatar
    Gold Member

    I made sure not to have ANY decimals. The integral itself is already scaled with only one input parameter that I put in as "2/5".

    Same thing.

    I'll attach a notebook with one of my integrals.


    The thing is there are no numbers really. Its basically a small section of phase space that I'm integrating over, and I wonder if its just too small, or what to change to get a consistent answer.
     

    Attached Files:

    Last edited: Aug 13, 2012
  6. Aug 13, 2012 #5

    Hepth

    User Avatar
    Gold Member

    I'm starting to think the integral is divergent, but then monte carlo integration shouldn't be getting smaller the larger the number of points I give it...

    Maybe its actually Zero, but in some roundabout, definite-integral way.

    EDIT: few minutes, I'll get the actual integrand.
    EDIT: the above notebook has been updated
     
    Last edited: Aug 13, 2012
  7. Aug 15, 2012 #6
    In your current notebook your integrand has a large repeated subexpression.

    If you look at this, and my scrape-n-paste hasn't failed,

    In[1]:= FN /. (2 S2^2+T1 (S1+T1) (-1+S1+T1+U2) - S2 (S1 (2+T1-U2) + 2 (-1+U2) + T1 (2+T1+U2)) - 2 Sqrt[S2 (-1+S1-S2+T1) (-S2+S1*T1) (S2 + S1 (-1+U2) + (-1+U2) (-1+T1+U2))] \[Zeta]2)/(-4 S2 + (S1+T1)^2) -> V1

    Out[1]= -(((S2*(-1+T1) - T1*(-1+S1+T1))*(-(S1*(-1+T1+U2-V1)) + (-U2+V1)*(1-T1-U2+V1)))/ ((29/25-T1-U2+V1)^2*Sqrt[(-S1^2 - (-1+S2)^2 + 2*S1*(1+S2))*(-1+\[Zeta]2^2)]))

    The numerator is a quadratic in that big subexpression V1.

    I'm thinking you might be able to get the numerator into a form a+(b+V1)^2 for some a and b and this might make it substantially more compact and easier to understand, and perhaps even see whether the integral is zero or not. I'm even questioning whether many of the S, T and U terms might get absorbed in that more compact result. Unfortunately I have not been able to coax Mathematica into doing that for me yet.

    Here is one attempt at the numerator

    -((S2*(-1+T1) - T1*(-1+S1+T1))^2*(1+S1-T1-2*U2)^2)/4 + (S2*(-1+T1) - T1*(-1+S1+T1))*(S1-U2)*(-1+T1+U2)
    +(((S2*(-1+T1) - T1*(-1+S1+T1))*(-1-S1+T1+2*U2))/2 + V1)^2

    but the S,T and U outside the squared term didn't mostly disappear inside the squared term as I had hoped they might. Perhaps I've made a mistake or perhaps there is another way to reorganize this to compact it further.

    Sometimes complicated integrals can be turned into a sum of substantially smaller integrals, sometimes using Apart to split apart the numerator. There have been examples where this resulted in a vast number of simpler integrals, most of which were able to be integrated exactly, not using NIntegrate, but I'm not sure whether your denominator would allow this much enhancement.
     
    Last edited: Aug 15, 2012
  8. Aug 15, 2012 #7

    Hepth

    User Avatar
    Gold Member

    Actually there was a "T2" that got replaced in terms of the {S1,S2,T1,U2,zeta}, as it and zeta are interchangeable. I think I can actually integrate out the zeta and just get a factor of "pi" overall, for some of the terms (The ones odd in "T2").

    I'll try and clean it up some. So does mathematica's numerical integration care how messy or clean the integral is? I thought it wouldnt matter since its getting evaluated at certain points, and the algebraic form doesnt matter.
     
  9. Aug 15, 2012 #8
    When you are doing more additions and subtractions and divisions of low precision numbers every one of those operations is contributing to the overall error of the result.

    Since you were subtracting low precision values that were close to each other I was actually thinking you were using Integrate, as you explicitly showed in your very first post in this thread, and sidestep NIntegrate, which deals in nothing but approximate values even after you substituted exact rational numbers as I suggested. That NIntegrate instantly threw away all the exact parameters I recommended and began fiddling with floating point and accumulated error.

    If you know more about the problem than the contents of your notebook then I would hope the problem might be simplified to the point where you could get good accurate answers quickly. Try to get rid of the redundancy in the subexpressions. That might help you and Mathematica see that the problem can be simplified even further.
     
  10. Aug 15, 2012 #9
    If you go to the help entry for Integrate (and NIntegrate...), the meaning of a multiple
    integral:

    NIntegrate[f[x,y],{x,xmin,xmax},{y,ymin,ymax}]
    is
    [tex]
    \int_{x_\mathrm{min}}^{x_\mathrm{max}}{dx \, \int_{y_\mathrm{min}}^{y_\mathrm{max}}{dy \, f(x,y)}}
    [/tex]
    i.e. the order of the integrations is reversed! Make sure you write the order of the integrands in the Mathematica command the way you want it performed.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook