In[]:=

deploy

Tue 6 Jul 2021 17:13:01

https://math.stackexchange.com/questions/4192183/solve-int-1n-exp-s-i-p-i-p-mathrm-di-zetap-epsilon

Exponential integral

Use approximation

s

(1-

-p

i

)

-p

i

≃exp(-

p

si

)

In[]:=

With[{p=1.1,s=100,n=1000},Print["exact: ",Sum[

s

(1-

-p

i

)

-p

i

,{i,1,n}]];Print["integral1: ",Integrate[

s

(1-

-p

i

)

-p

i

,{i,1,n}]];Print["integral2: ",Integrate[Exp[-s

-p

i

]

-p

i

,{i,1,n}]]]

exact: 1.28971

integral1: 1.28947

integral2: 1.29258

Solve problem in https://math.stackexchange.com/questions/4180259/behavior-of-inverse-of-fs-frac1h-n-sum-i-1n-left1-frac1i-right

In[]:=

solveExact[p_,n_,eps_]:=formula[s_]=

1

HarmonicNumber[n,p]

Sum

s

1-

1

p

i

1

p

i

,{i,1,n};invert[f_]:=FindRoot[f[s]eps,{s,2},WorkingPrecisionMachinePrecision+Log[n]];s/.invert[formula];solveIntegralExp[p_,n_,eps_]:=(*hardcodeoutputofIntegrate[Exp[-s

-p

i

]

-p

i

,{i,1,n}]*)f[s_]=

1-p

n

ExpIntegralE

1

p

,

-p

n

s-ExpIntegralE

1

p

,s

p

Zeta[p];invert[f_]:=FindRoot[f[s]eps,{s,100}];SetPrecision[s/.invert[f],5];

Compare exact and approximate solutions

In[]:=

eps=

-2

10

;n=1000;solveExact[2,n,eps]solveIntegralExp[2,n,eps]

Out[]=

2582.0880492818863065699

Out[]=

2579.7

In[]:=

plot1=With[{p=1.5,eps=

-2

10

},DiscretePlot[solveIntegralExp[p,n,eps],{n,100,100000,1000},FillingAxis,FillingStyleOpacity[.5],AxesLabel{"n","s"},PlotLabel{"ϵ"≐N[eps],"p"≐p}]]

Out[]=

In[]:=

plot2=With[{p=1.8,eps=

-2

10

},DiscretePlot[solveIntegralExp[p,n,eps],{n,100,100000,1000},FillingAxis,FillingStyleOpacity[.5],AxesLabel{"n","s"},PlotLabel{"ϵ"≐N[eps],"p"≐p},PlotRangeAll]]

Out[]=

In[]:=

GraphicsRow[{plot1,plot2}]

Out[]=

