# 1. Getting started¶

In [5]:
print("hello")   # this is a code cell; press Shift-Enter to run the code

hello


$$\frac{3}{4}$$

In [6]:
[x**2 for x in range(1,11)]  # The first 10 square numbers

Out[6]:
[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

# 5. Coding exercises.¶

(a) Iteration

In [10]:
import math
tol = 1e-13  # Iterate until |x_{n+1} - x_{n}| < tol  ...
nmax = 100   # ... or nmax iterations, whichever occurs soonest.
x1 = 1
for i in range(nmax):
x0 = x1
x1 = 1/2*(x0 + 2/x0)
if abs(x1-x0) < tol:
break
print("The iterative method finds sqrt(2) to be: %.12f" % x1)
print("The correct value is                    : %.12f" % math.sqrt(2))

The iterative method finds sqrt(2) to be: 1.414213562373
The correct value is                    : 1.414213562373

In [11]:
import decimal
from decimal import Decimal
decimal.getcontext().prec = 102
one, two = Decimal(1), Decimal(2)
x1 = one
tol = 1e-101
nmax = 100
for i in range(nmax):
x0 = x1
x1 = one/two*(x0 + two/x0)
if abs(x1-x0) < tol:
break
print("The value of sqrt(2) after %i iterations to 100d.p. is:\n" % i, x1)

The value of sqrt(2) after 8 iterations to 100d.p. is:
1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157274


(b) Pisano periods

In [18]:
def fib(nterms, N=0):
seq = [0,1]
for i in range(nterms-2):
f = seq[-1] + seq[-2]
if N > 0:
f = f % N
seq.append(f)
return seq
print(fib(12))
print(fib(12, 4))

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]
[0, 1, 1, 2, 3, 1, 0, 1, 1, 2, 3, 1]

In [30]:
def findperiod(N):
f = fib(10*N, N)
for i in range(2,len(f)):
if (f[i-1], f[i]) == (0,1):
break
return i-1

print("pi(4) = ", findperiod(4))
print("pi(10) = ", findperiod(10))
print("pi(144) = ", findperiod(144))

pi(4) =  6
pi(10) =  60
pi(144) =  24


Further reading: the discrete cat map (V.I. Arnold).

(c) roots of a cubic with the Newton-Raphson method

In [35]:
import cmath
def findroot(z1):
tol = 1e-12
nmax = 20
for i in range(nmax):
z0 = z1
z1 = z0 - (z1**3 - 1)/(3*z1**2)
if abs(z0-z1)<tol:
break
return z1

zs = [2.0, 2.0*1j, -2.0*1j]
for z0 in zs:
zroot = findroot(z0)
print("The root with starting guess (%f,%f) is (%f,%f)" % (z0.real, z0.imag, zroot.real, zroot.imag))

The root with starting guess (2.000000,0.000000) is (1.000000,0.000000)
The root with starting guess (0.000000,2.000000) is (-0.500000,0.866025)
The root with starting guess (-0.000000,-2.000000) is (-0.500000,-0.866025)


The (fractal) basins of attraction of the three roots of unity in the complex plane, under the Newton-Raphson method, are shown on a poster near office G18.