3
$\begingroup$

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

$\endgroup$
5
  • 3
    $\begingroup$ 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. $\endgroup$
    – rcollyer
    Mar 9, 2012 at 19:03
  • 1
    $\begingroup$ 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. $\endgroup$
    – Szabolcs
    Mar 9, 2012 at 19:18
  • 1
    $\begingroup$ 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. $\endgroup$ Mar 9, 2012 at 20:27
  • $\begingroup$ 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… $\endgroup$ Mar 10, 2012 at 17:20
  • 1
    $\begingroup$ Could you post a reference to this equation? It might be that you have made a typo while entering it. $\endgroup$ Mar 14, 2012 at 11:05

2 Answers 2

3
$\begingroup$

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                                                                                     
$\endgroup$
1
  • $\begingroup$ Thanks Rolf Mertig! It really helped! I didn't know about Linux speed up... Once again, thanks a lot! $\endgroup$ May 25, 2012 at 12:05
7
$\begingroup$

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})$.

$\endgroup$

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.