4
$\begingroup$

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]

enter image description here

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}, 
  Method -> "LocalAdaptive"]

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.

$\endgroup$
6
  • $\begingroup$ Up to the plot, 2.4951 is not correct. $\endgroup$
    – user64494
    Feb 18, 2017 at 6:43
  • $\begingroup$ @user64494 you mean the value should be larger? $\endgroup$ Feb 18, 2017 at 9:07
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Feb 18, 2017 at 17:25
  • $\begingroup$ It seems you have caused great confusion because the the function in the plot shown here is not the integrand. $\endgroup$
    – george2079
    Feb 18, 2017 at 19:02
  • 1
    $\begingroup$ @user64494 review the edit history, I have the fraction all the time. $\endgroup$ Feb 19, 2017 at 2:41

3 Answers 3

4
$\begingroup$

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}  *)
$\endgroup$
3
$\begingroup$

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

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

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

$\endgroup$
2
  • $\begingroup$ @ 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$. $\endgroup$
    – user64494
    Feb 18, 2017 at 7:33
  • $\begingroup$ @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 :) $\endgroup$
    – Sumit
    Feb 18, 2017 at 10:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.