# Plot Indefinite integral

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


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

• 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. Mar 9 '12 at 19:03
• 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. Mar 9 '12 at 19:18
• 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. 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… Mar 10 '12 at 17:20
• Could you post a reference to this equation? It might be that you have made a typo while entering it. Mar 14 '12 at 11:05

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)

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

• Thanks Rolf Mertig! It really helped! I didn't know about Linux speed up... Once again, thanks a lot! May 25 '12 at 12:05

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: 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: 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}]  Conclusion: the real part of your integrand diverges as$\large \exp(t^{-\alpha})\$.