Recurrence relations
A linear recurrence relation
You might have previously encountered a recurrence relation of the form
\[ u_{n+1}=au_n+b \tag{1}\] where \(a\) and \(b\) are constants.
Given numerical values for \(a\) and \(b\) and an initial condition, \(u_0\), a sequence can be computed that is a solution to Equation 1. This type of task is laborious and well suited to a computer (see Figure 1).
The population of Dundonian hobbits is observed to be declining by 5% per year. To increase the population, it is planned that 1000 of the species will be released at the end of May each year.
Let \(u_n\) represent the population of the hobbits at the beginning of June, \(n\) years after the first annual reintroduction into the population.
Suppose that \(u_n\) and \(u_{n+1}\) satisfy the recurrence relation \[ u_{n+1}=au_n+b, \] where \(a\) and \(b\) are constants.
- State the values of \(a\) and \(b\).
- Explain whether or not the population of the Dundonian hobbit will stabilise in the long term.
- The population of Dundonian hobbits at the beginning of the reintroduction programme was estimated at 5000. Explain whether or not the population will ever exceed 10000.
Explore how the computed solution depend on model parameters as follows:
- set \(b=0\) by varying the parameter \(a\) identify solutions that:
- tend to zero monotonically
- oscillate about 0
- blow up
- set b>0
- show that the solution converges to a non-zero value in the case where \(0<a<1\).
- show that the solution is oscillatory for \(-1<a<0\).
#| standalone: true
#| components: [viewer]
#| viewerHeight: 500
from shiny import App, Inputs, Outputs, Session, render, ui
from shiny import reactive
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
from scipy.integrate import odeint
app_ui = ui.page_fluid(
ui.layout_sidebar(
ui.sidebar(
ui.input_slider(id="a",label="a",min=-1.0,max=3,value=0.1,step=0.001),
ui.input_slider(id="b",label="b",min=0.0,max=15.0,value=10.0,step=0.01),
ui.input_slider(id="u0",label="u_0",min=0.0,max=20.0,value=5.0,step=1.0),
ui.input_slider(id="T",label="Number of iterations",min=0.0,max=60.0,value=20.0,step=1.0),
),
ui.output_plot("plot"),
),
)
def server(input, output, session):
@render.plot
def plot():
fig, ax = plt.subplots()
#ax.set_ylim([-2, 2])
# Filter fata
a=float(input.a())
b=float(input.b())
u0=float(input.u0())
T=int(input.T())
# Define rhs of LV ODEs
def rhs_pop_model(y,t,a,b):
rhs=a*y+b
return rhs
def DiscreteSol(rhs_pop_model,y_0,t,a,b):
y=np.zeros_like(t,dtype=float)
y[0]=y_0
for i in t:
if i>0:
y[i]=rhs_pop_model(y[i-1],t[i],a,b)
return y
# Define discretised t domain
t = np.arange(0, T, 1)
# define initial conditions
init_cond=u0
# Compute numerical solution of ODEs
sol1 = DiscreteSol(rhs_pop_model,init_cond,t,a,b)
# Plot results
y=sol1
ax.plot(t,y)
ax.set_xlabel('$n$')
ax.set_ylabel('$u_n$')
plt.grid()
plt.show()
app = App(app_ui, server)The logistic map
The recurrence relation explored in Figure 1 is linear (the right-hand-side is a linear function of \(u_n\)). When the model is generalised much richer dynamical behaviours can be observed. One famous example is the logistic map, where the governing equation can be written as \[ u_{n+1}=ru_n(1-u_n). \tag{2}\]
Note that the right-hand side is now a quadratic function of \(u_n\).
You can explore the solutions to Equation 2 using Figure 2.
- show that when \(0<r<1\) the solution converges monotonically to 0.0.
- show that when \(0<r<2\) the solution converges monotonically to a non-zero value.
- show that when \(2<r<3\) the solution is oscillatory and converges to a non-zero value.
- show that when \(r=3.2\) the solution is periodic and repeats every second step
- show that when \(r=3.47\) the solution is periodic and repeats every fourth step
- show that when \(r=3.7\) that the solution is neither periodic nor reaches a steady value.
- use the \(r \ \textrm{zoomed}\) and \(u \ \textrm{zoomed}\) sliders to magnify the third figure. Can you see self similarity (i.e. at fine scales the bifurcation structure looks similar to that at large scales?).
The logistic map provides one of the simplest mathematical formulations of a phenomenon known as chaos. Whilst a precise definition of chaos involves some technical concepts, chaotic systems are broadly characterised by:
- having non-periodic, non-steady solution
- sensitivity to initial conditions
- appearing to be unpredictable even through they are deterministic.
- self similarity
#| standalone: true
#| components: [viewer]
#| viewerHeight: 500
from shiny import App, Inputs, Outputs, Session, render, ui
from shiny import reactive
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
from scipy.integrate import odeint
app_ui = ui.page_fluid(
ui.layout_sidebar(
ui.sidebar(
ui.input_slider(id="r",label="r",min=0.0,max=5.0,value=0.1,step=0.001),
ui.input_slider(id="u0",label="u_0",min=0.0,max=1.0,value=0.5,step=0.01),
ui.input_slider(id="T",label="Number of iterations",min=0.0,max=60.0,value=20.0,step=1.0),
ui.input_slider(id="r_range",label="r zoomed",min=0.0,max=4.0,value=[0.0,4.0],step=0.001),
ui.input_slider(id="u_range",label="u zoomed",min=0.0,max=1.0,value=[0.0,1.0],step=0.01),
),
ui.output_plot("plot"),
),
)
def server(input, output, session):
@render.plot
def plot():
fig, ax = plt.subplots(3,1)
#ax.set_ylim([-2, 2])
# Filter fata
r=float(input.r())
u0=float(input.u0())
T=int(input.T())
r_min=float(input.r_range()[0])
r_max=float(input.r_range()[1])
u_min=float(input.u_range()[0])
u_max=float(input.u_range()[1])
# Define rhs of logistic map
def logistic_map(y,t,r):
rhs=r*y*(1-y)
return rhs
def DiscreteSol(rhs_pop_model,y_0,t,r):
y=np.zeros_like(t,dtype=float)
y[0]=y_0
for i in t:
if i>0:
y[i]=rhs_pop_model(y[i-1],t[i],r)
return y
# Define discretised t domain
t = np.arange(0, T, 1)
# define initial conditions
init_cond=u0
# Compute numerical solution of ODEs
sol1 = DiscreteSol(logistic_map,init_cond,t,r)
# Plot results
y=sol1
ax[0].plot(t,y)
ax[0].set_xlabel('$n$')
ax[0].set_ylabel('$u_n$')
plt.grid()
# Parameters
n_iterations = 1000 # total iterations for each value of r
n_last = 100 # number of iterations to plot (for steady state)
r_values = np.linspace(0.0, 4.0, 10000) # range of r values
u0 = 0.5 # initial population (seed)
delta=0.25
#r_min=r-delta
#r_max=r+delta
r_values2 = np.linspace(r_min, r_max, 10000) # range of r values
# Initialize plot
x = np.full_like(r_values, u0)
x2 = np.full_like(r_values2, u0)
# Iterate and plot bifurcation diagram
for _ in range(n_iterations):
x = r_values * x * (1 - x) # logistic map function
x2 = r_values2 * x2 * (1 - x2) # logistic map function
if _ >= (n_iterations - n_last): # plot only steady state
ax[1].plot(r_values, x, ',k', alpha=0.25)
ax[1].plot([r,r],[0,1],'r--')
ax[2].plot(r_values2, x2, ',k', alpha=0.25)
ax[2].plot([r,r],[0,1],'r--')
# Labels and display
ax[1].set_title("Bifurcation Diagram")
ax[1].set_xlabel("$r$")
ax[1].set_ylabel("$u^*$")
ax[2].set_title("Bifurcation Diagram (zoomed in) ")
ax[2].set_xlabel("$r$")
ax[2].set_ylabel("$u^*$")
ax[2].set_xlim([r_min,r_max])
ax[2].set_ylim([u_min,u_max])
plt.show()
app = App(app_ui, server)At Dundee, core concepts from calculus (e.g. differential equations) and algebra that are needed to study dynamical systems are introduced in the modules Maths 1A and Maths 1B and developed further in the modules Maths 2A and Maths 2B.
At Level 2 in the module Discrete Maths you would be introduced to discrete dynamical systems (e.g. recurrence relations, Markov chains). In the modules Introduction to Programming and Computer Algebra and Dynamical systems you would be introduced to techniques that enable you to numerically analyse difference equations.
At Level 3 in the module Mathematical Biology you would consider discrete dynamical systems model applied to Biological systems.
At Level 4 we offer a number of honours projects that investigate chaotic systems (e.g. the Lorenz equations, the double pendulum)
You can find out more about these modules here.