# Numerical integration converging too slowly

I must solve this integral which I suppose to be a very small number. How can I do? When I wrote this code:

r = (1.082)*10^8
h = (4.87)*10^(24)
Q[x_] := x - (0.206)*Sin[x] + (0.206)*(0.206)*Sin[x]*Cos[x]
R[x_] := ((57.91)*10^6)*(1 + 0.206*Cos[Q[x]])
B[x_] := (0.206)*Sin[2*x] - Sin[x]
A[x_] := ((0.62)*10^(-25))*((4.854)*Cos[Q[x]] - B[Q[x]]*Sin[Q[x]])
F[x_, y_] :=
R[x]*A[x] -
r*Cos[y]*((((57.91)*10^6)*A[x]*(1 - (11.93)*10^6)*
Cos[Q[x]]/R[x]) - ((-3.008)*10^(-25)) - ((57.91)*10^6)*
B[x]*((-0.052)*10^(-31))*Sin[Q[x]]) -
r*Sin[y]*(((0.6329)^10^(-25))*
Sin[Q[x]] + ((56.67)*10^6)*((-0.05196)*10^(-31))*B[x]*Cos[Q[x]])
G[x_, y_] := (1 - (2*
R[x]/r)*(Cos[y]*(((11.93)*10^6) + ((57.91)*10^6)*Cos[Q[x]])/
R[x]) + (R[x]/r)^2)^(-1.5)
u = ((0.01917)*10^(27))*((57.91)*10^6)
NIntegrate[(h/(((1.989)*10^(30)))*(r^3))*F[x, y]*G[x, y]*u, {x, 0,
2 Pi}, {y, 0, 2 Pi}, PrecisionGoal -> 50, MaxRecursion -> 50,
WorkingPrecision -> 100]


it gives me this error:**

"The precision of the argument function \
((3.44312*10^51\(<<1>>))/(1+0.286452\(<<1>>)^2-1.84843*10^-8\Cos[y]\(\
1.193*10^7+5.791*10^7\Cos[Plus[<<3>>]]))^1.5) is less than \
WorkingPrecision (100.)"


and this:

Numerical integration converging too slowly; suspect one of the \
following: singularity, value of the integration is 0, highly \
oscillatory integrand, or WorkingPrecision too small.

he global error of the strategy GlobalAdaptive has increased more \
than 2000 times. The global error is expected to decrease \
monotonically after a number of integrand evaluations. Suspect one of \
the following: the working precision is insufficient for the \
specified precision goal; the integrand is highly oscillatory or it \
is not a (piecewise) smooth function; or the true value of the \
integral is 0. Increasing the value of the GlobalAdaptive option \
MaxErrorIncreases might lead to a convergent numerical integration. \
NIntegrate \
obtained 8.\
2245804667122205619329698365928375872918387045667313082675253184241492\
7464246649828330519930578973435118997853216575354898819456603870372689\
813756795*10^42 and 3.\
8524159983189941460321064133370663169006624537848246120142620883724145\
3794096694780053810840247925397951565645812529096240278267666686786349\
959791842*10^48 for the integral and error estimates.

• The error explicitly suggests using higher WorkingPrecision. Have you tried doing just that? In particular, my understanding is that you must set WorkingPrecision to be at least as large as PrecisionGoal. Set WorkingPrecision -> 1000? – MarcoB Jun 11 '15 at 16:54
• i set PrecisionGoal -> 50, MaxRecursion -> 50, WorkingPrecision -> 100 – Deimos Jun 11 '15 at 17:10
• and he told me The precision of the argument function... is less than WorkingPrecision (100.) – Deimos Jun 11 '15 at 17:11
• i tried also PrecisionGoal -> 1000, MaxRecursion -> 1000, WorkingPrecision -> 2000 – Deimos Jun 11 '15 at 17:12
• OK, so what are the argument functions? Could you show how you defined them? – MarcoB Jun 11 '15 at 17:12

As noted in the Comments, this integration is plagued by precision problems. To proceed, factor the huge constant

(h/(((1.989)*10^(30)))*(r^3)) u
(* 3.44312*10^51 *)


from the integrand and then FullSimplify and Rationalize the functions F[x, y]*G[x, y]]

rat = Rationalize[FullSimplify[F[x, y]*G[x, y]], 0];


in order to achieve any possible cancellations symbolically (not many) and to increase the precision of the integrand to infinity. (The first error message in the Question was Mathematica complaining that the precision of the integrand was less (in fact, much less) than the requested WorkingPrecision -> 100.)

To gain some insight into the integrand,

