# Question about numerical integration of a double integral

I am trying to numerically integrate the following double integral in Mathematica for different values of t.

NIntegrate[(Xc[x]-I*Xs[x])*Exp[-I*x*t]*(Exp[I*y*t]/(y - x))*(g[y]*(I*y)*Exp[I*(x - y)*4]-g[x]*(I*x)), {y,1/10, 1}, {x, 1/10, 1}]

where

g[x_] := (23634.78881895757 + 7991.282962370622 I) - 1/x^12 (3.2553450760837677*^-10 - 1.7205915867327768*^-10 I) +1/x^10(1.6650398308338188*^-6 - 6.489989455356209*^-7 I) -1/x^9 (0.00009972422588178764 - 0.00003189730953401248 I) +1/x^8 (0.003083617887888503 - 0.0007696772987544534 I) -1/x^7 (0.06188400494926536 - 0.011082932555220979 I) +1/x^6 (0.8801611453125757 - 0.09543895809427982 I) -1/x^5 (9.282588939198732 - 0.35028669115085737 I) + 1/x^4 (74.49845808376074 + 2.458151790347552 I) -1/x^3 (461.69296523620704 + 48.05630698049144 I) +1/x^2 (2223.059410368775 + 392.556496017936 I) -1/x (8294.699448129079 +2097.6389759096696 I) - (49568.64394710974 +22060.93369851312 I) x + (69508.1296778248 +42984.2468547471 I) x^2 - (43449.9068787156 +53342.07572439036 I) x^3 - (50381.56033189651 -25301.12279548736 I) x^4 + (140618.56129587648 +37042.94278156334 I) x^5 - (90676.20127060641 +65871.11513129353 I) x^6 - (93845.15451151169 - 4017.80164549626 I) x^7 + (166126.617591954 + 73202.15594421528 I) x^8 + (24527.684407449655 -32143.939238180323 I) x^9 -(197006.64964445296 +78163.6737619927 I) x^10 + (39257.68452550401 +65252.70160939108 I) x^11 +195196.75053355526 +56540.798084673275 I) x^12 - (79145.20169999577 +  63057.97582026051 I) x^13 - (200259.17575548333 + 60782.92779231145 I) x^14 + (157866.57099768255 +86296.30791010859 I) x^15 + (152723.33096024842 + 38811.36333143129 I) x^16 - (308798.67979950836 + 129973.26458647089 I) x^17 + (209954.01964860896 +
100223.42774558903 I) x^18 - (69103.28262535635 + 35570.22891265897 I) x^19 + (9287.613868933946 +5052.39844141709 I) x^20

and

Xc[x_] :=-1.9822786347952313 +4.232100200840357*^-15/x^12-2.3091024542586374*^-11/x^10+1.4694235760892122*^-9/x^9-4.947763862180016*^-8/x^8+1.1138000174717181*^-6/x^7-0.000018407606213391123/x^6 + 0.0002351076696658846/x^5 - 0.0023966842802706245/x^4 +0.019917254972612142/x^3 -0.1367069136297532/x^2 + 0.7784581812138608/x+7.061989141343013*x-12.318273951718858*x^2+9.095575318222767*x^3+2.1169232598842935*x^4-13.127424409874147*x^5+9.759769792477377*x^6+7.149060936983623*x^7 - 14.143044062774418*x^8-1.4297912717457828*x^9 +15.30960928306342*x^10-3.2083738330017795*x^11-14.270695822044644*x^12+5.716658723315534*x^13+14.119588947432252*x^14-10.942024599532623*x^15-10.484868293647892*x^16 + 20.908398165073194*x^17-14.066767979105023*x^18+4.5905538331721605*x^19 -0.6127869681947318*x^20

and

