login

Reminder: The OEIS is hiring a new managing editor, and the application deadline is January 26.

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.
1

%I #31 Sep 20 2023 07:46:51

%S 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,

%T 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,

%U 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

%N 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.

%C 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.

%C 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.

%H Rok Cestnik, <a href="/A365125/a365125.png">Visualization</a>

%e 0.68784181032838926327134404409883348611583979...

%t p={0,1};

%t 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]]];];

%t A365125 = RealDigits[Floor[10^10*p[[-1]]]][[1]]

%o (Python)

%o from decimal import Decimal, getcontext, ROUND_UP, ROUND_DOWN

%o getcontext().prec = 100

%o def nm(f, df, x):

%o for i in range(10):

%o x -= f(x)/df(x)

%o return x

%o def flip_rounding():

%o if getcontext().rounding == ROUND_UP: getcontext().rounding = ROUND_DOWN

%o else: getcontext().rounding = ROUND_UP

%o def get_zero(vs, rounding):

%o getcontext().rounding = rounding

%o def p(x,v):

%o flip_rounding(); t = x-v

%o flip_rounding(); return 1/t

%o def dp(x,v):

%o flip_rounding(); t = x-v; t = t**2

%o flip_rounding(); return -1/t

%o def f(x): return sum(p(x,vs[n]) for n in range(len(vs)))

%o def df(x): return sum(dp(x,vs[n]) for n in range(len(vs)))

%o sign = -1 if rounding == ROUND_DOWN else 1

%o return nm(f, df, (vs[-1]+vs[-2])/2+sign*abs(vs[-1]-vs[-2])/3)

%o v_lo = [Decimal(0), Decimal(1)]

%o v_up = [Decimal(0), Decimal(1)]

%o for r in range(150):

%o v_lo.append(get_zero(v_lo, ROUND_DOWN))

%o v_up.append(get_zero(v_up, ROUND_UP))

%o lower_bounds = [v_lo[i] for i in range(0, len(v_lo), 2)]

%o upper_bounds = [v_up[i] for i in range(1, len(v_up), 2)]

%o right = True

%o A365125 = [int(l) for l, u in zip(str(lower_bounds[-1])[2:], str(upper_bounds[-1])[2:]) if right and (right := (l == u))]

%K nonn,cons

%O 0,1

%A _Rok Cestnik_, Aug 22 2023