Lecture 2

  • while, break, continue
  • cobweb plots
  • digits in double precison float
  • adding a small number to a large number
  • Using the Decimal package to increase precision
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

examples of while, break

In [9]:
x = 0
while(1==1):   # an infinite loop
    x+=1
    #print(x)
    if (x>100):
        break
        
print(x)

101

examples of continue

In [10]:
mysum=0
for i in range(5):
    print(mysum)
    mysum+=1
    if (mysum%2==0):  #my sum is even?
        continue
    mysum += 2   # this line is only reached if mysum is odd
    
        

0
3
4
7
8

maps of the line onto itself

cobweb plots

Cobweb plots are a way of visualizing iterates of a map.

Iterates of a map can go to $\pm \infty$, converge onto a fixed point, converge onto a periodic orbit, or give a chaotic orbit.

In [11]:
# define a function of variable x and parameter mu
def f(x,mu):
    return mu + x**3  #function defined here


# mu is argument to function f
# x0 is initial condition
# nint is the number of iterations
# return 2 vectors suitable for making a cobweb plot
def cobweb(mu,f,x0,nint):
    x = x0
    xvec = np.zeros(0)
    yvec = np.zeros(0)
    for i in range(nint):   #numbers of iteractions 
        y = f(x,mu)
        xvec = np.append(xvec,x) # add point (x,f(x)) 
        yvec = np.append(yvec,y)
        
        xvec = np.append(xvec,y) # add point (f(x),f(x)) on y=x line
        yvec = np.append(yvec,y)
        x = y   #to iterate!
        
    return xvec, yvec   # return coordinates

#############cobweb(mu,f,x0,nint)
xvec, yvec = cobweb(0.43,f,-0.8,10)
print(xvec)
print(yvec)

xv = np.linspace(-1,1.3,101)
yv = f(xv,0.43) #the function with my=0.43

plt.plot(xv,yv,'k-')# plot the function
plt.plot(xv,xv,'b-')  # plot y=x
plt.plot(xvec,yvec,'r-')  # plot the cobweb
plt.plot(xvec[0],yvec[0],'ro') # initial condition
    

[-0.8        -0.082      -0.082       0.42944863  0.42944863  0.50920155
  0.50920155  0.56202894  0.56202894  0.60753175  0.60753175  0.65423683
  0.65423683  0.71003026  0.71003026  0.78795677  0.78795677  0.91922334
  0.91922334  1.20671758]
[-0.082      -0.082       0.42944863  0.42944863  0.50920155  0.50920155
  0.56202894  0.56202894  0.60753175  0.60753175  0.65423683  0.65423683
  0.71003026  0.71003026  0.78795677  0.78795677  0.91922334  0.91922334
  1.20671758  1.20671758]
Out[11]:
[<matplotlib.lines.Line2D at 0x7fbaa09c4e80>]
In [13]:
# showing a period 2 orbit of a different function 

def g(x,mu):
    return mu*np.sin(np.pi*x)

mu=0.81
xvec, yvec = cobweb(mu,g,0.3,30)


xv = np.linspace(0,1.0,101)
yv = g(xv,mu)
plt.plot(xv,yv,'k-')
plt.plot(xv,xv,'b-')
plt.plot(xvec,yvec,'r-')
plt.plot(xvec[0],yvec[0],'ro')
Out[13]:
[<matplotlib.lines.Line2D at 0x7fbaa0ae97c0>]

showing a period 2 orbit

Precision of floating point numbers

How many digits are important?

In [23]:
for i in range(1,21):
    print(i, 1 + 10**-i)

1 1.1
2 1.01
3 1.001
4 1.0001
5 1.00001
6 1.000001
7 1.0000001
8 1.00000001
9 1.000000001
10 1.0000000001
11 1.00000000001
12 1.000000000001
13 1.0000000000001
14 1.00000000000001
15 1.000000000000001
16 1.0
17 1.0
18 1.0
19 1.0
20 1.0

A double precission float devotes a certain number of bits to the number and a certain number of bits to the exponent. The bit format is defined by the IEEE.

The accuracy of double precision numbers is $\sim 10^{-14}$ due to number of bits given to the number itself (not the exponent).

adding something small to something large gives an error

In [37]:
x=1e-13
y=1.0
z1 = x+y
print(z1)

1.0000000000001
In [38]:
x = 1.005e-13
z2 = x+y
print(z2)

1.0000000000001006
In [39]:
print(z2-z1)   # the answer should be 5e-16

6.661338147750939e-16

To improve precission add all the smallest things together first. Then add the larger things.

Comments on precission and Round-off error

Double precision floating point gives you and error of about $10^{-16}$ in the digits because only about 16 digits are stored in double precision format.

That is to say a number like 1 is actually $1 \pm 10^{-16}$.

A number like $10^{-14}$ is actually $ = 10^{-14} (1 \pm 10^{-16}) =10^{-14} \pm 10^{-30}$

Suppose we add $1 + \sqrt{2} \times 10^{-14}$? This is actually $$ 1 \pm 10^{16} + \sqrt{2} \times 10^{-14} ( 1 \pm 10^{-16} ) \sim 1 + \sqrt{2} \times 10^{-14}\pm 10^{16} $$ Let $z$ be this number and now let us subtract 1 from it. $$z - 1 \sim \sqrt{2} \times 10^{-14} \pm 10^{-16}$$

The error as a fraction is now $$\frac{{\rm err} (1-z)}{1-z} \sim \frac{10^{-16}}{10^{-14}} \sim 10^{-2} $$ or we have an error of 1\% which is pretty big.

Precision and using the Decimal package

You may need more precision than double precission floating point.

Your code will run slower if you use higher precision. This is because more computations are necessary. But also because chips dedicated to accelerating computing floating point operations are common and if you use higher precession then the computations will not be done on these accelerators.

In [40]:
from decimal import *
print("digits of precision ",getcontext().prec)

digits of precision  28
In [41]:
getcontext().prec = 28  # you can set the precision to whatever you want
In [42]:
x = Decimal(1.0)
y = Decimal(2.0)
print(x,y)

1 2
In [43]:
y.sqrt()
Out[43]:
Decimal('1.414213562373095048801688724')
In [44]:
z = x + Decimal(1)/Decimal(1000)
print(z)

1.001
In [45]:
small_14 = Decimal((0,(1,),-14)) 
small_14
Out[45]:
Decimal('1E-14')
In [46]:
Decimal.from_float(0.1)
Out[46]:
Decimal('0.1000000000000000055511151231257827021181583404541015625')

Not all arithmetic or numpy operations work on Decimal types.