App. A: Sequences and difference equations

Hans Petter Langtangen [1, 2]

[1] Simula Research Laboratory
[2] University of Oslo, Dept. of Informatics

Aug 21, 2016










Sequences

Sequences is a central topic in mathematics: $$ x_0,\ x_1,\ x_2,\ \ldots,\ x_n,\ldots, $$

Example: all odd numbers $$ 1, 3, 5, 7, \ldots, 2n+1,\ldots $$

For this sequence we have a formula for the \( n \)-th term: $$ x_n = 2n+1 $$ and we can write the sequence more compactly as $$ (x_n)_{n=0}^\infty,\quad x_n = 2n+1$$









Other examples of sequences

$$ 1,\ 4,\ 9,\ 16,\ 25,\ \ldots\quad (x_n)_{n=0}^\infty,\ x_n=n^2$$ $$ 1,\ {1\over 2},\ {1\over3},\ {1\over4},\ \ldots\quad (x_n)_{n=0}^\infty,\ x_n={1\over {n+1}}$$ $$ 1,\ 1,\ 2,\ 6,\ 24,\ \ldots\quad (x_n)_{n=0}^\infty,\ x_n=n!$$ $$ 1,\ 1+x,\ 1+x+{1\over2}x^2,\ 1+x+{1\over2}x^2+{1\over6}x^3,\ \ldots\quad (x_n)_{n=0}^\infty,\ x_n=\sum_{j=0}^n {x^j\over j!}$$









Finite and infinite sequences









Difference equations









Modeling interest rates

Problem:

Put \( x_0 \) money in a bank at year 0. What is the value after \( N \) years if the interest rate is \( p \) percent per year?

Solution:

The fundamental information relates the value at year \( n \), \( x_n \), to the value of the previous year, \( x_{n-1} \): $$ x_{n} = x_{n-1} + {p\over 100}x_{n-1} $$

How to solve for \( x_n \)? Start with \( x_0 \), compute \( x_1 \), \( x_2 \), ...









Simulating the difference equation for interest rates

What does it mean to simulate?

Solve math equations by repeating a simple procedure (relation) many times (boring, but well suited for a computer!)

Program for \( x_{n} = x_{n-1} + (p/100)x_{n-1} \):

from scitools.std import *
x0 = 100                      # initial amount
p = 5                         # interest rate
N = 4                         # number of years
index_set = range(N+1)
x = zeros(len(index_set))

# 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')









We do not need to store the entire sequence, but it is convenient for programming and later plotting

Thus, we could only store the two last values in memory:

x_old = x0
for n in index_set[1:]:
    x_new = x_old + (p/100.)*x_old
    x_old = x_new  # x_new becomes x_old at next step

However, programming with an array x[n] is simpler, safer, and enables plotting the sequence, so we will continue to use arrays in the examples









Daily interest rate

Just a minor change in the model: $$x_{n} = x_{n-1} + {r\over 100}x_{n-1}$$

How can we find the number of days between two dates?

>>> import datetime
>>> date1 = datetime.date(2007, 8, 3)  # Aug 3, 2007
>>> date2 = datetime.date(2008, 8, 4)  # Aug 4, 2008
>>> diff = date2 - date1
>>> print diff.days
367









Program for daily interest rate

from scitools.std import *
x0 = 100                           # initial amount
p = 5                              # annual interest rate
r = p/360.0                        # daily interest rate
import datetime
date1 = datetime.date(2007, 8, 3)
date2 = datetime.date(2011, 8, 3)
diff = date2 - date1
N = diff.days
index_set = range(N+1)
x = zeros(len(index_set))

# Solution:
x[0] = x0
for n in index_set[1:]:
    x[n] = x[n-1] + (r/100.0)*x[n-1]
print x
plot(index_set, x, 'ro', xlabel='days', ylabel='amount')









But the annual interest rate may change quite often...

Varying \( p \) means \( p_n \):

Modified program:

p = zeros(len(index_set))
# fill p[n] for n in index_set (might be non-trivial...)

r = p/360.0                      # daily interest rate
x = zeros(len(index_set))

x[0] = x0
for n in index_set[1:]:
    x[n] = x[n-1] + (r[n-1]/100.0)*x[n-1]









Payback of a loan

The fundamental relation from one month to the text: $$ x_n = x_{n-1} + {p\over 12\cdot 100}x_{n-1} - ({p\over 12\cdot 100}x_{n-1} + {L\over N})$$ which simplifies to $$ x_n = x_{n-1} - {L\over N}$$

(\( L/N \) makes the equation nonhomogeneous)









