Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Experiments with S(x) = sin(x) + sin(2x)/2 + sin(4x)/4 + ...
"""Experiments with the 'sumsin' function: maximization and plotting
S(x) = sin(x) + sin(2x)/2 + sin(4x)/4 + ...
it appears that S(x) reaches maximum at
xmax = 0.8905302010175791857059461558*7027*
whereas
36pi/127=0.8905302010175791857059461558*9025*
Plots of the function:
* One period of S(x): http://imgur.com/yDcTy9E
* Behavior of S(x) near 36pi/127, coordinates relative to (36pi/127,S(36pi/127)): http://imgur.com/Btkd1kH
* Same plot, zoomed to show that left peak is higher: http://imgur.com/hbpjufM
Language: Python 3
Requires: mpmath
"""
from mpmath import mp
def sumsin(x):
"""Calculate S(x) = sin(x) + sin(2x)/2 + ...."""
pi2 = mp.pi*2
if x < 0 or x > pi2:
raise ValueError("x must be in 0...2pi range, for convenience")
s = mp.mpf(0.0)
k = mp.mpf(1.0)
while k > mp.eps:
s += mp.sin(x) * k
x *= 2
if x > pi2:
x = x-pi2
k /= 2
return s
def _intervals_above_treshold( xs, ys, treshold ):
"""Returns list of (x0,x1) pairs, such that corresponding y is above treshold"""
a = None
xprev = None
for xi, yi in zip(xs, ys):
if a is None: #not yet entered range
if yi >= treshold:
#start range
a = xprev if xprev is not None else xi
else:
if yi <= treshold:
#end range
yield (a, xi)
a = None
xprev = xi
#final: end range if it has started
if a is not None:
yield (a, xi)
def search_intervals( func, a, b, eps, subdivisions=50, percent = 0.5, verbose=True, parallel=True ):
"""Extremely redundant maximizer.
Search for maximum of func(x) when x in [a,b]
With X tolerance eps.
On each step, determine set of intervals, where maximum is searched,
subdivide them, then take only sub-intervals above treshold.
Can use parallel processing.
"""
a,b,eps, percent = map(mp.mpf, (a,b,eps,percent))
if subdivisions < 2: raise ValueError("Subdivisions must be > 2")
if percent <= 0 or percent >=1: raise ValueError("percent must be in (0,1)")
print_ = print if verbose else lambda *x:None
print_ ("Start search fom {a} to {b}, subdivisions: {subdivisions}".format(**locals()))
if parallel:
from multiprocessing import Pool
pool = Pool()
map_ = pool.map
print_ ("Using multiprocessing")
else:
def map_(f,arg): return list(map(f,arg))
#list of intervals, where maximum is searched for.
intervals = [ (a,b) ]
step = 0
while True:
print_ ("Step: {step}, intervals: {nintervals}".format(nintervals=len(intervals), **locals()))
span = intervals[-1][1] - intervals[0][0]
if span <= eps:
break
print_ (" X span:", float(span))
xvalues = sum( (mp.linspace(ai,bi,subdivisions) for (ai, bi) in intervals),
[] )
yvalues = map_(func, xvalues)
ymin = min(yvalues)
ymax = max(yvalues)
print_(" Y span: {ymin} ... {ymax}".format(**locals()))
ytreshold = ymax - (ymax-ymin)*percent
print_ (" Y treshold:", ytreshold)
intervals = list(_intervals_above_treshold( xvalues, yvalues, ytreshold ))
step += 1
return intervals[-1][1], intervals[0][0]
def plot_near_extremum():
"""Make the plot of S(x), deeply zoomed around x=36pi/127"""
from matplotlib import pyplot as pp
N = 1500 #number of points
xcenter = mp.pi*36/127 #central point
dx = mp.mpf('-2e-29')*1.5 #horizontal span to plot
xs = mp.linspace(xcenter - dx, xcenter + dx, N)
ys = [sumsin(xi) for xi in xs]
ymin = min(ys)
ymax = max(ys)
yspan = ymax - ymin
dxs_float = list(map(float, mp.linspace(- dx, dx, N)))
dys_float = [float(yi-ymax) for yi in ys]
pp.plot( dxs_float, dys_float )
pp.grid()
pp.show()
def plot_overview():
"""Plot whole period of S(x)"""
from matplotlib import pyplot as pp
N = 1000
xs = [float(xi) for xi in mp.linspace(0, 2*mp.pi, N)]
ys = [float(sumsin(xi)) for xi in xs]
pp.plot( xs, ys )
pp.grid()
pp.show()
def search_maximum():
"""Search for maximum of S(x)"""
pi36_127 = mp.pi*36/127
print ("S(36pi/127={s0}".format(s0 = sumsin(pi36_127)))
#set disired tolerance 2 digits less than
xtol = mp.eps*100
print ("Searching argmax(S(x)) with tolerance {xtol}".format(xtol=float(xtol)))
a, b = search_intervals( sumsin, 0, mp.pi*2, xtol,
subdivisions=50, #on each step, divide each interval in this number of points.
#more - slower & reliable, less - faster & error-prone
percent = 0.2, #on each step, percent of Y span taken to the next step
verbose = True,
parallel = True)
c=(a+b)/2
print ("Finished search, xmax={c}".format(c=c))
print ("X tolerance = {xtol}".format(xtol=float(b-a)))
print ("xmax - 36pi/237 = {diff}".format(diff=float(c-pi36_127)))
print ("S(xmax) - S(36pi/127) = {df}".format( df = float(sumsin(c) - sumsin(pi36_127))))
if __name__ == "__main__":
#more then 31 to get interesting result
mp.dps=40
print ("Set precision to {d} digits".format(d=mp.dps))
#plot_all()
#plot_near_extremum()
search_maximum()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.