5mm.

Sequences

"Sequences" is a central topic in mathematics:

Sequences and Difference Equations
(Appendix A) x0 , x1 , x2 , . . . , xn , . . . ,

Hans Petter Langtangen Example: all odd numbers

1, 3, 5, 7, . . . , 2n + 1, . . .
Simula Research Laboratory
University of Oslo, Dept. of Informatics For this sequence we have a formula for the n-th term:

xn = 2n + 1

and we can write the sequence more compactly as

(xn )∞
n=0 , xn = 2n + 1

Sequences and Difference Equations (Appendix A) – p.1/?? Sequences and Difference Equations (Appendix A) – p.2/??

Other examples of sequences Finite and infinite sequences

2 Infinite sequences have an infinite number of terms (n → ∞)

1, 4, 9, 16, 25, . . . (xn )∞
n=0 , xn = n
In mathematics, infinite sequences are widely used
1 1 1 1
1, , , , ... (xn )∞
n=0 , xn = In real-life applications, sequences are usually finite: (xn )N
2 3 4 n+1 n=0

Example: number of approved exercises every week in INF1100

1, 1, 2, 6, 24, . . . (xn )∞
n=0 , xn = n!
n x0 , x1 , x2 , . . . , x15
1 1 1 X xj
1, 1 + x, 1 + x + x2 , 1 + x + x2 + x3 , . . . (xn )∞
n=0 , xn =
2 2 6 j=0
j! Example: the annual value of a loan

x0 , x1 , . . . , x20

Sequences and Difference Equations (Appendix A) – p.3/?? Sequences and Difference Equations (Appendix A) – p.4/??

Difference equations Modeling interest rates

For sequences occuring in modeling of real-world phenomena, Put x0 money in a bank at year 0. What is the value after N years
there is seldom a formula for the n-th term if the interest rate is p percent per year?
However, we can often set up one or more equations governing The fundamental information relates the value at year n, xn , to the
the sequence value of the previous year, xn−1 :
Such equations are called difference equations p
xn = xn−1 + xn−1
With a computer it is then very easy to generate the sequence by 100
solving the difference equations
Solution by simulation:
Difference equations have lots of applications and are very easy
to solve on a computer, but often complicated or impossible to start with x0 (known)
solve for xn (as a formula) by pen and paper! compute x1 = x0 + 0.01px0 = (1 + 0.01p)x0
The programs require only loops and arrays compute x2 = x1 + 0.01px1 = (1 + 0.01p)2 x0
...continue this boring procedure...
find that xn = (1 + 0.01p)n x0
...which is what you learned in high school
That is: we can solve this difference equation by hand
Sequences and Difference Equations (Appendix A) – p.5/?? Sequences and Difference Equations (Appendix A) – p.6/??

Simulating the difference equation for interest rates Note about the use of arrays

Simulate = solve math equations by repeating a simple procedure We store the xn values in a NumPy array
(relation) many times (boring, but well suited for a computer) To compute xn , we only need one previous value, xn−1
Let us make a program for Thus, for the computations we do not need to store all the
p previous values, i.e., we do not need any array, just two numbers:
xn = xn−1 + xn−1 x_old = x0
100 for n in index_set[1:]:
x_new = x_old + (p/100.)*x_old
from scitools.std import * x_old = x_new # x_new becomes x_old at next step
x0 = 100 # initial amount
p = 5 # interest rate However, programming with an array x[n] is simpler, safer, and
N = 4 # number of years
index_set = range(N+1)
enables plotting the sequence, so we will continue to use arrays in
x = zeros(len(index_set)) the examples
# solution:
x[0] = x0
for n in index_set[1:]:
x[n] = x[n-1] + (p/100.0)*x[n-1]
print x
plot(index_set, x, ’ro’, xlabel=’years’, ylabel=’amount’)

Sequences and Difference Equations (Appendix A) – p.7/?? Sequences and Difference Equations (Appendix A) – p.8/??
 Daily interest rate Program for daily interest rate

