# How to improve the speed of 'NIntegrate'?

Problem 1
My code is given as follows

Clear["Global*"]
d = 0.8;
p = 0.26;
a = 2;
b = 0.03;
int1[s_?NumericQ, a_] :=
NIntegrate[(2 (k^4 - 1) (1 + a b (k^2 + k^-2 - 2)^(a - 1)))/(
k^3 (k^2 - 1)), {k, ((d + s^2)/(1 + d))^(1/2), s}]
int2[x_?NumericQ, a_] := NIntegrate[int1[s, a], {s, 1, x}]
Plot3D[int2[x, a] - (1/2) p*(x^2 - 1) + (1/4) x^2*
Log[(d + x^2)/x^2] y^2, {x, 0.5, 9.5}, {y, -1.5, 1.5},
PlotRange -> {-0.3, 0.15}] // AbsoluteTiming


It runs about 209s in my computer. How to improve the speed? Any suggestion is much appreciated!
Problem 2
When we vary a, the code is

Clear["Global*"]
d = 0.8;
p = 0.26;
b = 0.03;

ContourPlot3D[
NIntegrate[
NIntegrate[(2 (k^4 - 1) (1 + a b (k^2 + k^-2 - 2)^(a - 1)))/(
k^3 (k^2 - 1)), {k, ((d + s^2)/(1 + d))^(1/2), s}], {s, 1, x}] -
1/2 p*(x^2 - 1) + 1/4 x^2*Log[(d + x^2)/x^2] y^2 == 0.029744, {x,
0.5, 13}, {y, -1.5, 1.5}, {a, 1.95, 2.08}]


For this problem, I can't get a result. How can we deal with this problem?

• It takes so much time, because int1 does integration for overlapping k ranges multiple times for different s. If you want to do it with NIntegrate, form a double integral like doubleint[x_, a_] := NIntegrate[(2 (k^4 - 1) (1 + a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)), {s, 1, x}, {k, ((d + s^2)/(1 + d))^(1/2), s}]  . Attention, the integration which depends on s must be more outside. This is 20 times faster. Oct 8 '20 at 13:04

## 2 Answers

If you want to do or you must do integration numericaly (e.g. a is not an integer), you can get antiderivatives as interpolation functions with NDSolve.

Clear["Global*"];
d = 0.8;
p = 0.26;
a = 2;
b = 0.03;
f[k_] = (2 (k^4 - 1) (1 +
a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)
) //   Simplify;


Get the first antiderivative i1sol. The intial conditions are not decisive, only differences count as shown in the answer of @cvgmt.

(i1sol =
i1 /. First@
NDSolve[{i1'[k] == f[k], i1[1/2] == 0}, i1, {k, 1/2, 10}];

i2sol =
i2 /. First@
NDSolve[{i2'[s] == i1sol[s] - i1sol[((d + s^2)/(1 + d))^(1/2)],
i2[1] == 0}, i2, {s, 1/2, 10}];

Plot3D[((i2sol[x] - i2sol[1]) - (1/2) p*(x^2 - 1) + (1/4) x^2*
Log[(d + x^2)/x^2] y^2), {x, 0.5, 9.5}, {y, -1.5, 1.5},
PlotRange -> {-0.3, 0.15}, ImageSize -> 300]
) // Timing


This is even faster than the analytical approach.

Edit to Problem 2

In higher versions of MMA you would do this, when using the above method, with ParametricNDSolve  with a  as parameter. Since i use version 8.0, i have to do it this way:

d = 0.8;
p = 0.26;
b = 0.03;
f[k_, a_] = (2 (k^4 - 1) (1 +
a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)) // Simplify

i1sol = i1 /.
First@NDSolve[{Derivative[1, 0][i1][k, a] == f[k, a],
i1[1/2, a] == 0}, i1, {k, 1/3, 13}, {a, 1.95, 2.08}]

i2sol = i2 /.
First@NDSolve[{Derivative[1, 0][i2][s, a] ==
i1sol[s, a] - i1sol[((d + s^2)/(1 + d))^(1/2), a],
i2[1, a] == 0}, i2, {s, 1/2, 13}, {a, 1.95, 2.08}]

ContourPlot3D[(i2sol[x, a] - i2sol[1, a]) - 1/2 p*(x^2 - 1) +
1/4 x^2*Log[(d + x^2)/x^2] y^2 == 0.029744, {x, 0.5, 13}, {y, -1.5,
1.5}, {a, 1.95, 2.08}]


• Thanks, It is a smart method! And do you have any ideas of a similar for Problem 2 I added in my question when a is not fixed? Oct 9 '20 at 0:36

Use Newton-Leibniz Formula $$\int_a^b f(t)\,dt=F(b)-F(a)$$

We use F[x_]=Integrate[f[x],x] and subtract F[b] and F[a] instead of Integrate[f[x],{x,a,b} to calculate $$\int_a^b f(t)\,\mathrm{d}t$$ since Integrate[f[x],x] is easy to calculate.

Clear["Global*"]
d = 0.8;
p = 0.26;
a = 2;
b = 0.03;
intA[k_] =
Integrate[(2 (k^4 - 1) (1 +
a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)), k]
int1[s_, a_] = intA[s] - intA[((d + s^2)/(1 + d))^(1/2)];
intB[s_] = Integrate[int1[s, a], s];
int2[x_, a_] = intB[x] - intB[1];
Plot3D[int2[x, a] - (1/2) p*(x^2 - 1) + (1/4) x^2*
Log[(d + x^2)/x^2] y^2, {x, 0.5, 9.5}, {y, -1.5, 1.5},
PlotRange -> {-0.3, 0.15}] // AbsoluteTiming


• It is a useful method. But I wonder if there are other ways to involve the 'NIntegrate', since that when I evaluate any points in the figure, for instance, 'int2[1.5, a] // AbsoluteTiming', it gives the result very fast. Oct 8 '20 at 12:14
• This method seems not to work if a is not fixed. The code is given in Problem2, do you have any suggestions? Thanks! Oct 9 '20 at 0:33