ls30a.sagews

LS 30A related items

• iteration

• solving differential equations (desolve_odeint)

• plotting a timeseries

• plotting a trajectory

• plotting a vector field

Iteration

The word iterate loosely means "to do something repeatedly". Its mathematical definition is slightly more precise: to iterate means to repeatedly perform the same calculation or procedure, and to use the result of the previous calculation as the input to the next one.

The following is a quick example of that for a system that grows exponentially, doubling each time:

$$x_{i} = f(x_{i-1}) = 2x$$

f(x) = 2*x

x_i = 1
for i in srange(20):
    print x_i, 
    xnew = f(x_i)
    x_i = xnew # overwriting variables

1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536 131072 262144 524288

This worked well, but what if we wanted to view a time series of the result? Instead of just overwriting the variable each time, we can just append each entry to a list and then plot them afterwards.

# appending to list
f(x) = 2*x

x_i = 1
data = [x_i] # initialize list with initial condition
for i in srange(20):
    xnew = f(x_i)
    x_i = xnew
    data.append(x_i) # append data here

print data
list_plot(data)

[1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576]

if we plot the data, we can easily see that this is an exponential function.

Solving a Differential Equation

In LS30, we will encounter two types of differential equations: ordinary differential equations (ODEs) and delay differential equations. This section will concentrate on solving ODEs. ODEs have the mathematical form $$\dot{x} = f(x),$$ where $\dot{x}$ is the rate of change of the states $x$.

Generally, we will go through these steps, in about this order:

1. Define the differential equation, the initial conditions and the timespan of interest.

2. Solve the differential equation

3. Display the results as either a time series or a trajectory.

Let's illustrate this process using the Lotka-Volterra (Shark-Tuna) model. We'll use the convention that the data will be represented in the order (N,P).

The Lotka-Volterra model has the pair of differential equations

$$N' = 0.5N - 0.01NP$$$$P' = 0.5(0.01NP) - .2P$$

and we want to simulate a trajectory of this system with an initial condition of (50, 75) and with a timespan from (0,100) time units and a stepsize of 0.1.

# Step 1

# Define the differential equation
var('N','P')
DEs = [0.5*N - 0.01*N*P, 0.5*0.01*N*P - .2*P] # N' is first, followed by P'.

# set the initial conditions
ic = [50, 75] # N=50, P=75

# set the timespan
t = srange(0,100,.1) # 1000 data points with a step size of 0.1 usually works pretty good for most simulations. Try this first.

(N, P)

# Step 2

# solve the differential equation
sol = desolve_odeint(DEs, dvars=[N,P], ics=ic, times=t)

There are a couple of common pitfalls in step 2:

1. switching the order of the variables in dvars

2. starting with invalid or inappropriate initial conditions for ic

3. using a timespan that either is too short to see the behavior we're interested in, has so many points that it takes too long to calculate, or has a timestep that is too large so the solution is not very accurate.

NB: Many other people prefer to write the differential equations straight into the differential equation, like this:

sol = desolve_odeint([0.5*N - 0.01*N*P, 0.5*0.01*N*P - .2*P], dvars=[N,P], ics=ic, times=t)

This is fine. It just means that if you need to change parameters or write something different, you may have to track down changes over a lot of cells.

# Step 3

# plotting the time series
p = list_plot(zip(t, sol[:,0]), legend_label="prey") + list_plot(zip(t, sol[:,1]), color="red", axes_labels=['time','N,P'], legend_label="predator")
p.set_legend_options(loc='lower right') # optional, sets the legend location to be in the bottom right hand corner.
show(p)

# plotting trajectory
list_plot(sol, axes_labels=['N','P']) # it's especially important to label the axes in a trajectory.

Plotting a Timeseries

A time series is a series of data plotted on the vertical axis with corresponding times on the horizontal axis. Usually to do this, you'll need a list/vector of times, and a list (or multiple lists) of data and then you'll zip them together and then plot them.

times = srange(10)
data = [1,6,3,8,3,7,2,8,-3,0]
list_plot(zip(times, data)) # list_plot(times, data)) returns an error.

In the previous cell, list_plot(times, data) returns an error because list_plot only accepts lists of coordinates, not individual lists.

In situations where you just want to the view the elements in a list, you can just use the list_plot function with the list as the only argument. This plots the elements of the list on the y axis and the index of each element on the x axis.

