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 having a hard time to find a solution to a numerical limitation on mathematica:

PwD = 0.005*Integrate[(Exp[-(0.25/\[Tau])]/\[Tau])*
  Integrate[Exp[-((Tan[45*Degree]^2 - (z + 30)^2)/(4*\[Tau]))]*
         (1 + 
      2*Sum[Exp[-((n^2*Pi^2*\[Tau])/10000)]*Cos[0.8*n*Pi]*
         Cos[(n*Pi)/2 - (n*Pi*z)/100], {n, 1, 4}]), {z, -50, 
    50}], {\[Tau], 0, t}]

Mathematica graphics

After solving the integral in space {z,-50,50} I need to evaluate the indefinite integral in time {tau,0,t}, so I eventually generate a plot of PwD as function of Log[t]. But it seems that it does not converge.

Can someone provide an alternative way?


Hi everybody thanks for the comments... They sound about right! In the past 2 months I dedicated my time to derive these solutions again - once I had no help from the author. Well, I finally get it! The published expression has a type mistake... The right equation should be:

PwD = (1/200)*Integrate[
    (1/(E^(0.25/\[Tau])*\[Tau]))*Integrate[
      (1 + 2*Sum[(Cos[0.8*n*Pi]*Cos[(n*Pi)/2 - 
              (n*Pi*z)/100])/E^((n^2*Pi^2*\[Tau])/
             10000), {n, 1, 100}])/
       E^((Tan[\[Psi]*Degree]^2*(z + 30)^2)/(4*\[Tau])), 
      {z, -50, 50}], {\[Tau], 0, t}]

Now I'm able to generate the graphics. For different angles \[Psi] this is what I did:

Pw1 = (Exp[-0.25/t]/t)*
   Integrate[
    Exp[-(((Tan[0 Degree]^2)*(z + 30)^2)/(4*t))]*(1 + 
       2*Sum[Exp[-((n^2)*(Pi^2)*t)/(100^2)]*Cos[0.8*n*Pi]*
          Cos[(n*Pi/2) - (n*Pi*z/100)], {n, 1, 100}]), {z, -50, 50}];
Pw2 = (Exp[-0.25/t]/t)*
   Integrate[
    Exp[-(((Tan[30 Degree]^2)*(z + 30)^2)/(4*t))]*(1 + 
       2*Sum[Exp[-((n^2)*(Pi^2)*t)/(100^2)]*Cos[0.8*n*Pi]*
          Cos[(n*Pi/2) - (n*Pi*z/100)], {n, 1, 100}]), {z, -50, 50}];
Pw3 = (Exp[-0.25/t]/t)*
   Integrate[
    Exp[-(((Tan[60 Degree]^2)*(z + 30)^2)/(4*t))]*(1 + 
       2*Sum[Exp[-((n^2)*(Pi^2)*t)/(100^2)]*Cos[0.8*n*Pi]*
          Cos[(n*Pi/2) - (n*Pi*z/100)], {n, 1, 100}]), {z, -50, 50}];

PD1[y_] := 0.005*NIntegrate[Pw1, {t, 0, y}, MaxRecursion -> 20];
PD2[y_] := 0.005*NIntegrate[Pw2, {t, 0, y}, MaxRecursion -> 20];
PD3[y_] := 0.005*NIntegrate[Pw3, {t, 0, y}, MaxRecursion -> 20];

T1 = Table[{y, 
    PD1[y]}, {y, {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 
     2.2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
     18, 19, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 
     90, 100, 110, 120, 130, 140, 150, 200, 250, 300, 350, 400, 450, 
     500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1100, 
     1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200,
      2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3500, 4000, 
     4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500,
      10000, 11000, 13000, 15000, 17000, 19000, 20000, 22000, 24000, 
     26000, 28000, 30000, 35000, 40000, 45000, 50000, 55000, 60000, 
     65000, 70000, 75000, 80000, 90000, 95000, 100000, 150000, 200000,
      250000, 300000, 350000, 400000, 450000, 500000, 550000, 600000, 
     700000, 800000, 900000, 1000000}}];
T2 = Table[{y, 
    PD2[y]}, {y, {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 
     2.2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
     18, 19, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 
     90, 100, 110, 120, 130, 140, 150, 200, 250, 300, 350, 400, 450, 
     500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1100, 
     1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200,
      2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3500, 4000, 
     4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500,
      10000, 11000, 13000, 15000, 17000, 19000, 20000, 22000, 24000, 
     26000, 28000, 30000, 35000, 40000, 45000, 50000, 55000, 60000, 
     65000, 70000, 75000, 80000, 90000, 95000, 100000, 150000, 200000,
      250000, 300000, 350000, 400000, 450000, 500000, 550000, 600000, 
     700000, 800000, 900000, 1000000}}];
T3 = Table[{y, 
    PD3[y]}, {y, {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 
     2.2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
     18, 19, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 
     90, 100, 110, 120, 130, 140, 150, 200, 250, 300, 350, 400, 450, 
     500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1100, 
     1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200,
      2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3500, 4000, 
     4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500,
      10000, 11000, 13000, 15000, 17000, 19000, 20000, 22000, 24000, 
     26000, 28000, 30000, 35000, 40000, 45000, 50000, 55000, 60000, 
     65000, 70000, 75000, 80000, 90000, 95000, 100000, 150000, 200000,
      250000, 300000, 350000, 400000, 450000, 500000, 550000, 600000, 
     700000, 800000, 900000, 1000000}}];

PwD1 = Interpolation[T1];
PwD2 = Interpolation[T2];
PwD3 = Interpolation[T3];

