1
$\begingroup$

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

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

3 Answers 3

4
$\begingroup$

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!

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Feb 24, 2023 at 8:08
  • $\begingroup$ @PooyaAzizi You should modify your question. Only multiplying or power operations keep the possibilty of separation I think. $\endgroup$ Feb 24, 2023 at 9:41
  • $\begingroup$ You are right. I apologize. Could you please help with this conditions? $\endgroup$ Feb 24, 2023 at 11:40
  • $\begingroup$ Please show the modified integrand! $\endgroup$ Feb 24, 2023 at 12:14
  • $\begingroup$ I reviewed and corrected the post $\endgroup$ Feb 24, 2023 at 17:25
3
$\begingroup$

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}}    *)
$\endgroup$
1
  • $\begingroup$ 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 $\endgroup$ Feb 24, 2023 at 17:44
0
$\begingroup$

How about ParallelTable?

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

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.