Ensuring appropriate accuracy in integration over infinite sums

Question

I would like to be able to evaluate, numerically, the function $f(n, a, w)$ shown below so that I can:

• integrate $f(n, a, w)$ over the range $[0 \le p \le 1]$ for various positive integer $n < 50$ real $0 < a < 5$ (I'm particularly interested in the range $1/50 < a < 3/2$ );

• graphically explore the relationship between the value of the integral and $w$ for the those same values of $n$ and $a$;

• estimate, within about $0.005$, the value of $w$ that minimizes the integral for the same $n$ and $a$.

I have more experience with Mathematica's symbolic capabilities than with the numerical side . Given the oscillatory nature of $f$, and the infinite sums (which are fortunately moderated by the fact that the negative exponential is going to decay rapidly), I'd appreciate guidance as to the overall approach (in terms of functions and options) to ensure the accuracy I'm looking for.

Function

$$ \begin{align} f(n, a, w) &= \frac{1}{n} \left(w^2+\frac{a^2}{12}+4 \sum _{k=1}^{\infty } \left((-1)^k \exp \left(-\frac{2 k^2 \pi ^2 w^2}{a^2}\right) \left(w^2+\frac{a^2}{4 k^2 \pi ^2}\right) \cos \left(\frac{2 \pi p k}{a}\right)-\frac{a^2}{\pi ^2} S^2\right)\right) \\ &-\left(p+\frac{a}{\pi } S \right)^2-\frac{2 p }{n} \left(p+\frac{a S}{\pi }\right)+p^2 \end{align} $$

where $$ S = \sum _{j=1}^{\infty } \frac{(-1)^j }{j} \: \exp \left(-\frac{2 j^2 \pi ^2 w^2}{a^2}\right) \sin \left(\frac{2 j \pi p}{a}\right) $$

Mathematica code

F[n_, a_, w_] := (1/n) * (w^2 + a^2/12 + 4 * Sum[(-1)^k * Exp[-((2 * k^2 * Pi^2 * w^2)/a^2)] * (w^2 + a^2/(4 * k^2 * Pi^2)) * Cos[(2 * Pi * p * k)/a] - (a^2/Pi^2) * S^2, {k, 1, Infinity}]) - (p + (a/Pi) * S)^2 - ((2 * p)/ n) * (p + (a/Pi) * S) + p^2

with

S = Sum[((-1)^j/j) * Exp[-((2 * j^2 * Pi^2 * w^2)/a^2)] * Sin[(2 * j * Pi * p)/a],{j, 1, Infinity}]

Edited: Following the suggestion of @Bill, I have changed the variable of summation in $S$ to $j$ for the sake of clarity.

Answer 1

To solve this problem, we can use numerical function defined as follows

int[km_?NumericQ, jm_?NumericQ, n_?NumericQ, w_?NumericQ, 
  a_?NumericQ] := 
 Module[{p}, 
  S = Sum[((-1)^j/j)*Exp[-((2*j^2*Pi^2*w^2)/a^2)]*
     Sin[(2*j*Pi*p)/a], {j, 1, jm}];
  F = (1/n)*(w^2 + a^2/12 + 
       4*Sum[(-1)^k*
           Exp[-((2*k^2*Pi^2*w^2)/a^2)]*(w^2 + a^2/(4*k^2*Pi^2))*
           Cos[(2*Pi*p*k)/a] - (a^2/Pi^2)*S^2, {k, 1, km}]) - (p + (a/
          Pi)*S)^2 - ((2*p)/n)*(p + (a/Pi)*S) + p^2; 
  NIntegrate[F, {p, 0, 1}, AccuracyGoal -> 5, PrecisionGoal -> 6, 
   Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
     "SymbolicProcessing" -> 0}]];

Example of usage same as Mariusz

int[10, 10, 10, 2, 2]
(*0.366667*)

Visualization

Quiet@Plot[
  Evaluate@Table[int[10, 10, 10, w, a], {a, {1/50, 1, 3/2}}], {w, -1, 
   1}, PlotLegends -> Table[Row[{"a = ", a}], {a, {1/50, 1, 3/2}}], 
  PlotRange -> All]

As we can see the minimal value int is around w=0, so we use FindMinimum to catch it, for example,

FindMinimum[int[10, 10, 10, w, 1], {w, 0}]

Out[]= {-0.542059, {w -> -5.80511*10^-10}}

Answer 2

For example You can start with:

ClearAll["`*"]; Remove["`*"]

P[n_, a_, w_, M_] := NIntegrate[(1/n)*(w^2 + a^2/12 + 
                     4*Sum[(-1)^k*Exp[-((2*k^2*Pi^2*w^2)/a^2)]*(w^2 + a^2/(4*k^2*Pi^2))*
                     Cos[(2*Pi*p*k)/a] - (a^2/Pi^2)*S^2, {k, 1, M}]) - (p + (a/
                     Pi)*S)^2 - ((2*p)/n)*(p + (a/Pi)*S) + p^2 /. 
                     S -> Sum[((-1)^j/j)*Exp[-((2*j^2*Pi^2*w^2)/a^2)]*
                     Sin[(2*j*Pi*p)/a], {j, 1, M}], {p, 0, 1}, 
                     Method -> "GaussKronrodRule", WorkingPrecision -> 20, PrecisionGoal -> 10];

P[10, 2, 2, 10](*For example: n=10, a=2, w=2, M=10 terms.*)

(*0.366667*)

Plot[P[10, 2, w, 10], {w, -5, 5}, PlotRange -> All] // Quiet

(*Nice plot *)