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 have this equation:

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}]

Mathematica graphics Now, to generate 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, due to the summation limit (1 to 100)! ... Could anyone propose an alternative way to computing this plot?

share|improve this question
    
Is there any way you can simplify the problem? Is there, by chance, an analytical integral you could pre-compute? –  Verbeia May 23 '12 at 12:36
    
When posting blocks of code, please prefix it with four spaces. The easiest way is to select the block of code and press <kbd>Ctrl</kbd>-k. Only use backward quotes for inline code. –  Heike May 23 '12 at 12:38
    
Verbeia, unfortunately I believe it doesn’t. I already have made simplification to express like this... –  Bruno Rangel May 23 '12 at 12:56
    
Check out Rolf Mertig's answer on your other post. In summary, he maps Integrate across each term in the sum which can then be done in parallel. Additionally, he calculates them as indefinite integrals, substituting the end-points after the integrals have been performed. All together, he sees a dramatic speed up. –  rcollyer May 23 '12 at 13:39
3  
Please post a short and more to the point question. As far as I can tell, everything above the second code block isn't relevant. Second, you want to speed up PD1. It's much easier to work out what's going on for us if you ask that, instead of posting your notebook and expecting us to work out where the bottleneck is. I did it because I was stuck waiting for someone, but it took me 30min to reorganize your code, 5min to speed it up (although not enough I'm afraid). –  acl May 23 '12 at 13:41
add comment

2 Answers

up vote 1 down vote accepted

Here's a possible approach. The basic idea is to Integrate a single term of the sum over z, and do the sum afterwards. I use the undefined symbol tps to represent Tan[psi] as this seems to speed the integration up.

For the numerical integration, I apply Re to the function since we know the answer is real and there's no point having NIntegrate spend time worrying over tiny imaginary parts which shouldn't even be there. I've also split the integrations into individual sections rather than starting at 0 every time.

\[Psi]=60Degree;
(* the z integral is faster if we keep Tan[psi] undefined,
but we need to deal with the Tan[psi]=0 case separately*)
tanpart=If[Tan[\[Psi]]==0,1,Exp[-(tps^2(z+30)^2)/(4\[Tau])]];
(* this is the term in the sum *)
term=Cos[4n Pi/5]Cos[n Pi/2-n Pi z/100]Exp[-n^2 Pi^2 \[Tau]/10000]tanpart;
(* integrate over z *)
zint=Integrate[term,{z,-50,50}];
(* this is the n=0 term *)
zint0=Integrate[tanpart,{z,-50,50}];
(* now do the sum and insert the value of Tan[psi] *)
summed=(0.005/\[Tau])Exp[-1/(4\[Tau])]*(zint0+2Sum[zint,{n,1,100}])/.tps->Tan[\[Psi]];
(* y values *)
yvals={0,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};
(* integration limits *)
intervals=Most@Transpose@{yvals,RotateLeft[yvals]};
(* do the numerical integrations *)
sections=NIntegrate[Re@summed,{\[Tau],#[[1]],#[[2]]},PrecisionGoal->5]&/@intervals;
(* final results *)
values=Accumulate[sections];   
share|improve this answer
    
Hey Simon, thanks! Last week I was trying to do exactly that... But I wasn't able to find the proper codes for it. This really helps. Also, I'm trying to run mathematica in Linux, it seem to speed up a lot... thanks! –  Bruno Rangel May 25 '12 at 12:13
add comment

Have you tried defining your functions like this:

ClearAll[f];
f[t_, z_] = 
  Total@Table[
    Exp[-((n^2)*(Pi^2)*t)/(100^2)]*Cos[0.8*n*Pi]*
     Cos[(n*Pi/2) - (n*Pi*z/100)], {n, 1, 100}];

Pw1[t_?NumericQ] := (Exp[-0.25/t]/t)*
   NIntegrate[
    Exp[-(((Tan[0 Degree]^2)*(z + 30)^2)/(4*t))]*(1 + 
       2*f[t, z]), {z, -50, 50}];

Pw2[t_?NumericQ] := (Exp[-0.25/t]/t)*
   NIntegrate[
    Exp[-(((Tan[30 Degree]^2)*(z + 30)^2)/(4*t))]*(1 + 
       2*f[t, z]), {z, -50, 50}];

Pw3[t_] := (Exp[-0.25/t]/t)*
   NIntegrate[
    Exp[-(((Tan[60 Degree]^2)*(z + 30)^2)/(4*t))]*(1 + 
       2*f[t, z]), {z, -50, 50}];

then if PD1[y_] := 0.005*NIntegrate[Pw1[t], {t, 0, y}, MaxRecursion -> 20, PrecisionGoal -> 5], PD[0.1] takes 40s on my machine.

That's all I have time for right now I'm afraid.

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.