2
$\begingroup$

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

2 Answers 2

4
$\begingroup$

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

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}}
]]

minimum s versus exponent p

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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