A more relevant model is to add the interest every day from scitools.std import *
x0 = 100 # initial amount
The interest rate per day is r = p/D if p is the annual interest rate p = 5 # annual interest rate
r = p/360.0 # daily interest rate
and D is the number of days in a year import datetime
date1 = datetime.date(2007, 8, 3)
A common model in business applies D = 360, but n counts exact date2 = datetime.date(2011, 8, 3)
(all) days diff = date2 - date1
N = diff.days
New model: index_set = range(N+1)
r x = zeros(len(index_set))
xn = xn−1 + xn−1
100 # solution:
x[0] = x0
How can we find the number of days between two dates? for n in index_set[1:]:
>>> import datetime x[n] = x[n-1] + (r/100.0)*x[n-1]
>>> date1 = datetime.date(2007, 8, 3) # Aug 3, 2007 print x
>>> date2 = datetime.date(2008, 8, 4) # Aug 4, 2008 plot(index_set, x, ’ro’, xlabel=’days’, ylabel=’amount’)
>>> diff = date2 - date1
>>> print diff.days
367

Sequences and Difference Equations (Appendix A) – p.9/?? Sequences and Difference Equations (Appendix A) – p.10/??

But the annual interest rate may change quite often... Payback of a loan

This is problematic when computing by hand A loan L is paid back with a fixed amount L/N every month over
In the program, a varying p is easy to deal with N months + the interest rate of the loan

Just fill an array p with correct annual interest rate for day no. n, Let p be the annual interest rate and p/12 the monthly rate
n=0,...,N (this can be a bit challenging) Let xn be the value of the loan at the end of month n
Modified program: The fundamental relation from one month to the other:
p = zeros(len(index_set))
# fill p[n] for n in index_set p p L
xn = xn−1 + xn−1 − ( xn−1 + )
r = p/360.0 # daily interest rate 12 · 100 12 · 100 N
x = zeros(len(index_set))

x[0] = x0
which simplifies to
for n in index_set[1:]: L
x[n] = x[n-1] + (r[n-1]/100.0)*x[n-1] xn = xn−1 −
N
(The constant term L/N makes the equation nonhomogeneous,
while the previous interest rate equation was homogeneous (all
terms contain xn or xn−1 ))
The program is left as an exercise
Sequences and Difference Equations (Appendix A) – p.11/?? Sequences and Difference Equations (Appendix A) – p.12/??

How to make a living from a fortune (part 1) How to make a living from a fortune (part 2)

We have a fortune F invested with an annual interest rate of p Assume I percent inflation per year and that c0 is q percent of the
percent interest the first year
Every year we plan to consume an amount cn (n counts years) cn then develops as money with interest rate I, and xn develops
Let xn be the development of our fortune with rate p but with a loss cn every year:

A fundamental relation from one year to the other is p pq

xn = xn−1 + xn−1 − cn−1 , x0 = F, c0 = F
100 104
p
xn = xn−1 + xn−1 − cn I
100 cn = cn−1 + cn−1
100
Simplest possibility: keep cn constant
This is a coupled system of two difference equations
Drawback: inflation demands cn to increase...
The programming is still simple: we update two arrays (x[n],
c[n]) inside the loop

Sequences and Difference Equations (Appendix A) – p.13/?? Sequences and Difference Equations (Appendix A) – p.14/??

Fibonacci numbers; mathematics Fibonacci numbers; program

No programming or math course is complete without an example Program:

on Fibonacci numbers! N = int(sys.argv[1])
from numpy import zeros
Fibonacci derived the sequence by modeling rat populations, but x = zeros(N+1, int)
x[0] = 1
the sequence of numbers has a range of peculiar mathematical x[1] = 1
properties and has therefore attracted much attention from for n in range(2, N+1):
x[n] = x[n-1] + x[n-2]
mathematicians print n, x[n]
The difference equation reads

xn = xn−1 + xn−2 , x0 = 1, x1 = 1

This is a homogeneous difference equation of second order (three

levels: n, n − 1, n − 2) – this classification is important for
mathematical solution technique, but not for simulation in a
program

