# NIntegrate a very sharp function

I have a very sharp function s[t] which is I want to numerically integrate, here it is:

u[j_, k_, k1_,
a_] := (SphericalBesselJ[j + 1/2, k a] SphericalBesselJ[j - 1/2,
k1 a] - SphericalBesselJ[j - 1/2, k a] SphericalBesselJ[j + 1/2,
k1 a])/(SphericalHankelH1[j + 1/2, k a] SphericalBesselJ[
j - 1/2, k1 a] -
SphericalHankelH1[j - 1/2, k a] SphericalBesselJ[j + 1/2, k1 a])
lmax = 60;
s[t_, k_, k1_, a_] :=
2/(1 + Cos[t]) Norm[
Sum[(l + 1) (LegendreP[l + 1, Cos[t]] + LegendreP[l, Cos[t]]) u[
l + 1/2, k, k1, a], {l, 0, lmax}]]^2
k = 0.5; q = 1; a = 100; k1 = k - q;
Plot[s[t, k, k1, a], {t, -Pi, Pi}, PlotRange -> All]


The integral I want to do is:

NIntegrate[Sin[t] s[t, k, k1, a], {t, 0, Pi}]/
NIntegrate[Sin[t] (1 - Cos[t]) s[t, k, k1, a], {t, 0, Pi},


I've received a lot of errors when doing this integral and I am not sure if the result is correctly or not.

### Edit:

The possible result is 2.4951, I derived it semi-analytically, but I am not sure if it is right or not.

• Up to the plot, 2.4951 is not correct. Feb 18, 2017 at 6:43
• @user64494 you mean the value should be larger? Feb 18, 2017 at 9:07
• I mean the unclaimed edit in your question. There was no fraction in the original question and the integrand in the numerator was changed too. Compare with the Sumit's answer dropbox.com/s/0agjojzr5etq10z/screen18.02.17.docx?dl=0. Feb 18, 2017 at 17:25
• It seems you have caused great confusion because the the function in the plot shown here is not the integrand. Feb 18, 2017 at 19:02
• @user64494 review the edit history, I have the fraction all the time. Feb 19, 2017 at 2:41

To supplement @buzhidao's answer, there is a way to predict the needed WorkingPrecision for a given PrecisionGoal. (AccuracyGoal doesn't matter much in this example, since the integrals are so large.)

First a measure of the conditioning:

k = 1/2; k1 = k - q; (* to make code exact, for high precision *)
logcond =  With[{wp = 100},
wp - Precision[Sin[t] s[t, k, k1, a] /. t -> N[Subdivide[0, Pi, 16][[2 ;; -2]], wp]]
]

(*  Out[106]= 22.2319  *)


There's a loss of 22+ digits. The working precision (WP) needs to be at least this large, plus whatever is needed to achieve a given precision goal (PG). The heuristic Mathematica uses is WP = 2 * PG. In poorly conditioned problems, WP may need to have a higher ratio. In this case 2*PG is enough, but NIntegrate runs a bit faster with 4*PG or higher, up to some point. (This suggests that the ill-conditioning is showing up in the error estimate.)

Block[{wp, pg},
pg = 8; (* precision goal *)
wp = logcond + 4 pg; (* +2pg  is normal, seems to work; +4pg is faster *)
NIntegrate[Sin[t] s[t, k, k1, a], {t, 0, Pi},
WorkingPrecision -> wp, PrecisionGoal -> pg]/
NIntegrate[Sin[t] (1 - Cos[t]) s[t, k, k1, a], {t, 0, Pi},
WorkingPrecision -> wp, PrecisionGoal -> pg] // AbsoluteTiming
]
(*  Out[114]= {7.50135, 2.49510186134610385390934586152274152298365069672002498}  *)


Using the combination of WorkingPrecision and AccuracyGoal, I am able to eliminate the warnings and get the supposed value, i.e. 2.4951. Here is the code:

NIntegrate[Sin[t] s[t, k, k1, a], {t, 0, Pi}, MaxRecursion -> 20,
WorkingPrecision -> 30, PrecisionGoal -> 6]/
NIntegrate[Sin[t] (1 - Cos[t]) s[t, k, k1, a], {t, 0, Pi},
MaxRecursion -> 20, WorkingPrecision -> 30, PrecisionGoal -> 6]


2.49510186134610350267992744002

• Both the numerator (= 2755.39038618878269174009269118) and the denominator (= 1104.31979907315439966452460777) contradict the plot in the question. Feb 18, 2017 at 10:04
• @user64494 I don't see why, please elaborate it, maybe write an answer about it? Feb 18, 2017 at 10:08

In general, you can try different NIntegrate Integration Strategies. For diverging functions, try to use different Singularity Handling techniques.

(I am going with lmax=30 to save time)

You must be careful about what is causing the sharpness. If it is a divergence then you can specify that in your NIntegrate. For example, if I mention that 0 is a singularity

NIntegrate[s[t, k, k1, a], {t, -Pi, Pi, 0}]


13145.6

and without that

NIntegrate[s[t, k, k1, a], {t, -Pi, Pi}]


26291.1

If you are sure that your function does not have any divergence (which I think is true in this case) try to increase MaxRecursion to see if the result changes.

NIntegrate[s[t, k, k1, a], {t, -Pi, Pi}, MaxRecursion -> 20]


26291.1

NIntegrate[s[t, k, k1, a], {t, -Pi, Pi}, MaxRecursion -> 50]


26291.1

• @ Sumit : Your answer 26291.1 is not in accordance with the plot for lmax= 30. Also s[0, k, k1, a] equals 225195. in this case. Therefore, there is no singularity at $t=0$. Feb 18, 2017 at 7:33
• @user64494, I agree there is no singularity ("If you are sure that your function does not have any divergence (which I think is true in this case)..."). I give here a general prescription so I mentioned both possibilities. About lmax=30, I have to check the calculation again :) Feb 18, 2017 at 10:04