How to make a living from a fortune with constant consumption

A fundamental relation from one year to the other is $$ x_n = x_{n-1} + {p\over 100}x_{n-1} - c_n $$

Simplest possibility: keep \( c_n \) constant, but inflation demands \( c_n \) to increase...









How to make a living from a fortune with inflation-adjusted consumption

\( x_n \) develops with rate \( p \) but with a loss \( c_n \) every year: $$ \begin{align*} x_n &= x_{n-1} + {p\over 100}x_{n-1} - c_{n-1}, \quad x_0=F,\ c_0 = {pq\over 10^4}F\\ c_n &= c_{n-1} + {I\over100}c_{n-1} \end{align*} $$

This is a coupled system of two difference equations, but the programming is still simple: we update two arrays, first x[n], then c[n], inside the loop (good exercise!)









The mathematics of Fibonacci numbers

No programming or math course is complete without an example on Fibonacci numbers: $$ x_n = x_{n-1} + x_{n-2},\quad x_0=1,\ x_1=1$$

Mathematical classification.

This is a homogeneous difference equation of second order (second order means three levels: \( n \), \( n-1 \), \( n-2 \)). This classification is important for mathematical solution technique, but not for simulation in a program.

Fibonacci derived the sequence by modeling rat populations, but the sequence of numbers has a range of peculiar mathematical properties and has therefore attracted much attention from mathematicians.









Program for generating Fibonacci numbers

N = int(sys.argv[1])
from numpy import zeros
x = zeros(N+1, int)
x[0] = 1
x[1] = 1
for n in range(2, N+1):
    x[n] = x[n-1] + x[n-2]
    print n, x[n]









Fibonacci numbers can cause overflow in NumPy arrays

Run the program with \( N=50 \):

2 2
3 3
4 5
5 8
6 13
...
45 1836311903
Warning: overflow encountered in long_scalars
46 -1323752223

Note:









No overflow when using Python int types









Program with Python's int type for integers

The program now avoids arrays and makes use of three int objects (which automatically changes to long when needed):

import sys
N = int(sys.argv[1])
xnm1 = 1                           # "x_n minus 1"
xnm2 = 1                           # "x_n minus 2"
n = 2
while n <= N:
    xn = xnm1 + xnm2
    print 'x_%d = %d' % (n, xn)
    xnm2 = xnm1
    xnm1 = xn
    n += 1

Run with \( N=200 \):

x_2 = 2
x_3 = 3
...
x_198 = 173402521172797813159685037284371942044301
x_199 = 280571172992510140037611932413038677189525
x_200 = 453973694165307953197296969697410619233826

Limition: your computer's memory









New problem setting: exponential growth with limited environmental resources

The model for growth of money in a bank has a solution of the type $$ x_n = x_0C^n \quad (= x_0e^{n\ln C})$$

Note:









Modeling growth in an environment with limited resources

Initially, when there are enough resources, the growth is exponential: $$ x_n = x_{n-1} + {r\over 100}x_{n-1}$$

The growth rate \( r \) must decay to zero as \( x_n \) approaches \( M \). The simplest variation of \( r(n) \) is a linear: $$ r(n) = \varrho \left(1 - {x_n\over M}\right) $$

Observe: \( r(n)\approx \varrho \) for small \( n \) when \( x_n\ll M \), and \( r(n) \rightarrow 0 \) as \( x_n\rightarrow M \) and \( n \) is big

Logistic growth model:

$$ x_n = x_{n-1} + {\varrho\over 100} x_{n-1}\left(1 - {x_{n-1}\over M}\right)$$ (This is a nonlinear difference equation)









The evolution of logistic growth

In a program it is easy to introduce logistic instead of exponential growth, just replace

x[n] = x[n-1] + p/100.0)*x[n-1]

by

x[n] = x[n-1] + (rho/100.0)*x[n-1]*(1 - x[n-1]/float(M))













The factorial as a difference equation

The factorial \( n! \) is defined as $$ n(n-1)(n-2)\cdots 1,\quad 0!=1$$

The following difference equation has \( x_n=n! \) as solution and can be used to compute the factorial: $$ x_n = nx_{n-1},\quad x_0 = 1 $$









Difference equations must have an initial condition









Have you ever though about how \( \sin x \) is really calculated?









Would you expect these fantastic mathematical results?

Amazing result by Gregory, 1667:

$$ \sin x = \sum_{k=0}^\infty (-1)^k\frac{x^{2k+1}}{(2k+1)!}$$

Even more amazing result by Taylor, 1715:

$$ f(x) = \sum_{k=0}^\infty \frac{1}{k!}(\frac{d^k}{dx^k} f(0))x^k $$ For "any" \( f(x) \), if we can differentiate, add, and multiply \( x^k \), we can evaluate \( f \) at any \( x \) (!!!)









Taylor polynomials

Practical applications works with a truncated sum:

$$ f(x) \approx \sum_{k=0}^{{\color{red} N}} \frac{1}{k!}(\frac{d^k}{dx^k} f(0))x^k $$

\( N=1 \) is very popular and has been essential in developing physics and technology

Example:

$$ \begin{align*} e^x &= \sum_{k=0}^\infty \frac{x^k}{k!}\\ &\approx 1 + x + \frac{1}{2}x^2 + \frac{1}{6}x^3\\ &\approx 1 + x \end{align*} $$









Taylor polynomials around an arbitrary point

The previous Taylor polynomials are most accurate around \( x=0 \). Can make the polynomials accurate around any point \( x=a \): $$ f(x) \approx \sum_{k=0}^N \frac{1}{k!}(\frac{d^k}{dx^k} f(a))(x-a)^k $$









Taylor polynomial as one difference equation

The Taylor series for \( e^x \) around \( x=0 \) reads $$ e^x= \sum_{n=0}^\infty {x^n\over n!} $$

Define $$ e_{n}=\sum_{k=0}^{n-1} \frac{x^k}{k!} = \sum_{k=0}^{n-2} \frac{x^k}{k!} + \frac{x^{n-1}}{(n-1)!}$$

We can formulate the sum in \( e_n \) as the following difference equation: $$ e_n = e_{n-1} + \frac{x^{n-1}}{(n-1)!},\quad e_0=0 $$









More efficient computation: the Taylor polynomial as two difference equations

Observe: $$ \frac{x^n}{n!} = \frac{x^{n-1}}{(n-1)!}\cdot \frac{x}{n} $$

Let \( a_n = x^n/n! \). Then we can efficiently compute \( a_n \) via $$ a_n = a_{n-1}\frac{x}{n},\quad a_0 = 1$$

Now we can update each term via the \( a_n \) equation and sum the terms via the \( e_n \) equation: $$ \begin{align*} e_n &= e_{n-1} + a_{n-1},\quad e_0 = 0,\ a_0 = 1\\ a_n &= \frac{x}{n} a_{n-1} \end{align*} $$

See the book for more details









Nonlinear algebraic equations

Generic form of any (algebraic) equation in \( x \):

\[ f(x)=0 \]

Examples that can be solved by hand:

\[ ax + b = 0 \] \[ ax^2 + bx + c = 0 \] \[ \sin x + \cos x = 1 \]









Newton's method for finding zeros; illustration













Newton's method for finding zeros; mathematics

Newton's method.

Simpson (1740) came up with the following general method for solving \( f(x)=0 \) (based on ideas by Newton): $$ x_n = x_{n-1} - {f(x_{n-1})\over f'(x_{n-1})},\quad x_0 \hbox{ given} $$

Note:









A program for Newton's method

Quick implementation:

def Newton(f, x, dfdx, epsilon=1.0E-7, max_n=100):
    n = 0
    while abs(f(x)) > epsilon and n <= max_n:
        x = x - f(x)/dfdx(x)
        n += 1
    return x, n, f(x)

Note:









An improved function for Newton's method

Only one \( f(x) \) call in each iteration, optional storage of \( (x,f(x)) \) values during the iterations, and ensured float division:

def Newton(f, x, dfdx, epsilon=1.0E-7, max_n=100,
           store=False):
    f_value = f(x)
    n = 0
    if store: info = [(x, f_value)]
    while abs(f_value) > epsilon and n <= max_n:
        x = x - float(f_value)/dfdx(x)
        n += 1
        f_value = f(x)
        if store: info.append((x, f_value))
    if store:
        return x, info
    else:
        return x, n, f_value









Application of Newton's method

$$ e^{-0.1x^2}\sin ({\pi\over 2}x) =0\] Solutions: $x=0, \pm 2, \pm 4, \pm 6, \ldots$ $$

Main program:

from math import sin, cos, exp, pi
import sys

def g(x):
    return exp(-0.1*x**2)*sin(pi/2*x)

def dg(x):
    return -2*0.1*x*exp(-0.1*x**2)*sin(pi/2*x) + \
           pi/2*exp(-0.1*x**2)*cos(pi/2*x)

x0 = float(sys.argv[1])
x, info = Newton(g, x0, dg, store=True)
print 'Computed zero:', x

