# Plotting a double integral with variable limits

I need to plot the following double integral

$\sigma(t)=\sum_{p=2}^N\int_0^td\tau_1\int_0^td\tau_2e^{-|\tau_1-\tau_2|/\tau_A}\tau_1^{\alpha-1}\tau_2^{\alpha-1}E_{\alpha,\alpha}[-(\tau_1/\tau_p)^{\alpha}]E_{\alpha,\alpha}[-(\tau_2/\tau_p)^{\alpha}]$

where the E is the Mittag-Leffler function.

I tried to do it the following way:

s[t_?NumericQ] := Sum[NIntegrate[
Exp[-Abs[t2 - t1]/tA] * t1^(a-1) * t2^(a-1) *
MittagLefflerE[a, a, -(t1/tp[[p]])^a] * MittagLefflerE[a, a, -(t2/tp[[p]])^a],
{t1, 0, t}, {t2, 0, t}
], {p, 2, n}];

LogLogPlot[s[t], {t,0.001,1000}]


where

a = 0.7;
tA = 0.01;
tp = {0.1, 0.2, 0.3, 0.4, 0.5};
n = Length[tp];


Unfortunately, this doesn't produce any result for me, Mathematica just keeps calculating. Does anybody know how to get a plot of this function? Thanks a lot!

• How are a, tA, tp and n defined? – Michael E2 Oct 25 '16 at 10:59
• a=0.7, tA=0.01, tp is list of eigenvalues of a matrix and n is the length of this list. For simplicity we can just use tp={0.1,0.2,0.3,0.4,0.5} or something like that – jones Oct 25 '16 at 11:24
• Should tp be tp[[p]]? – Michael E2 Oct 25 '16 at 11:41
• Yes, sorry, I fixed it. – jones Oct 25 '16 at 12:08

Let me tell you something you already know. Mathematica doesn't work on magic.

Running the code for just one value of t:

AbsoluteTiming[
Sum[NIntegrate[
Exp[-Abs[t2 - t1]/tA]*t1^(a - 1)*t2^(a - 1)*
MittagLefflerE[a, a, -(t1/tp[[p]])^a]*
MittagLefflerE[a, a, -(t2/tp[[p]])^a], {t1, 0, 1}, {t2, 0,
1}], {p, 2, n}]]


{124.122, 0.028238}

That's over two minutes, for a range of $\left[0,1\right]$, you can imagine how long this takes for $\left[0,1000\right]$.

As @Michael-E2 suggested transposing the Integral[] and the Sum[] could tweak performance, as well as standard methods such as reducing AccuracyGoal, but this will not make the crunching time manageable.

How would I plot this? Book time at a supercomputer.

• You can save some time by lowering AccuracyGoal (at least for plotting), but I couldn't set it less than AccuracyGoal -> 6 and get a seemingly accurate result. It cut the time only in 1/2 or 1/3. Still quite slow. -- (Oh, and I transposed the integral and the sum.) – Michael E2 Oct 25 '16 at 14:11