# Speeding up the time of integration calculation

I am going to integrate from Chebyshev orthogonal functions that is in the form of:

  u1[i_, j_] = (1 + 2*x/a)*(1 - 2*x/a)*(1 + 2*y/b)*(1 - 2*y/b)*
Cos[(i - 1)*ArcCos[2*x/a - 1]]*Cos[(j - 1)*ArcCos[2*y/b - 1]];
u2[K_, P_] = Sum[u1[i, j]^2, {i, K}, {j, P}];


Where I and j are in the range of 1 to 100. Because i*j=10^4, ten thousand integrations should be computed, so, it is very important to speed up the integrations.

The integration time is very long especially when the amount of i and j increase. For example, when i and j are from 1 to 5, we have

    a = 0.5;
b = 0.7;
n = 5;
z = 5;

Table[NIntegrate[u2[K, P], {x, 0, a}, {y, 0, b}], {K, 1, n}, {P, 1,
z}] // Timing


The result is=

  {154.063, {{0.822889, 1.32889, 1.59296, 1.95275, 2.33988}, {1.32889,
2.14603, 2.57249, 3.15352, 3.77868}, {1.59296, 2.57249, 3.08369,
3.78018, 4.52957}, {1.95275, 3.15352, 3.78018, 4.63398,
5.55264}, {2.33988, 3.77868, 4.52957, 5.55264, 6.65341}}}


The elapsed time is 154 second and during the computations. My computer is old, maybe in another computer, fewer time is displayed.

you see when i and j are from 1 to 5, the elapsed time is too long, Now consider how long it will take if the i and j are from 1 to 100. Any suggestion for speeding up the integrations will be very appreciated. Thank you

• Look up NumericalDifferentialEquationAnalysis/ref/GaussianQuadratureWeights in the help. Or on the web: https://reference.wolfram.com/language/NumericalDifferentialEquationAnalysis/ref/GaussianQuadratureWeights.html Feb 23, 2023 at 16:59
• You might separate the x and y part of the intergration Feb 23, 2023 at 17:05
• In your modified code you define scalar u=... but later use array u[k,p] ??? Feb 24, 2023 at 18:28
• Sorry, I corrected it. Feb 24, 2023 at 18:48
• Your "modification" gives a new question. The original content, my answer is related to, cannot be seen anymore Feb 24, 2023 at 21:25

The integral you are looking for might be separated into x-part and y-part. Additionally introducing new variables \[Xi]=x/a, \[Eta]=y/b gives the following expression for the integral

a b Integrate[(1 -(2 \[Xi])^2) Cos[(i - 1) ArcCos[2 \[Xi] - 1]], {\[Xi], 0, 1}]*Integrate[(1 -(2 \[Eta]) ^2) Cos[(i - 1) ArcCos[2 \[Eta] - 1]], {\[Eta], 0, 1}]

It's only necessary to evaluate one part

int=Table[NIntegrate[(1 -(2 \[Xi])^2) Cos[(i - 1) ArcCos[2 \[Xi] - 1]], {\[Xi], 0, 1},Method -> "LocalAdaptive"],{i,1,1000}]; //AbsoluteTiming


Evaluation only spents 20s for 1000x1000 integrals!

The Matrix of integrals i,j follow to

intij=a b Outer[Times ,int,int];


Hope it helps!

• Dear Ulrich Neumann, Thank you very much for your answer. Unfortunately, I cannot separate the integration into x- and y- parts manually. because, I am doing some Algebric manipulation on the original function (like multiplying, power etc. ) and the result should be integrated. 1- Is there any command that separate into x- and y- parts? 2- With these conditions, how should we change your suggested code? Thank you. Feb 24, 2023 at 8:08
• @PooyaAzizi You should modify your question. Only multiplying or power operations keep the possibilty of separation I think. Feb 24, 2023 at 9:41
• You are right. I apologize. Could you please help with this conditions? Feb 24, 2023 at 11:40
• Please show the modified integrand! Feb 24, 2023 at 12:14
• I reviewed and corrected the post Feb 24, 2023 at 17:25

Combining @UlrichNeumann's solution with analytic integration we get instantaneous results (0.6 seconds for 1000×1000 integrals) that remain accurate for very large $$i$$ and $$j$$:

g[i_?EvenQ] = 2/((i + 1) (i - 3));
g[i_?OddQ] = (i^2 - 2 i - 2)/((i + 2) i (i - 2) (i - 4));

a = 0.5;
b = 0.7;
G = Table[g[i], {i, 0, 10}];
intij = a * b * KroneckerProduct[G, G]

(*    {{0.155556, 0.0777778, 0.155556, 0.0155556, -0.0933333, -0.0288889, -0.0222222, -0.00814815, -0.0103704, -0.00410774, -0.00606061},
{0.0777778, 0.0388889, 0.0777778, 0.00777778, -0.0466667, -0.0144444, -0.0111111, -0.00407407, -0.00518519, -0.00205387, -0.0030303},
...
{-0.00606061, -0.0030303, -0.00606061, -0.000606061, 0.00363636, 0.00112554, 0.000865801, 0.00031746, 0.00040404, 0.000160042, 0.000236128}}    *)

• Dear roman, Thank you very much for your answer. because my function I integrate in not exactly the mentioned function, I change the post because the function couldn't separate to x- and y- parts. please see the function again and suggest your opinion. Thank you very much Feb 24, 2023 at 17:44

ParallelTable[NIntegrate[u2[K, P], {x, 0, a}, {y, 0, b}], {K, 1, n}, {P, 1, z}] // Timing