Problem set #2 PHY256 Problems with Solutions

In [7]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
# try to get a pretty picture for the logistic-like map biffurcation plot
%config InlineBackend.figure_format = 'retina'  # if you have a retina screen

Problem 1. On precision

In [9]:
# a solution 
from scipy.special import zeta

def myzeta(s):
    zsum = 0
    Nmax = 1e7
    # use lots of terms in the sum
    # sum down to try to improve accuracy
    for i in range(int(Nmax),0,-1): 
        #print(i) 
        zsum += i**-s
    return zsum
 
# how accurate is it?
def diff_zeta(x):
    diff = myzeta(x) - zeta(x)
    print(np.abs(diff))

diff_zeta(2)  # satisfies the requested accuracy
diff_zeta(3)

9.99999949513608e-08
4.884981308350689e-15

Convergence is very slow as the terms are close to $1/n$ as $s$ gets near 1 and a sum of $1/n$ does not converge.

Summing small things first is more accurate than summing large things first because of round off error.

I found I need at least $10^7$ terms in the sum to achieve an accuracy of $10^{-7}$ at $s=2$.

In [31]:
# another solution
# let's try comparing fsum to sum on 2 lists
from math import fsum
def myzeta2(s):
    Nmax = 2e7 #  use lots of terms in the sum
    pvec_up = np.arange(1,int(Nmax),dtype=float)  #increasing
    # sum down to try to improve accuracy?
    pvec_down = np.arange(int(Nmax-1),0,-1,dtype=float) #decreasing
    # list of integers (but a floats) starting with Nmax, going to 1
    pvec_up = pvec_up**-s  # only works if array is float
    pvec_down = pvec_down**-s
    
    sumf_up = fsum(pvec_up)  #  try different ways of summing
    sum_up = np.sum(pvec_up)
    sumf_down = fsum(pvec_down) 
    sum_down = np.sum(pvec_down)
    return sumf_up,sum_up,sumf_down,sum_down

def diff_zeta2(x):
    sumf_up,sum_up,sumf_down,sum_down=myzeta2(x)
    zz = zeta(x)
    df_up = np.abs(sumf_up - zz)
    d_up = np.abs(sum_up - zz)
    df_down = np.abs(sumf_down - zz)
    d_down = np.abs(sum_down - zz)
    print(df_up)
    print(d_up)   # this one seems the best!
    print(df_down)
    print(d_down) # then this one
    

diff_zeta2(2)  # satisfies the requested accuracy
#diff_zeta2(3)

5.0000001250438686e-08
4.999999969612645e-08
5.0000001250438686e-08
5.000000102839408e-08

Problem 2. Iterates of a function like the logistic map, Feigenbaum plot

In [3]:
# solution

# like the logistic map but with x^alpha, here alpha=3
def myf(x,mu):
    alpha = 3.0
    y = mu*x**alpha*(1.0-x**alpha)
    if (np.abs(y) > 10):  # a cutoff to keep things from getting too large
        y=0.0
        
    return y

# x0 is initial condition
# j is the number of iterations to discard
# n is the number of iterations to use
# mu is the parameter to pass to function
# returns an array filled with iterations of the function myf
def myit(j,n,x0,mu):
    y= x0  #initial condition
    for i in range(j):
        y = myf(y,mu)   # discard i iterates
    
    output = np.zeros(0);  # array starter
    for i in range(n):   # plot n subsequent iterates
        #plt.plot(mu,y, 'ro', ms=1) too slow
        output = np.append(output,[y])
        y = myf(y,mu)
        
    return output
In [4]:
mu0 = 3.367  # lowest value of mu
mumax=3.75  # highest value of mu
nmu=1000   # number of mus to calculate and plot
dmu = (mumax-mu0)/nmu  # spacing between mu values
nit= 1000  # number of iterations to plot
nitskip =300  # number of iterations to skip before plotting

mu=mu0; # iterating with this as initial mu value
mutot=np.zeros(0)  # trying to append to arrays in case this helps with scatter plot
ztot=np.zeros(0)

# loop over all the different values of mu
for i in range(nmu):
    z = myit(nitskip,nit,0.8,mu) # this is an array of iterated values
    marr = z*0.0 + mu;   # an array of mu values
    mutot = np.append(mutot,marr);  # append mu values to the mu array
    ztot = np.append(ztot,z);  #append iterates to the z array
    mu = mu + dmu; # increase mu and do it again
    
# we now should have two arrays, a ztot array containing iterates
# and an array the same length mutot that has mu values for each iterate
In [50]:
# plot our results
fig,ax = plt.subplots(figsize=(3,3),dpi=300)
plt.rcParams['lines.markersize']=5
#plt.scatter(mutot,ztot,c='b', s=0.1,marker='o', alpha=0.1) #plotting all iterates at once
#ax.ylim=([0.3,1.]) # ignored, backend problem somehow
ax.plot(mutot,ztot,'b,',ms=1,alpha=0.7) #plotting all iterates at once
axes = plt.gca()  # finally I got the overide for the y limits!
axes.set_ylim([0.55,0.95])
#ax.set_ylim=([0.25,1.])  # didn't work!
ax.set_xlabel('$\mu$', fontsize=22);
ax.set_ylabel('1000 iterates after 300')
ax.text(3.38,0.58,r'$f(x) = \mu x^3 (1-x^3)$',fontsize=9)

ax.plot([3.39,3.545],[0.73,0.73],'r-')
ax.annotate('period \ndoubling', (3.4,0.7),fontsize=6)
ax.annotate('period 3', (3.63,0.68),xytext=(3.6,0.62),\
            fontsize=6,arrowprops=dict(arrowstyle='->'))
Out[50]:
Text(3.6, 0.62, 'period 3')
In [ ]: