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.

Here's m problem simpler in terms of codes

Clear["Global`*"]
e = 0.454648; A = 1.6297; rs = 3.68888; kb = 0.001872041; Na = 
 6.02214129*10^23;
ep = 0.2961; p = 3.78;
pot[r_] := -e (1 - (1 - Exp[-A (r - rs)])^2)

B[T_] := -2 \[Pi] Na NIntegrate[(Exp[-pot[r]/(kb T)] - 1) r^2, {r, 0, 
    Infinity}] 10^-24
C2[T_] := 8 \[Pi]^2 /3 Na ^2 10^-48 NIntegrate[
   NIntegrate[
    NIntegrate[
     r12^2 (1 - Exp[-pot[r12]/(kb T)]) r13^2 (1 - 
        Exp[-pot[r13]/(kb T)]) (1 - 
        Exp[-pot[Sqrt[r12^2 + r13^2 - 2 r12 r13 n]]/(kb T)]), {r13, 0,
       150}], {r12, 0, 150}], {n, 1, 150}]

My problem is that second integration function C2[T_] keeps giving me zero with 3 warnings of

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>

and if I use one NIntegrate function

C2[T_] := 
 8 \[Pi]^2 /3 Na ^2 NIntegrate[
   r12^2 (1 - Exp[-pot[r12]/(kb T)]) r13^2 (1 - 
      Exp[-pot[r13]/(kb T)]) (1 - 
      Exp[-pot[Sqrt[r12^2 + r13^2 - 2 r12 r13 n]]/(kb T)]), {r13, 0, 
    150}, {r12, 0, 150}, {n, 1, 150}]

I get 3 of this warning

NIntegrate::errprec: Catastrophic loss of precision in the global error estimate due to insufficient WorkingPrecision or divergent integral. >>

any ideas how to fix it?

The original integration is

enter image description here

share|improve this question
    
Some suggestions: 1) please don't use capital C because it is already used in Mathematica. Maybe try small c. 2) T is not defined in your equation, neither as a parameter nor as a variable for integration. –  Eric Brown Jul 1 '13 at 5:19
    
sorry I had aleady noticed the capital c error but I have defned T all exponentials are divided by (kb T). –  Raymond Ghaffarian Shirazi Jul 1 '13 at 5:31
    
Another issue: 1) I would combine all the NIntegrates into one 2) Is it possible for the argument to the Sqrt to become negative? Plot3D[r12^2 + r13^2 - 2 r12 r13 n /. n -> 10, {r13, 0, 150}, {r12, 0, 150}] –  Eric Brown Jul 1 '13 at 5:37
    
Another suggestion: if you need the 10^-48 stuff for scaling, then you may want to move that inside of the NIntegrate so that it does not blow up or even vanish. –  Eric Brown Jul 1 '13 at 5:43
    
You are probably calculating a partition function there, the integral is dominated (for low T) by the minimum value of the exponent integrand (ground state energy), therefore it is better to first find the minimum value (and location if you can) to factor out the exponent of the minimum and concentrate the integral only around the location of the minimum (see Laplace and Stationary Phase approximation). –  alfC Jul 1 '13 at 8:30

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.