A365125 Put a positive charge at 0 and a negative charge at 1, then keep adding alternating charges at points of zero potential; this is the decimal expansion of the limit.

6, 8, 7, 8, 4, 1, 8, 1, 0, 3, 2, 8, 3, 8, 9, 2, 6, 3, 2, 7, 1, 3, 4, 4, 0, 4, 4, 0, 9, 8, 8, 3, 3, 4, 8, 6, 1, 1, 5, 8, 3, 9, 7, 9, 4, 8, 7, 6, 6, 8, 9, 5, 4, 1, 1, 7, 4, 7, 5, 8, 6, 6, 9, 4, 4, 1, 0, 7, 8, 5, 2, 8, 1, 7, 2, 1, 2, 4, 7, 5, 3, 8, 9, 1, 0, 8, 7, 9, 1, 2, 6, 5, 7, 8, 9, 7, 8, 5, 3, 6

Offset: 0

Author

Rok Cestnik, Aug 22 2023

Keywords

• nonn
• cons

Comments

The potential is inversely proportional to the distance: a positive charge at p has potential 1/abs(x-p). Starting from a positive charge at 0 and a negative charge at 1, the potential is zero at the midpoint 1/2, so another positive charge is added there. Then the potential is zero at 0.7886... and a negative charge is added there, and so on. The first few zero potential points are: 0.5, 0.7886..., 0.6325..., 0.7179..., 0.6714..., ... The limit is the constant of this sequence.

After initial points p_0 = 0 and p_1 = 1, each new zero point fits between previous two points and is the solution to Sum_{n>=0} 1/(x-p_n) = 0.

Examples

			0.68784181032838926327134404409883348611583979...
		

Links

• Rok Cestnik, Visualization

Programs

Mathematica

p={0,1};
For[r=1,r<39,++r,p=Append[p,x/.NSolve[{Sum[1/(x-p[[k]]),{k,1,Length[p]}]==0,(x-p[[-1]])*(x-p[[-2]])<0},x,WorkingPrecision -> 30][[1,1]]];];
A365125 = RealDigits[Floor[10^10*p[[-1]]]][[1]]

Python

from decimal import Decimal, getcontext, ROUND_UP, ROUND_DOWN
getcontext().prec = 100
def nm(f, df, x):
    for i in range(10):
        x -= f(x)/df(x)
    return x
def flip_rounding():
    if getcontext().rounding == ROUND_UP: getcontext().rounding = ROUND_DOWN
    else: getcontext().rounding = ROUND_UP
def get_zero(vs, rounding):
    getcontext().rounding = rounding
    def p(x,v):
        flip_rounding(); t = x-v
        flip_rounding(); return 1/t
    def dp(x,v):
        flip_rounding(); t = x-v; t = t**2
        flip_rounding(); return -1/t
    def f(x): return sum(p(x,vs[n]) for n in range(len(vs)))
    def df(x): return sum(dp(x,vs[n]) for n in range(len(vs)))
    sign = -1 if rounding == ROUND_DOWN else 1
    return nm(f, df, (vs[-1]+vs[-2])/2+sign*abs(vs[-1]-vs[-2])/3)
v_lo = [Decimal(0), Decimal(1)]
v_up = [Decimal(0), Decimal(1)]
for r in range(150):
    v_lo.append(get_zero(v_lo, ROUND_DOWN))
    v_up.append(get_zero(v_up, ROUND_UP))
lower_bounds = [v_lo[i] for i in range(0, len(v_lo), 2)]
upper_bounds = [v_up[i] for i in range(1, len(v_up), 2)]
right = True
A365125 = [int(l) for l, u in zip(str(lower_bounds[-1])[2:], str(upper_bounds[-1])[2:]) if right and (right := (l == u))]