Sequences and Difference Equations (Appendix A) – p.15/?? Sequences and Difference Equations (Appendix A) – p.16/??
 Fibonacci numbers and overflow (part 1) Fibonacci numbers and overflow (part 2)

Run the program with N = 50: Best: use Python scalars of type int – these automatically
2 2 changes to long when overflow in int
3 3
4 5 The long type in Python has arbitrarily many digits (as many as
5 8
6 13 required in a computation)
...
45 1836311903 Note: long for arrays is 64 bit integer (int64), while scalar
Warning: overflow encountered in long_scalars long in Python is an integer with as "infinitely" many digits
46 -1323752223

Can change int to long or int64 for array elements - now we

can generate numbers up to N = 91 before we get overflow and
garbage
Can use float96 despite the fact that xn are integers (float
gives only approximatively correct numbers) – now N up to 23600
is possible

Sequences and Difference Equations (Appendix A) – p.17/?? Sequences and Difference Equations (Appendix A) – p.18/??

Program with Python’s long type for integers Exponential growth with limited resources

The program now avoids arrays and makes use of three int The model for growth of money in a bank has a solution of the type
objects (which automatically changes to long when needed):
import sys xn = x0 C n (= x0 en ln C )
N = int(sys.argv[1])
xnm1 = 1 # "x_n minus 1"
xnm2 = 1 # "x_n minus 2" This is exponential growth in time (n)
n = 2
while n <= N: Populations of humans, animals, and cells also exhibit the same
xn = xnm1 + xnm2
print ’x_%d = %d’ % (n, xn) type of growth as long as there are unlimited resources (space
xnm2 = xnm1 and food)
xnm1 = xn
n += 1 The environment can only support a maximum number M of
Run with N = 200: individuals
x_2 = 2 How can we model this?
x_3 = 3
... We shall introduce a logistic model
x_198 = 173402521172797813159685037284371942044301
x_199 = 280571172992510140037611932413038677189525
x_200 = 453973694165307953197296969697410619233826

Can simulate until the computer’s memory becomes too small to

store all the digits...
Sequences and Difference Equations (Appendix A) – p.19/?? Sequences and Difference Equations (Appendix A) – p.20/??

Modeling logistic growth The evolution of logistic growth

Initially, when there are enough resources, the growth is In a program it is easy to introduce logistic instead of exponential
exponential: growth, just replace
r x[n] = x[n-1] + p/100.0)*x[n-1]
xn = xn−1 + xn−1
100 by
x[n] = x[n-1] + (rho/100.0)*x[n-1]*(1 - x[n-1]/float(M))
The growth rate r must decay to zero as xn approaches M
A very simple r(n) function with this behavior is 500

450
 xn 
r(n) = ̺ 1 − 400

M 350
no of individuals

Observe that r(n) ≈ ̺ for small n when xn << M , and r(n) → 0 300

as xn → M and n is big 250

200
The model for limited growth, called logistic growth, is then
150

̺  xn−1  100
xn = xn−1 + xn−1 1 − 0 50 100 150 200

100 M time units

Sequences and Difference Equations (Appendix A) – p.21/?? Sequences and Difference Equations (Appendix A) – p.22/??

The factorial as a difference equation Taylor series as difference equations

The factorial n! is defined as n(n − 1)(n − 2) · · · 1 (0! = 1) The Taylor series for ex reads
The following difference equation has xn = n! as solution and can ∞
be used to compute the factorial:
X xn
ex =
n=0
n!
xn = nxn−1 , x0 = 1
We can formulate this series as two coupled difference equations
(and solving these difference equations is (probably) the most
efficient way to compute the Taylor series!):
x
an = an−1 , a0 = 1
n
en = en−1 + an , e0 = 1

See the book for how to solve the difference equations by hand
and show that the solution is the Taylor series for ex

Sequences and Difference Equations (Appendix A) – p.23/?? Sequences and Difference Equations (Appendix A) – p.24/??
 Newton’s method for finding zeros A program for Newton’s method