# Print the evolution of the difference equation
# (i.e., the search for the root)
for i in range(len(info)):
    print 'Iteration %3d: f(%g)=%g' % (i, info[i][0], info[i][1])









Results from this test problem

\( x_0=1.7 \) gives quick convergence towards the closest root \( x=0 \):

zero: 1.999999999768449
Iteration  0: f(1.7)=0.340044
Iteration  1: f(1.99215)=0.00828786
Iteration  2: f(1.99998)=2.53347e-05
Iteration  3: f(2)=2.43808e-10

Start value \( x_0=3 \) (closest root \( x=2 \) or \( x=4 \)):

zero: 42.49723316011362
Iteration  0: f(3)=-0.40657
Iteration  1: f(4.66667)=0.0981146
Iteration  2: f(42.4972)=-2.59037e-79









What happened here??

Try the demo program src/diffeq/Newton_movie.py with \( x_0=3 \), \( x\in [-2,50] \) for plotting and numerical approximation of \( f'(x) \):

Terminal> python Newton_movie.py "exp(-0.1*x**2)*sin(pi/2*x)" \
          numeric 3 -2 50

Lesson learned:

Newton's method may work fine or give wrong results! You need to understand the method to interpret the results!









First step: we're moving to the right (\( x=4 \)?)













Second step: oops, too much to the right...













Third step: disaster since we're "done" (\( f(x)\approx 0 \))













Programming with sound

Tones are sine waves:

A tone A (440 Hz) is a sine wave with frequency 440 Hz: $$ s(t) = A\sin\left( 2\pi f t\right),\quad f = 440 $$

On a computer we represent \( s(t) \) by a discrete set of points on the function curve (exactly as we do when we plot \( s(t) \)). CD quality needs 44100 samples per second.









Making a sound file with single tone (part 1)

Sampled sine function for this tone: $$ s_n = A\sin\left( 2\pi f {n\over r}\right), \quad n=0,1,\ldots, m\cdot r $$

Code (we use descriptive names: frequency \( f \), length \( m \), 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









Making a sound file with single tone (part 2)

data = note(440, 2)
data = data.astype(numpy.int16)
max_amplitude = 2**15 - 1
data = max_amplitude*data
import scitools.sound
scitools.sound.write(data, 'Atone.wav')
scitools.sound.play('Atone.wav')









Reading sound from file

# echo: e[n] = beta*s[n] + (1-beta)*s[n-b]

def add_echo(data, beta=0.8, delay=0.002,
             sample_rate=44100):
    newdata = data.copy()
    shift = int(delay*sample_rate)  # b (math symbol)
    for i in xrange(shift, len(data)):
        newdata[i] = beta*data[i] + (1-beta)*data[i-shift]
    return newdata

Load data, add echo and play:

data = scitools.sound.read(filename)
data = data.astype(float)
data = add_echo(data, beta=0.6)
data = data.astype(int16)
scitools.sound.play(data)









Playing many notes

# put data1, data2, ... after each other in a new array:
data = numpy.concatenate((data1, data2, data3, ...))

The start of "Nothing Else Matters" (Metallica):

E1 = note(164.81, .5)
G = note(392, .5)
B = note(493.88, .5)
E2 = note(659.26, .5)
intro = numpy.concatenate((E1, G, B, E2, B, G))
...
song = numpy.concatenate((intro, intro, ...))
scitools.sound.play(song)
scitools.sound.write(song, 'tmp.wav')









Summary of difference equations

Solution of difference equations by simulation:

index_set = <array of n-values: 0, 1, ..., N>
x = zeros(N+1)
x[0] = x0
for n in index_set[1:]:
    x[n] = <formula involving x[n-1]>

Can have (simple) systems of difference equations:

for n in index_set[1:]:
    x[n] = <formula involving x[n-1]>
    y[n] = <formula involving y[n-1] and x[n]>

Taylor series and numerical methods such as Newton's method can be formulated as difference equations, often resulting in a good way of programming the formulas









Summarizing example: music of sequences

We will study two sequences: $$ x_n = e^{-4n/N}\sin(8\pi n/N)$$ and $$ x_n = x_{n-1} + q x_{n-1}\left(1 - x_{n-1}\right),\quad x = x_0 $$

The first has values in \( [-1,1] \), the other from \( x_0=0.01 \) up to around 1

Transformation from "unit" \( x_n \) to frequencies: $$ y_n = 440 + 200 x_n $$ (first sequence then gives tones between 240 Hz and 640 Hz)









Module file: soundeq.py

Try it out in these examples:

Terminal> python soundseq.py oscillations 40
Terminal> python soundseq.py logistic 100

Try to change the frequency range from 200 to 400.