Using `FindRoot` to solve over positive integers?

Question

I need to solve for the smallest positive integer $s$ such that the following is satisfied for $p>1, n\approx \infty, \epsilon>0$

$$\frac{\sum_{i=1}^n \left(1-i^{-p}\right)^si^{-p}}{\sum_{i=1}^n i^{-p}}\le \epsilon$$

Using FindRoot for this problem sometimes fails because it tries negative $s$. There's an 8-year old answer for adding dynamic constraints to FindRoot but looks a bit awkward, is it still the best way to do it?

solve[p_, n_, eps_] := (
   formula[s_] = 
    1/ HarmonicNumber[n, p] Sum[(1 - 1/i^p)^s 1/i^p, {i, 1, n}];
   invert[f_] := 
    FindRoot[f[s] == eps, {s, 2}, AccuracyGoal -> 4, 
     PrecisionGoal -> 4];
   invert[formula]
   );
solve[6, 1000, 10^-2] (* works *)
solve[7, 1000, 10^-2] (* fails *)

Answer 1

Here is one way to force s to be non-negative by replacing s with Exp[logs]:

solve[p_, n_, eps_] := (formula[s_] = 1/HarmonicNumber[n, p] Sum[(1 - 1/i^p)^s 1/i^p, {i, 1, n}];
  invert[f_] := FindRoot[f[Exp[logs]] == eps, {logs, Log[2]}, AccuracyGoal -> 4, PrecisionGoal -> 4];
  invert[formula]);
logs /. solve[6, 1000, 10^-2] // Exp // Ceiling
(* 39 *)
logs /. solve[7, 1000, 10^-2] // Exp // Ceiling
(* 1 *)

Answer 2

The expression is a ratio of two sums, and is very insensitive to $n$. Any $n\ge 10$ gives results equivalent to $n= \infty$. Therefore, redefine the expression for the limiting value of $n\rightarrow\infty$.

The sum in the denominator is the Harmonic number of order $p$, or Zeta[p] when $n\rightarrow\infty$.

The sum in the numerator is a linear combination of Harmonic numbers of order $k*p$, for $k=1,2,...,s+1$. When $n\rightarrow\infty$, the sum becomes a linear combination of Zeta[k*p], for $k=1,2,...,s+1$. The coefficients of the linear combination are (-1)^Range[0,s]*Binomial[s,Range[0,s]].

Define the ratio of the two sums to be fz[s,p] with 20 digits of precision.

fz[s_Integer, p_] :=
    N[(((-1)^Range[0,s] Binomial[s, Range[0,s]]) . Zeta[Range[s + 1]*p])/Zeta[p], 20]

Minimizing over integer $s$ gives a solution in a reasonable time. The value of $s$ returned is the one minimizing the distance of fz from $\epsilon$, and is usually 1 less than the value of $s$ for which fz first becomes $\le\epsilon$.

Minimize[{Abs[fz[s, 6] - 0.005], 100 >= s >= 1}, s, Integers]

{0.0000146715, {s -> 94}}

Table[{s, fz[s, 6]}, {s, 91, 95}]

{{91, 0.0051888315300105947735}, {92, 0.0051298860275994943465}, {93, 0.0050718373295032803372}, {94, 0.0050146714552903004091}, {95, 0.0049583746429293128651}}

A brute force search is usually much faster than Minimize.

search[p_, e_] := Block[{s = 1}, While[fz[s, p] > e, s += 1]; {s, fz[s, p]}]

For example,

search[6, 0.005]

{95, 0.0049583746429293128651}

The value of fz[s=1,p] is larger than any other fz[s>1,p].

The value of fz[s=1,p]$\approx 10^{-3 p / 10}$. If this maximum of fz is less than $\epsilon$, then the minimum positive value of s is $s_{min}=1$. That is, $s_{min}=1$ when $p>\rm{Ceiling}[-10*Log[10,\epsilon]/3]$.

$MaxExtraPrecision = 10^4;
Block[{s},
   ListLogPlot[
      Table[
         s = 1;
         While[fz[s, p] > 0.002, s += 1]; {p, s},
      {p, 4, 10, 1/4}],
      Frame -> True, FrameLabel -> {"Exponent  p", "Minimum Integer  s"},
      PlotRange -> {{4, 10}, Automatic}, ImageSize -> 600, 
      BaseStyle -> {FontSize -> 15},
      FrameTicks -> {Automatic, {1, 2, 5, 10, 20, 50, 100, 200, 500, 1000}},
      PlotLabel -> "\[Epsilon] = 0.002", 
      GridLines -> {Range[4, 12], {1, 2, 5, 10, 20, 50, 100, 200, 500, 1000}}
]]