Consider some dataset corresponding to the grid x1, x2, function[x1,x2]:

function[x1_, x2_] = Exp[-x2/25]*Cos[20*x1]^2*Sin[x1];
dataset = 
  Flatten[Table[{x1, x2, function[x1, x2]}, {x1, 
      x, {x, -5., Log10[1/20. Pi/2], (Log10[1/20. Pi/2] + 5.)/
       20}]}, {x2, 
      x, {x, Log10[5.3], Log10[200.], (Log10[200.] - Log10[5.3])/
       20}]}], {1, 2}];
{x1min, x1max} = MinMax[dataset[[All, 1]]]
{x2min, x2max} = MinMax[dataset[[All, 2]]]

I interpolate it in two ways:

int1[x1_, x2_] = 10^(Interpolation[{Log10[#[[1]]], Log10[#[[2]]], Log10[#[[3]] + 10^-90]}&/@dataset, InterpolationOrder -> 1][Log10[x1], Log10[x2]]);
int2[x1_, x2_] = Interpolation[{#[[1]], #[[2]], #[[3]]}&/@dataset, InterpolationOrder -> 1][x1, x2];

The integration of the int1 is much slower than int2:

testfunction[x1_, x2_] = Sqrt[(1 - Cos[x1])]*Exp[5/x2];
NIntegrate[int1[x1, x2], {x1, x1min, x1max}, {x2, x2min, x2max}, 
  Method -> "InterpolationPointsSubdivision"] // AbsoluteTiming
NIntegrate[int2[x1, x2], {x1, x1min, x1max}, {x2, x2min, x2max}, 
  Method -> "InterpolationPointsSubdivision"] // AbsoluteTiming
  int1[x1, x2]*testfunction[x1, x2], {x1, x1min, x1max}, {x2, x2min, 
   x2max}, Method -> 
   "InterpolationPointsSubdivision"] // AbsoluteTiming
  int2[x1, x2]*testfunction[x1, x2], {x1, x1min, x1max}, {x2, x2min, 
   x2max}, Method -> 
   "InterpolationPointsSubdivision"] // AbsoluteTiming

{0.960714, 0.0146339}

{0.0011156, 0.0189933}

{1.01974, 0.000414336}

{0.0718938, 0.000659965}

On the other hand, int1 has better behavior in case if function changes fast between the neighboring values of x1, x2 in the grid.

Switching to "AdaptiveMonteCarlo" integration method somehow improves the situation, but non-negligible difference in timing (up to a few times) remains (it's funny that it slows down much the integration of int2 without testfunction).

Could you please tell me how to speed up the integration with int1?

  • 1
    $\begingroup$ What does int1 have to do with int2 ? $\endgroup$ Commented Mar 15, 2023 at 17:53
  • $\begingroup$ Of course int1 is slower: In addition to evaluating an interpolation function, it has also to evaluate 10^... and Log10. $\endgroup$ Commented Mar 15, 2023 at 18:04
  • $\begingroup$ What do you actually try to achieve here? $\endgroup$ Commented Mar 15, 2023 at 18:05
  • $\begingroup$ This question seems ill-defined, because you're relying on Mathematica to give you an interpolation of a function, and it could be computing the two interpolations totally differently. One might have $1000$ times as many points as the other one. That would be my first guess as to why one is slow vs. fast. I think you need to provide more context for the question to be meaningful to others. What real data are you given to work with, and what are your options for processing it? I assume that you're finding one option is a lot slower than the other, but otherwise preferable. $\endgroup$
    – Myridium
    Commented Mar 15, 2023 at 19:43
  • $\begingroup$ @HenrikSchumacher : I would like to get the accurate interpolation in the domain where the distribution function decreases fastly (large x1 and x2), keeping also good performance. $\endgroup$ Commented Mar 15, 2023 at 21:13

1 Answer 1


I suggest that you try to use the substitution formula:

$$ \int_a^b \int_c^d f( \varphi(x),\psi(y)) \, \varphi'(x) \, \psi'(y) \, \mathrm{d}x \, \mathrm{d}y = \int_{\varphi(a)}^{\varphi(a)} \int_{\psi(c)}^{\psi(d)} f(u,v) \, \mathrm{d}u \, \mathrm{d}v. $$

So you could apply this to $\varphi(t) = \psi(t) = \log_{10}(t)$ or $\varphi(t) = \psi(t) = 10^{t}$. (I don't know at the moment which direction would be suitable for you.)

If $\varphi'(x) \neq 0$ and $\psi'(y) \neq 0$ for all $x \in [a,b]$ and $y \in [c,d]$, then you may define

$$ g(u,v) := f(u,v) \, \varphi'(\varphi^{-1}(u)) \, \psi'(\psi^{-1}(v)), $$ which leads to the following identity: $$ \int_a^b \int_c^d g( \varphi(x),\psi(y)) \, \mathrm{d}x \, \mathrm{d}y = \int_{\varphi(a)}^{\varphi(a)} \int_{\psi(c)}^{\psi(d)} g(u,v) \, \frac{1}{\varphi'(\varphi^{-1}(u))} \, \frac{1}{\psi'(\psi^{-1}(v)))} \, \mathrm{d}u \, \mathrm{d}v. $$

The point is to apply this to the following function $g$:

ip = Interpolation[Transpose[{Log10[x], Log10[y], Log10[z + 10^-90]}],
InterpolationOrder -> 1];
g[x1_, x2_] := 10^ip[x1, x2];

Method -> "InterpolationPointsSubdivision" should allow for fast integration of the integral on the right. And if I understand correctly, it is the integral on the left that you want to compute in the end.

  • 2
    $\begingroup$ Glückwunsch zu den 100k! $\endgroup$
    – user21
    Commented Mar 16, 2023 at 5:23
  • $\begingroup$ Hey, danke! War 'ne Menge Arbeit... =) $\endgroup$ Commented Mar 16, 2023 at 6:25

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.