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
InputForm
, and post the result. $\endgroup$