Newton’s method for solving f (x) = 0 reads A simple program can be

def Newton(f, x, dfdx, epsilon=1.0E-7, N=100):
n = 0
f (xn−1 ) while abs(f(x)) > epsilon and n <= N:
xn = xn−1 − , x0 given
f ′ (xn−1 ) x = x - f(x)/dfdx(x)
n += 1
return x, n, f(x)
This is a (nonlinear!) difference equation
Note: f (x) is evaluated twice in each pass of the loop – only one
As n → ∞, we hope that xn → xs , where xs solves f (xs ) = 0 evaluation is strictly necessary (can store the value in a variable
Now we will not simulate N steps, because we do not know how and reuse it)
large N must be in order to have xn as close to the exact solution Note: f(x)/dfdx(x) can give integer division
xs as we want
Note: it could be handy to have an option for storing the x and
The program is therefore a bit different: we simulate the difference f(x) values in each iteration (for plotting or printing a
equation as long as f (x) > ǫ, where ǫ is small convergence table)
However, Newton’s method may (easily) diverge, so to avoid
simulating forever, we stop when n > N

Sequences and Difference Equations (Appendix A) – p.25/?? Sequences and Difference Equations (Appendix A) – p.26/??

A better program for Newton’s method Application of Newton’s method

2
Only one f (x) call in each iteration, optional storage of (x, f (x)) Example: solve e−0.1x sin( π2 x) = 0
values during the iterations, and float division: from math import sin, cos, exp, pi
import sys
def Newton(f, x, dfdx, epsilon=1.0E-7, N=100,
store=False): def g(x):
f_value = f(x) return exp(-0.1*x**2)*sin(pi/2*x)
n = 0
if store: info = [(x, f_value)] def dg(x):
while abs(f_value) > epsilon and n <= N: return -2*0.1*x*exp(-0.1*x**2)*sin(pi/2*x) + \
x = x - float(f_value)/dfdx(x) pi/2*exp(-0.1*x**2)*cos(pi/2*x)
n += 1
f_value = f(x) x0 = float(sys.argv[1])
if store: info.append((x, f_value)) x, info = Newton(g, x0, dg, store=True)
if store: print ’zero:’, x
return x, info for i in range(len(info)):
else: print ’Iteration %3d: f(%g)=%g’ % \
return x, n, f_value (i, info[i][0], info[i][1])

Sequences and Difference Equations (Appendix A) – p.27/?? Sequences and Difference Equations (Appendix A) – p.28/??

Results from this test problem Programming with sound

Start value 1.7: Sound on a computer = sequence of numbers

zero: 1.999999999768449
Iteration 0: f(1.7)=0.340044 Example: A 440 Hz
Iteration 1: f(1.99215)=0.00828786
Iteration 2: f(1.99998)=2.53347e-05 This tone is a sine wave with frequency 440 Hz:
Iteration 3: f(2)=2.43808e-10
This works fine! s(t) = A sin (2πf t) , f = 440
Start value 3:
zero: 42.49723316011362 On a computer we represent s(t) by a discrete set of points on the
Iteration 0: f(3)=-0.40657 function curve (exactly as we do when we plot s(t))
Iteration 1: f(4.66667)=0.0981146
Iteration 2: f(42.4972)=-2.59037e-79 How many points (samples) along the curve do we need to use?
WHAT???? CD quality needs 44100 samples per second
Lesson learned: Newton’s method may work fine or give wrong
results! You need to understand the method to interpret the
results!

Sequences and Difference Equations (Appendix A) – p.29/?? Sequences and Difference Equations (Appendix A) – p.30/??

Making a sound file with single tone (part 1) Making a sound file with single tone (part 2)

r: sampling rate (samples per second, default 44100) We have data as an array with float and unit amplitude
f : frequency of the tone Sound data in a file should have 2-byte integers (int16) as data
m: duration of the tone (seconds) elements and amplitudes up to 215 − 1 (max value for int16
data)
Sampled sine function for this tone: data = note(440, 2)
data = data.astype(numpy.int16)
 n max_amplitude = 2**15 - 1
