-4
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ Is there any way you can simplify the problem? Is there, by chance, an analytical integral you could pre-compute? $\endgroup$
    – Verbeia
    Commented May 23, 2012 at 12:36
  • $\begingroup$ 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. $\endgroup$
    – Heike
    Commented May 23, 2012 at 12:38
  • $\begingroup$ Verbeia, unfortunately I believe it doesn’t. I already have made simplification to express like this... $\endgroup$
    – Bruno Rangel
    Commented May 23, 2012 at 12:56
  • $\begingroup$ 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. $\endgroup$
    – rcollyer
    Commented May 23, 2012 at 13:39
  • 3
    $\begingroup$ 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). $\endgroup$
    – acl
    Commented May 23, 2012 at 13:41

2 Answers 2

Reset to default
1
$\begingroup$

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];
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – Bruno Rangel
    Commented May 25, 2012 at 12:13
1
$\begingroup$

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.

Improve this answer
$\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.