Why does NSum fail here? - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-08-20T10:49:36Z https://mathematica.stackexchange.com/feeds/question/16269 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://mathematica.stackexchange.com/q/16269 1 Why does NSum fail here? fpghost https://mathematica.stackexchange.com/users/1882 2012-12-13T16:01:40Z 2012-12-13T18:35:11Z <p>I want to evaluate a sum of integrals; each integral has a pole on the real axis and I handle this via the Cauchy Principal Value</p> <p>$f(E)=-\frac{a^2}{2q}\sum^{\infty}_{n=-\infty}\int^{\infty}_{-\infty}dr \exp{(-2iEr/q)}\frac{1}{\left[\sinh^2{(r)}-(\frac{rp}{q}+\frac{npP}{2})^2\right]}$</p> <p>$a,p,q, E$ are constants. I should say first of all the following applies for $n\neq 0$, for zero $n$ I handle differently. For each $n$ and $E$ the denominator of the integrand has two poles on the real axis (amongst others in the complex plane I shall not be concerned with). I will use Mathematica's <code>FindRoot</code> to get the locations of these poles, since the equation whose root we need is transcendental. I will then use <code>NIntegrate</code> with the <code>CauchyPrincipalValue</code> method to evaluate the principle parts of my integral. To get the additional contribution from the $i\pi \times \text{residue}$ at each pole, I first use partial fractions on my integrand to obtain</p> <p>$f(E)=-\frac{a^2}{4q}\sum^{\infty}_{n=-\infty}\int^{\infty}_{-\infty}dr \frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\left(\frac{1}{\left[\sinh{(r)}-(\frac{rp}{q}+\frac{npP}{2})\right]}-\frac{1}{\left[\sinh{(r)}+(\frac{rp}{q}+\frac{npP}{2})\right]}\right).$</p> <p>Now the two poles for each $n,E$ are simple poles and providing we know their location (numerically via <code>FindRoot</code>) the residue is easy to compute for each. For example, considering the first term in the integrand</p> <p>$-\frac{a^2}{4q}\frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\frac{1}{\left[\sinh{(r)}-(\frac{rp}{q}+\frac{npP}{2})\right]}$</p> <p>and say we have numerically found the root to be at $r_0$ for this piece, then</p> <p>$\text{Res}=\lim_{r\to r0}\left[-\frac{a^2}{4q}\frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\frac{(r-r_0)}{\left[\sinh{(r)}-(\frac{rp}{q}+\frac{npP}{2})\right]}\right]$ $\text{Res}=\lim_{r\to r0}\left[-\frac{a^2}{4q}\frac{\exp{(-2iEr/q)}}{\left(\frac{rp}{q}+\frac{npP}{2}\right)}\frac{1}{\left[\cosh{(r)}-\frac{p}{q}\right]}\right]$</p> <p>the second equality following from L'Hopital's rule. Note the zeroes of $\left(\frac{rp}{q}+\frac{npP}{2}\right)$ are not really poles when the full expression is considered together as further above. There will be a similar residue from the other piece, and we must add $i\pi\times$ both of these residues to our PV.</p> <p><strong>Mathematica Implementation</strong></p> <pre><code>(*some constants and definitions*) M = 1; R = 10 M; wp = 40; ag = wp - 8; pg = wp/2; a = 1/(4 M) Sqrt[R/(R - 2 M)]; p = 1/(4 M) Sqrt[(M R)/((R - 2 M) (R - 3 M))]; q = 1/(4 M) Sqrt[R/(R - 3 M)]; P = 2 Pi Sqrt[(R^2 (R - 2 M))/M]; (*find the poles for each n,E*) findNroots[n_?IntegerQ, x0_?NumericQ] := Block[{rM, rP, roots}, rP = r /. FindRoot[Sinh[r] + ((r p)/q + (n p P)/2) == 0, {r, x0}, PrecisionGoal -&gt; pg, WorkingPrecision -&gt; wp][]; rM = r /. FindRoot[Sinh[r] - ((r p)/q + (n p P)/2) == 0, {r, x0}, PrecisionGoal -&gt; pg, WorkingPrecision -&gt; wp][]; roots = {rP, rM}; roots ] (*integrate for given n, E, summing PVs and residues*) integrator[n_ /; n != 0, En_?NumericQ] := Block[{rtM, rtP, lowRoot, highRoot, mid, tmp1, tmp2, tmp3, tmp4, sum}, rtP = findNroots[n, 0][]; rtM = findNroots[n, 0][]; lowRoot = Sort[{rtM, rtP}][]; highRoot = Sort[{rtM, rtP}][]; mid = (highRoot + lowRoot)/2; (*const somewhere in middle of two poles on real line*) tmp1 = -a^2/(2 q) NIntegrate[ Exp[-2 I En r /q] (1/((Sinh[r])^2 - ((r p)/q + (n p P)/2)^2)), {r, -Infinity, lowRoot, mid}, Method -&gt; {"PrincipalValue"}]; tmp3 = -a^2/(2 q) NIntegrate[Exp[-2 I En r /q] (1/((Sinh[r])^2 - ((r p)/q+ (n p P)/2)^2)), {r, mid, highRoot, Infinity}, Method -&gt; {"PrincipalValue"}]; (*add residues too:*) tmp2 = I Pi ( -a^2/(4 q) Exp[-2 I En rtM /q] 1/((rtM p)/q + (n p P)/2) ) Limit[1/( Cosh[r] - p/q), r -&gt; rtM]; tmp4 = I Pi ( -a^2/(4 q) Exp[-2 I En rtP /q] 1/((rtP p)/q + (n p P)/2) ) Limit[-1/( Cosh[r] + p/q), r -&gt; rtP]; sum = tmp1 + tmp2 + tmp3 + tmp4; {sum} ] </code></pre> <p>This seems to work OK, I can do <code>Table[integrator[n,1],{n,1,100}]</code> and Mathematica seems to compute this list just fine. However if I do just<code>NSum[integrator[n,1],{n,1,2}]</code> then Mathematica will give all manner of complaints such as <code>NIntegrate::deorela</code>. Why is this when it is perfectly happy to compute the values of <code>integrator</code> in a list individually? Is there another way I should be performing this sum? </p> <p>A separate issue: I also find that at certain values of $n,E$ for example <code>integrator[1, 24]</code> I also get errors like <code>NIntegrate::deodiv: DoubleExponentialOscillatory returns a finite integral estimate, but the integral might be divergent</code> and I can't even compute the result individually let alone in <code>NSum</code>. I suspect this is because things become too oscillatory at higher values of $E$.</p> https://mathematica.stackexchange.com/questions/16269/-/16280#16280 4 Answer by Jens for Why does NSum fail here? Jens https://mathematica.stackexchange.com/users/245 2012-12-13T18:35:11Z 2012-12-13T18:35:11Z <p>If I replace <code>NSum</code> by <code>Sum</code>, your code works fine. With <code>NSum</code>, you'd have to first find a way to suppress forming the derivate with respect to <code>n</code>, which is part of its evaluation method. </p> <p>But reading your comment, I would suggest doing something slightly different anyway. As you noted, the <code>Table</code> gives no problems, but you would also like to adjust the number of terms in the sum manually by looking at the intermediate results. Then why not do this:</p> <pre><code>t = Table[integrator[n, 1], {n, 1, 100}]; sums = Accumulate[t] </code></pre> <p>This outputs a long list, and each entry <code>sums[[m]]</code> is equal to the sum $\sum_{n=1}^{m} \text{integrator}(n,1)$</p> <p>By the way, you may also want to replace the <code>Limit</code> by <a href="http://reference.wolfram.com/mathematica/ref/Residue.html" rel="nofollow"><code>Residue</code></a> (but that's not related to your actual question).</p>