Plot3D[rat, {x, 0, 2 Pi}, {y, 0, 2 Pi}, PlotRange -> All] Again as noted in the comments, the integrand is highly symmetric in both dimensions. This suggests that the four quadrants be summed symbolically before integration.

Plot3D[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx,
y -> yy}) + (rat /. {x -> xx,
y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx,
y -> 2 Pi - yy}), {xx, 0,  Pi}, {yy, 0, Pi}, PlotRange -> All] Observe that the modified integrand is some 15 orders of magnitude smaller (a huge cancellation) and also that it is highly irregular. It now can be integrated.

NIntegrate[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx,
y -> yy}) + (rat /. {x -> xx, y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx,
y -> 2 Pi - yy}), {xx, 0,  Pi}, {yy, 0, Pi},
MaxPoints -> 500000, WorkingPrecision -> 30, PrecisionGoal -> 9]
(* 3.94313*10^-9 *)


without error messages. Note that MaxPoints -> 500000, WorkingPrecision -> 30 is necessary to obtain reasonable accuracy, at the cost of a few minutes of computational time. Multiplying the answer by the huge constant factored out earlier gives 1.35767*10^43.

Improved computation

The question below prompted me to consider more carefully the second plot above and the subsequent integration. Applying FullSimplify to the argument of the second plot reduces its LeafCount from 887 to 192, and provides a much more accurate plot.

sim = FullSimplify[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx, y -> yy})
+ (rat /. {x -> xx, y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx, y -> 2 Pi - yy})];
Plot3D[sim, {xx, 0, Pi}, {yy, 0, Pi}, PlotRange -> All] Then,

NIntegrate[sim, {xx, 0, Pi}, {yy, 0, Pi}, MaxPoints -> 1000000,
WorkingPrecision -> 30, PrecisionGoal -> 9]
(* -2.25620530195937290197090585866*10^-18 *)


Multiplying the answer by the huge constant factored out earlier gives -7.76838*10^33.

• Anyone kindly explain this " Plot3D[(rat /. {x -> xx, y -> yy}) + (rat /. {x -> 2 Pi - xx, y -> yy}) + (rat /. {x -> xx, y -> 2 Pi - yy}) + (rat /. {x -> 2 Pi - xx, y -> 2 Pi - yy}), {xx, 0, Pi}, {yy, 0, Pi}, PlotRange -> All] " .Please if possible kindly explain what is happening here? – Gummala Navneeth Nov 25 '19 at 13:01
• @GummalaNavneeth Basically, I separated the integrand into four quadrants and combined them symbolically before performing the integration. Doing so is very effective in this case, because of the substantial symbolic cancellations that occur. Thanks for asking this question, because it caused me to realize that I should have applied FullSimplify to the new integrand to achieve much greater numerical precision. Please see the new section of my answer. – bbgodfrey Nov 26 '19 at 22:09

I think you have a couple of typos in your formulas. Once these are corrected, the integral evaluates fine:

1. There should be no backtick at the end of the definition of Q[x], and
2. in the definition of F[x,y] I think you should replace (0.6329)^10^(-25) with (0.6329)*10^(-25).

Once these are fixed, the integral comes out as $$1.35767\times 10^{43}$$:

r = (1.082)*10^8;
h = (4.87)*10^(24);
Q[x_] = x - (0.206)*Sin[x] + (0.206)*(0.206)*Sin[x]*Cos[x];
R[x_] = ((57.91)*10^6)*(1 + 0.206*Cos[Q[x]]);
B[x_] = (0.206)*Sin[2*x] - Sin[x];
A[x_] = ((0.62)*10^(-25))*((4.854)*Cos[Q[x]] - B[Q[x]]*Sin[Q[x]]);
F[x_, y_] = R[x]*A[x] - r*Cos[y]*((((57.91)*10^6)*A[x]*(1 - (11.93)*10^6)*
Cos[Q[x]]/R[x]) - ((-3.008)*10^(-25)) - ((57.91)*10^6)*
B[x]*((-0.052)*10^(-31))*Sin[Q[x]]) -
r*Sin[y]*(((0.6329)*10^(-25))*
Sin[Q[x]] + ((56.67)*10^6)*((-0.05196)*10^(-31))*B[x]*
Cos[Q[x]]);
G[x_, y_] = (1 - (2*R[x]/r)*(Cos[y]*(((11.93)*10^6) + ((57.91)*10^6)*Cos[Q[x]])/
R[x]) + (R[x]/r)^2)^(-1.5);
u = ((0.01917)*10^(27))*((57.91)*10^6);
NIntegrate[(h/(((1.989)*10^(30)))*(r^3))*F[x, y]*G[x, y]*u,
{x, 0, 2 Pi}, {y, 0, 2 Pi}]

(*    1.35767*10^43    *)