I am trying to integrate the following function,

f(r)=(80.50305990274495 - 35.03089630981622 r)/(80.50305990274495 - 79.5030598897815 r - r^2) where r =[0,1]

The correct answer could be obtained by using 'NIntegrate' command in Mathematica and it gave 13.027. The 'NIntegrate' command always gave me 13.027 no matter what method I specify to integrate the integrand.

My main aim is to use Gauss Quadrature or other quadrature rules to get the correct result. Mathematica has a built-in command 'GaussianQuadratureWeights' which generates the sampling points and the weights in an interval. Using following code I tried to integrate the function, f(r) using Gauss quadrature rule but i am not getting the correct result and the Reason is presence of singularity at one of the end points (r=1) due to which the convergence is very slow.

 nof=1500;
f= (80.50305990274495 - 35.03089630981622 r)/(80.50305990274495 -
79.50305988978158 r - 1. r^2) ;
xi = gq[[All, 1]];
wi = gq[[All, 2]];
sol = Sum[wi[[y]]*(f /. {r -> xi[[y]]}), {y, 1, nof}];


Ans: 9.24379 (which is not correct) Can anyone help me with this? Thanks ashu

• How can I prove that ? Is there a way to calculate this using any other quadrature rules? – ashu sharma Jan 9 '17 at 15:23
• Your integrand has a value of 3.50774*10^9 at r -> 1, which contributes to the slow convergence of a straightforward rule application. NIntegrate overcomes this by subdividing the interval several times near r == 1.. You seem to appreciate this, but the implication is that the sampling has to be very high near r == 1. You might also need to increase the working precision to avoid roundoff error. – Michael E2 Jan 9 '17 at 16:27
• Be aware the result is extremely sensitive to the difference between 80.50305990274495 and (1+79.50305988978158), which is ~10^-8. You are really pushing things trying to do this with machine precision. – george2079 Jan 9 '17 at 16:29
• @MichaelE2 The slow convergence is due to singularity and NIntegrate is applying either some transformation or some other quadrature rule to avoid this singularity. I didn't specify the precision while using NIntegrate so I do not think that rounding off error would play major role. – ashu sharma Jan 9 '17 at 16:43
• (1) If you mean that the singularity is a little outside the interval, then, yes, that is what I was talking about. You have to expect slow convergence in that case. The further outside the singularity (relatively), the better the convergence, which is why in part recursive subdivision succeeds. (2) If you mean that the true singularity is at r == 1 (exactly) and not a little outside the interval, then you cannot assign a numeric value to the integral. It is a simple pole at the end point of the interval, and AFAIK, no regularization or principal value can be applied to yield a value. – Michael E2 Jan 9 '17 at 17:12

note this has an analytic solution:

res = Integrate[(a - b r)/(a - (a - 1 - eps)  r - r^2), {r, 0, 1},
Assumptions -> {a > 0, b > 0, eps > 0}]


 res/.{a -> 80.50305990274495,
b -> 35.03089630981622,
eps -> 1.2963454310011002*^-8}


13.027

Note the integrand is not singular, having a finite value (a-b)/eps (~10^9) at r=1, it is however not approximately polynomial which is I guess why Gauss integration doesn't work so well.

Edit:

you may find this useful to see how NIntegrate does its adaptive sampling:

rules = {a -> 80.50305990274495, b -> 35.03089630981622,
eps -> 1.2963454310011002*^-8}
xx = Reap[
NIntegrate[
Evaluate[(a - b r)/(a - (a - 1 - eps)  r - r^2) /. rules], {r, 0,
1}, EvaluationMonitor :> Sow[r]]][[2, 1]];
Histogram[xx, {0, 1, .01}]


you see the vast majority of the sample points are taken near the (almost) singularity.

• I appreciate your effort. But I want to get the same result using quadrature rules because I want a polynomial expression in terms of a, b and eps. May be using your methodology i could get something. Let me try. thanks – ashu sharma Jan 9 '17 at 16:56
• your quadrature scheme needs to sample points spaced on the order 10^-8 near r=1. Any more or less uniform sampling scheme (like Gauss) with then need on the order 10^8 sample points... (Not to mention I expect numerical precision errors if you try to go that route) – george2079 Jan 9 '17 at 17:11
• Gauss quadrature works if I divide the interval into three sub-intervals (0 to 0.999, 0.999 to 0.99999 and 0.99999 to 1.0) I used 1000 nodes in each sub-interval except the first one where I used only 500 nodes. Is there a way to find out the number of nodes required for the convergence? – ashu sharma Jan 9 '17 at 18:55

The function f goes to infinity when it gets close to r=1. So you will have some error integrating it numericaly.

f = (80.5031 - 35.0309 r)/(80.5031 - 79.5031 r - 1. r^2);
Plot[f, {r, 0, 0.99}]
Limit[f, r -> 1]
NIntegrate[f, {r, 0, 0.99999999999947}]


What you can do to improve your solution is concetrate the fire int the most nonlinear region of the function. Here is the code:

<< NumericalDifferentialEquationAnalysis

• Here NIntegrate[f, {r, 0,1.}] gives 54.774, and fails to converge. How do you know the rigth answer is 13.027? – Diogo Jan 9 '17 at 16:06