P1 = LogLogPlot[{PwD1[y], y*PwD1'[y]}, {y, 0.1, 1000000}, 
   PlotRange -> {0.01, 10}, PlotStyle -> {{Black}, {Dashed, Black}}, 
   Frame -> True, FrameLabel -> {tD, "PD e tD*PD'"}, 
   BaseStyle -> {FontSize -> 12}];
P2 = LogLogPlot[{PwD2[y], y*PwD2'[y]}, {y, 0.1, 1000000}, 
   PlotRange -> {0.01, 10}, PlotStyle -> {{Brown}, {Dashed, Brown}}];
P3 = LogLogPlot[{PwD3[y], y*PwD3'[y]}, {y, 0.1, 1000000}, 
   PlotRange -> {0.01, 10}, PlotStyle -> {{Purple}, {Dashed, Purple}}];

Show[P1, P2, P3]

The problem is that it takes a very long time to compute these codes... Since I don't have good programming skills, could anyone propose an alternative way to computing this plot?

thanks

share|improve this question
3  
Please do not post code that has been typeset. It is unreadable, and this case uninterpretable by Mathematica. Wrap the above expression in InputForm, and post the result. –  rcollyer Mar 9 '12 at 19:03
1  
The expression you posted got corrupted and cannot be interpreted correctly by Mathematica any more. It is important to post as input form (select cell, right click, copy as -> input text), otherwise we will not be able to read your code correctly. –  Szabolcs Mar 9 '12 at 19:18
1  
The inner integral has extreme high values for low values of tau. They are imaginary as well, so you can't plot them. I don't think there's an alternative. –  Sjoerd C. de Vries Mar 9 '12 at 20:27
    
That’s the problem... This is a well know equation published in 1975 within several graphics (with different limits of the integral in space: 50, 100, 1000) but it seems that nobody can find out a way to do it… –  Bruno Rangel Mar 10 '12 at 17:20
1  
Could you post a reference to this equation? It might be that you have made a typo while entering it. –  Sjoerd C. de Vries Mar 14 '12 at 11:05
add comment

2 Answers

up vote 3 down vote accepted

You could do the indefinite integral and put in the limits afterwards. E.g., on Linux I get a speed-up of about 14 for one integral (I have no time to do more right now ...)

Mathematica 8.0 for Linux x86 (64-bit)
Copyright 1988-2011 Wolfram Research, Inc.

In[1]:= !!i
Print["timing 1: ", 
First @ AbsoluteTiming[
  Pw3 = ExpandAll[(Exp[-0.25/t]/t)*
     Integrate[
      Exp[-((Tan[60*Degree]^2*
           (z + 30)^2)/(4*t))]*
       (1 + 2*Sum[Exp[-(n^2*Pi^2*t)/
            100^2]*Cos[0.8*n*Pi]*
           Cos[n*(Pi/2) - n*Pi*
            (z/100)], {n, 1, 100}]), 
      {z, -50, 50}]]; ]
];
Print["timing 2 : ", 
First @ AbsoluteTiming[
  nPw3 = Expand[(Exp[-0.25/t]/t)*
     ((#1 /. z -> 50) - (#1 /. 
         z -> -50) & )[
      ParallelMap[(ExpandAll[Integrate[#1, z]] & ),
       (TrigToExp[Expand[
          Exp[-((Tan[60*Degree]^2*
            (z + 30)^2)/(4*t))]*
           (1 + 2*Sum[Exp[-(n^2*Pi^2*
            t)/100^2]*Cos[0.8*n*Pi]*
            Cos[n*(Pi/2) - n*Pi*
            (z/100)], {n, 1, 
            100}])]])]] ]; ]
]
Print["check ", Chop[Pw3 - nPw3]]

In[1]:= <<i
timing 1: 75.121517
Launching kernels...       Mathematica 8.0 for Linux x86 (64-bit) Copyright 1988-2011 Wolfram Research, Inc.                                                                                                            
timing 2 : 5.255553                                                                                                                    
check 0                                                                                     
share|improve this answer
    
Thanks Rolf Mertig! It really helped! I didn't know about Linux speed up... Once again, thanks a lot! –  Bruno Rangel May 25 '12 at 12:05
add comment

Sjoerd’s comment on your question is spot on, but let’s explain a bit more. We will have a look at the integrand, i.e. the expression inside the outer integral:

$Assumptions = \[Tau] > 0 && n > 0;
f = (Exp[-(1/(4*\[Tau]))]/\[Tau])*
  Integrate[
   Exp[-((Tan[45*Degree]^2 - (z + 30)^2)/(4*\[Tau]))]*(1 + 
      2*Sum[Exp[-((n^2*Pi^2*\[Tau])/10000)]*Cos[8*n*Pi/10]*
         Cos[(n*Pi)/2 - (n*Pi*z)/100], {n, 1, 4}]), {z, -50, 50}]

I have removed the numerical coefficient 0.005, and rewritten your floating-point value 0.8 by the fraction 8/10. This gives a long expression for the integrand:

enter image description here

If you check the values of a few terms in the [0.1,2.0] range, you'll see that they are very big (and in fact overflow Mathematica’s precision limit) and have non-zero imaginary parts:

enter image description here

In fact, it only takes reasonable values for large values of tau:

In[28]:= f /. \[Tau] -> 20.
Out[28]= -9.16719*10^32 + 2.55438*10^17 I

but even there the imaginary part is not zero.


In short, the expression is complex, so you cannot plot it easily. Moreover, it diverges very fast near zero, so its evaluation is very hard. We can even quantify how fast its real part diverges by plotting Log@Log@Abs@Re@f vs. Log[tau]:

ListPlot@Table[{Log@\[Tau], Log@Log@Abs@N@Re[f]} /. \[Tau] -> 10^x, {x, -1, -5, -0.2}]

enter image description here

Conclusion: the real part of your integrand diverges as $\large \exp(t^{-\alpha})$.

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.