Xs[x_] := -0.5543449775870041 + 5.152295531964746*^-15/x^12 -3.0396613988535765*^-11/x^10 +
1.989603080424189*^-9/x^9 - 6.814705925440601*^-8/x^8 +1.537223625582891*^-6/x^7 -
0.000024950767687590638/x^6 + 0.00030468848164531885/x^5 -0.0028637592671235373/x^4 +
0.020840985554216367/x^3 - 0.11537446279155367/x^2 +0.44531977062138367/x +
0.9647558248886902*x + 4.23106215830098*x^2 -9.162586748844152*x^3 +
7.397388545365722*x^4 + 1.5118748118502208*x^5 -8.872844269214584*x^6 +
3.788492865641893*x^7 + 7.16726040461682*x^8 - 5.715862319293517*x^9 -
6.616070309864808*x^10 + 8.294276652580344*x^11 + 3.2609579031597744*x^12 -
6.457044690769427*x^13 - 3.480270133569723*x^14 +7.221666205344835*x^15 +
1.619574214112665*x^16 - 8.89035897402073*x^17 + 7.133762722587951*x^18 -
2.549042217230423*x^19 + 0.36068707813056955*x^20

I fitted some data points to get these polynomials. It seems that there is a singularity at line y=x, but if we look at the entire integral with y=x, the integral is finite. The problem is that Mathematica is very slow in computing this integral for large t, say t=150 or 200. I tried different methods of integration like LevinRule but it is still slow and also gives me the error "NIntegrate converging slowly ...". I was wondering if you can help me with this. Thank you!

• It seems that all of your terms are multiples of $\frac{x^{\text{kx}} x^{\text{ky}} e^{-i t (x-y)}}{y-x}$. If so, then one can use Integrate rather than NIntegrate as in Integrate[(E^(-I t (x - y)) x^kx x^ky)/(-x + y) /. kx -> 3 /. ky -> -3, {x, 1/10, 1}, {y, 1/10, 1}, Assumptions -> t > 0] which results in (I (-20 - 9 \[Pi] t + 20 Cos[(9 t)/10] + 18 t SinIntegral[(9 t)/10]))/(10 t).
– JimB
May 30, 2023 at 19:57
• @JimB I don't think they are multiples of $x^{kx} y^{kx}e^{-it(x-y)}/y-x$. May 30, 2023 at 21:07

This is just an extended comment: Maybe a symbolic integration might speed things up as all terms to be integrated are proportional to a small number of patterns:

g[x_] := cg[0] + Sum[cg[-i]/x^i, {i, 1, 12}] + Sum[cg[i] x^i, {i, 1, 20}]
Xc[x] := cXc[0] + Sum[cXc[-1]/x^i, {i, 1, 12}] + Sum[cXc[i] x^i, {i, 1, 20}]
Xs[x] := cXs[0] + Sum[cXs[-1]/x^i, {i, 1, 12}] + Sum[cXs[i] x^i, {i, 1, 20}]
integrands = List @@ ((Xc[x] - I*Xs[x])*Exp[-I*x*t]*(Exp[I*y*t]/(y - x))*(g[y]*(I*y)*Exp[I*(x - y)*4] - g[x]*(I*x)) // ExpandAll);

Now "condense" the patterns of the individual terms:

terms = integrands /. cg[k_] -> 1 /. cXc[k_] -> 1 /. cXs[k_] -> 1 /. I Exp[z_] -> Exp[z];
terms = terms // FullSimplify;
terms = terms /. Times[Complex[0, 1], Power[E, z1_], z2_, z3_] -> Times[Power[E, z1], z2, z3];
terms = terms /. x^k_ -> x^kx /. y^k_ -> y^ky /. Times[-1, z_] -> Times[z] // DeleteDuplicates

I'll make further reductions by inspection:

Now some of the terms when given specific values of kx and ky do not integrate over {x, 1/10, 1} and {y, 1/10, 1} so that issue might be hidden when using numerical integration.

• Thanks a lot for your detailed answer! Jun 3, 2023 at 2:07