sn = A sin 2πf , n = 0, 1, . . . , m · r data = max_amplitude*data
r import scitools.sound
scitools.sound.write(data, ’Atone.wav’)
Code (we use descriptive names: frequency=f , length=m, scitools.sound.play(’Atone.wav’)

amplitude=A, sample_rate=r):
import numpy
def note(frequency, length, amplitude=1,
sample_rate=44100):
time_points = numpy.linspace(0, length,
length*sample_rate)
data = numpy.sin(2*numpy.pi*frequency*time_points)
data = amplitude*data
return data

Sequences and Difference Equations (Appendix A) – p.31/?? Sequences and Difference Equations (Appendix A) – p.32/??
 Reading sound from file Playing many notes

Let us read a sound file and add echo Each note is an array of samples from a sine with a frequency
Sound = array s[n] corresponding to the note

Echo means to add a delay of the sound Assume we have several note arrays data1, data2, ...:
# put data1, data2, ... after each other in a new array:
# echo: e[n] = beta*s[n] + (1-beta)*s[n-b]
data = numpy.concatenate((data1, data2, data3, ...))
def add_echo(data, beta=0.8, delay=0.002, The start of "Nothing Else Matters" (Metallica):
sample_rate=44100):
newdata = data.copy() E1 = note(164.81, .5)
shift = int(delay*sample_rate) # b (math symbol) G = note(392, .5)
for i in xrange(shift, len(data)): B = note(493.88, .5)
newdata[i] = beta*data[i] + (1-beta)*data[i-shift] E2 = note(659.26, .5)
return newdata intro = numpy.concatenate((E1, G, B, E2, B, G))
...
Load data, add echo and play: song = numpy.concatenate((intro, intro, ...))
scitools.sound.play(song)
data = scitools.sound.read(filename)
scitools.sound.write(song, ’tmp.wav’)
data = data.astype(float)
data = add_echo(data, beta=0.6)
data = data.astype(int16)
scitools.sound.play(data)

Sequences and Difference Equations (Appendix A) – p.33/?? Sequences and Difference Equations (Appendix A) – p.34/??

Summary of difference equations Summarizing example: music of sequences

Sequence: x0 , x1 , x2 , . . ., xn , . . ., xN Given a x0 , x1 , x2 , . . ., xn , . . ., xN
Difference equation: relation between xn , xn−1 and maybe xn−2 Can we listen to this sequence as "music"?
(or more terms in the "past") + known start value x0 (and more Yes, we just transform the xn values to suitable frequencies and
values x1 , ... if more levels enter the equation) use the functions in scitools.sound to generate tones
Solution of difference equations by simulation: We will study two sequences:
index_set = <array of n-values: 0, 1, ..., N>
x = zeros(N+1)
x[0] = x0 xn = e−4n/N sin(8πn/N )
for n in index_set[1:]:
x[n] = <formula involving x[n-1]>
and
Can have (simple) systems of difference equations:
for n in index_set[1:]:
xn = xn−1 + qxn−1 (1 − xn−1 ) , x = x0
x[n] = <formula involving x[n-1]>
y[n] = <formula involving y[n-1] and x[n]> The first has values in [−1, 1], the other from x0 = 0.01 up to
around 1
Taylor series and numerical methods such as Newton’s method
can be formulated as difference equations, often resulting in a Transformation from "unit" xn to frequencies:
good way of programming the formulas
yn = 440 + 200xn

Sequences and Difference Equations (Appendix A) – p.35/??

(first sequence then gives tones between 240 Sequences
Hz and 640 Hz)
and Difference Equations (Appendix A) – p.36/??

Three functions: two for generating sequences, one for the sound

Module file: soundeq.py

Look at files/soundeq.py for complete code. Try it out in these

examples:
Unix/DOS> python soundseq.py oscillations 40
Unix/DOS> python soundseq.py logistic 100

Try to change the frequency range from 200 to 400.

Sequences and Difference Equations (Appendix A) – p.37/??