list_plot(data)

Plotting a Trajectory

A trajectory (sometimes referred to outside of this course as a phase portrait) is a plot of the state-data from a dynamical system. In contrast to a timeseries, trajectories do not explicitly show what time the data is from.

If the trajectory that you want to plot came from the desolve_odeint command, then you're all set. Just use list_plot(sol) and everything will come out alright. Since the output of desolve_odeint is a numpy array with each of the columns corresponding to the time series of each state, sol already has the "zipped list" characteristic built in.

# Romeo and Juliet state space model
var('R')
var('J')

DEs = [R -J, R + 1*J]
ic = [1,2]
t = srange(0,3,.1)

sol = desolve_odeint(DEs, dvars=[R,J], ics=ic, times=t)
list_plot(sol, axes_labels=['R','J'])

R J

If you print out the content of sol itself, it is a little more clear what is going on:

sol

array([[ 1. , 2. ], [ 0.87898366, 2.30963231], [ 0.71174545, 2.63676729], [ 0.49174825, 2.97804927], [ 0.21217373, 3.32906695], [ -0.13398915, 3.68421713], [ -0.5538318 , 4.03656472], [ -1.0543872 , 4.37770114], [ -1.64246138, 4.69760392], [ -2.32443278, 4.98450092], [ -3.10601662, 5.22474315], [ -3.99199168, 5.40269097], [ -4.98588737, 5.50061896], [ -6.08962912, 5.49864596], [ -7.30314068, 5.37469722], [ -8.62390262, 5.10450668], [-10.04646714, 4.66166813], [-11.56193004, 4.01774508], [-13.15736237, 3.14244976], [-14.81520505, 2.00390273], [-16.5126319 , 0.56898513], [-18.22088803, -1.1962037 ], [-19.90461277, -3.32577003], [-21.52115869, -5.85331732], [-23.01992099, -8.81109219], [-24.34169434, -12.22897127], [-25.41807715, -16.13327606], [-26.17094754, -20.54540406], [-26.51203795, -25.48026558], [-26.34264002, -30.94451792]])

Notice that up at the top, 1 and 2 are the initial conditions. Accordingly, they are at the beginning of the array columns. (Remember that each column is a timeseries for a single state)

list_plot(zip(t, sol[:,0]), axes_labels=['t','R'], figsize=4)
list_plot(zip(t, sol[:,1]), axes_labels=['t','J'], figsize=4)

Likewise, since each column of sol is a time series, we can also zip them together and achieve the same result as if we had just two separate lists zipped together.

list_plot(zip(sol[:,0],sol[:,1]), axes_labels=['R','J'], figsize=4) # plots trajectory

Plotting a Vector Field

A differential equation has a function that maps each coordinate of the state space to a change vector. (or the tangent bundle for the mathematicians) We can visualize 2D vector fields easily using the command plot_vector_field. The syntax is very similar to the plot command.

var('x','y')
des = [2*x*y +y, 5+y/x ]
plot_vector_field(des, (x,-5,5), (y,-5,5), axes_labels=['x','y'])

(x, y)

Since plot_vector_field returns a Graphics object, (plot, list_plot are also Graphics objects) we can overlay these together to get some nice results combining trajectories and the vector field. Notice how the trajectory seems to "follow" the path hinted at by the green change vectors.

# solve differential equation
t = srange(0,1.4,.01)
ic = (-4,-4)
sol = desolve_odeint(des, dvars=[x,y], ics=ic, times=t)

# build up figure
p = plot_vector_field(des, (x,-5,5), (y,-5,5), axes_labels=['x','y'], color="green")
p += list_plot(sol, color="red") # overlay list plot

show(p)

Plotting tricks

• limiting size of graph with xmin, xmax, ymin, ymax

• creating empty figure with list_plot([])

# Limiting size of graph. Useful for systems with exponential growth.
var('x','y')
sol = desolve_odeint([.2*x, -.3*y], dvars=[x,y], times=srange(0,100,.1), ics=[.1, 1])

(x, y)

# unrestricted plot has poorly scaled axes
list_plot(sol)

# restricted plot lets us focus on a particular area
list_plot(sol, xmin=-5, xmax=5, ymin=-5, ymax=5)

# plotting extra