README for WINPP
CONTENTS
INTRODUCTION
This version of WINPP comes as a winzipped file witha bunch of ODE examples
and some very minimal documentation. I don't have that much time to write
this stuff and never was very good at it. XPP which begat this program has
more extensive documentation although there are some substantial differences.
However, the basic structure of the input files is identical.
WINPP is a program designed to solve dynamical systems problems that take
several different forms:
- ordinary differential equations (ODEs)
- delay differential equations (DDEs)
- differential-algebraic equations (DAEs)
- Volterra integro-differential equations (IDEs)
- discrete dynamical systems (MAPs)
- boundary value problems (BVP)
The basic idea is to create an input file with any text editor. This file
tells WINPP the equations, parameters, functions, etc that are to be used
in the numerical simulation. In addition, it is possible to define plot parameters,
etc within the file. The interface is based on Win95/NT and will not run
under Windows 3.*. All aspects of the simulation can be changed within the
program using the interface, but the names of variables and the dimension
of the system cannot be changed nor can a new file be loaded in without exiting
first.
To run WINPP, just click on the icon for it and then use the file
selector to choose an ODE file. Alternatively, you can type it from the
DOS prompt followed by an ODE file name:
winpp pend.ode
LOGFILE
Every time you run WINPP, the program creates a file called WPP.LOG in the
directore where WINPP is found if called by clicking on the icon or in
the initiating directory if called from the DOS Prompt. This can be
used to see any comments etc that the program produces such as error
messages etc. If you run with a bad ODE file, this will tell you what
was wrong.
SUBWINDOWS
When you run WINPP, several windows come up including the main window
seen at the top of this document. There is the BROWSER
which has numerical data and several small windows which allow you
to change parameters, initial data, and boundary conditions. Typing
numbers into the window and clicking OK will set them. Clicking Go will
set them and run the simulation.
SIMPLE COMMANDS
- To integrate the equations, click on Run Go. Or just type Ctrl G.
- Change parameters by editing the parameter windows and then clicking
on OK in the window or clicking on Go. Change initial conditions this
way as well.
- To change the View, Click on Graphics View or Ctrl-V. Check the
Autofit box to have WINPP choose the limits. Choose 2 or 3 dimensional
plots.
- A quick way to plot a quantity versus time is to click Graphics XvsT
or Ctrl-X.
- Ctrl-E erases the window and Ctrl-R redraws it.
- File Exit will get you out as will Ctrl-Q.
- Ctrl-L will use the endpoint of the most recent integration as the
starting point of a new one.
- To change integration parameters, click on Numerics Int. Pars. or
type Ctrl I.
- Click on Edit Copy to copy the picture to the clipboard and then
paste it into your own documents or click Ctrl-C.
- Use the Browser to look at the numbers from your simulation, save
the data in space delimited ascii format (to import with EXCEL, etc)
and read in similar data to plot.
- The only hardcopy currently avaiable is POSTSCRIPT - click File
Print to create a postscript black and white or color picture.
INPUT FILES
The standard convention for WINPP is to name files with the extension
ODE. I will call such files ODE files. They all have the following
format. Note that there should be no spaces between a number and an
equals sign in the definions of parameters and initial values of
phase-space variables.
ODE File Format
# comment line - name of file, etc
d<name>/dt=<formula>
<name>'=<formula>
<name>(t)=<formula>
volt <name>=<formula>
<name>(t+1)=<formula>
x[n1..n2]' = ...[j] [j-1] ... <-- Arrays
markov <name> <nstates>
{t01} {t02} ... {t0k-1}
{t10} ...
...
{tk-1,0} ... {tk-1 k-1}
aux <name>=<formula>
!<name>=<formula>
<name>=<formula>
parameter <name1>=<value1>,<name2>=<value2>, ...
wiener <name1>, <name2>, ...
number <name1>=<value1>,<name2>=<value2>, ...
<name>(<x1>,<x2>,...,<xn>)=<formula>
table <name> <filename>
table <name> % <npts> <xlo> <xhi> <function(t)>
global sign {condition} {name1=form1;...}
init <name>=<value>,...
<name>(0)=<value>
bdry <expression>
0= <expression> <--- For DAEs
solv <name>=<expression> <------ For DAEs
special <name>=conv(type,npts,ncon,wgt,rootname)
fconv(type,npts,ncon,wgt,rootname,root2,function)
sparse(npts,ncon,wgt,index,rootname)
fsparse(npts,ncon,wgt,index,rootname,root2,function)
# comments
@ <name>=<value>, ...
done
INTEGRAL EQUATIONS
The general integral equation:
u(t)=f(t)+int(0,t) K(t,s,u(s))ds
becomes
u = f(t) + int{K(t,t',u)}
The convolution equation:
v(t) = exp(-t) + int(0,t) e^{-(t-s)^2}v(s) ds
would be written as:
v(t) = exp(-t) + int{exp(-t^2)#v}
If one wants to solve, say,
u(t) = exp(-t) + int(0,t) (t-t')^{-mu} K(t,t',u(t'))dt'
where 0 <= mu < 1 , the form is:
u(t)= exp(-t) + int[mu]{K(t,t',u}
and for convolutions, use the form:
u(t)= exp(-t) + int[mu]{w(t)#u}
OPTIONS
The format for changing the options is:
@ name1=value1, name2=value2, ...
where { name} is one of the following and { value} is either an
integer, floating point, or string. (All names can be upper or lower
case). The first option can only be set outside the program.
- MAXSTOR={ integer} sets the total number of time steps that
will be kept in memory. The default is 5000. If you want to perform
very long integrations change this to some large number.
The remaining options can be set from within the program. They are
- XP=name sets the name of the variable to plot on the x-axis.
The default is { T}, the time-variable.
- YP=name sets the name of the variable on the y-axis.
- ZP=name sets the name of the variable on the z-axis (if the plot
is 3D.)
- AXES={2,3} determine whether a 2D or 3D plot will be
displayed.
- TOTAL=value sets the total amount of time to integrate the
equations (default is 20).
- DT=value sets the time step for the integrator (default is 0.05).
- NJMP={ integer} tells XPP how frequently to output the
solution to the ODE. The default is 1, which means at each
integration step.
- T0=value sets the starting time (default is 0).
- TRANS=value tells XPP to integrate until { T=TRANS} and then
start plotting solutions (default is 0.)
- NMESH={ integer} sets the mesh size for computing nullclines
(default is 40).
- METH={ discrete,euler,modeuler,rungekutta,adams,gear,volterra, backeul,
qualrk,stiff,cvode } sets the integration method (see below; default is
Runge-Kutta.
- DTMIN=value sets the minimum allowable timestep for the Gear
integrator.
- DTMAX=value sets the maximum allowable timestep for the Gear
integrator
- { JAC_EPS=value, NEWT_TOL=value, NEWT_ITER=value} set
parameters for the root finders.
- ATOLER=value sets the absolute tolerance for CVODE.
- TOLER=value sets the error tolerance for the Gear, adaptive RK,
and stiff integrators. It is the relative tolerance for CVODE.
- BOUND=value sets the maximum bound any plotted variable can
reach in magnitude. If any plottable quantity exceeds this, the
integrator will halt with a warning. The program will not stop
however (default is 100.)
- DELAY=value sets the maximum delay allowed in the integration
(default is 0.)
- XLO=value,YLO=value,XHI=value,YHI=value set the limits for
two-dimensional plots (defaults are 0,-2,20,2 respectively.) Note that
for three-dimensional plots, the plot is scaled to a cube with
vertices that are +/-1 and this cube is rotated and projected onto
the plane so setting these to +/-2 works well for 3D plots.
- XMAX=value, XMIN=value, YMAX=value, YMIN=value, ZMAX=value, ZMIN=value
set the scaling for three-d plots.
- POIMAP={ section,maxmin } sets up a Poincare map for
either sections of a variable or the extrema.
- POIVAR=name sets the variable name whose section you are
interested in finding.
- POIPLN=value is the value of the section; it is a floating
point.
- POISGN={ 1, -1, 0 } determines the direction of the
section.
- POISTOP=1 means to stop the integration when the section is
reached.
- RANGE=1 means that you want to run a range integration (in batch
mode).
- RANGEOVER=name, RANGESTEP=number, RANGELOW=number, RANGEHIGH=number,
RANGERESET={ Yes,No}, RANGEOLDIC={ Yes,No} all correspond to
the entries in the range integration option.
- TOR_PER=value, defined the period for a toroidal phasespace and
tells WINPP that there will be some variables on the circle.
- FOLD=name, tells WINPP that the variable <name> is to be
considered modulo the period. You can repeat this for many variables.
RESERVED WORDS
You should be aware of the following
keywords that should not be used in your ODE files for anything other
than their meaning here.
sin cos tan atan atan2 sinh cosh tanh
exp delay ln log log10 t pi if then else
asin acos heav sign ceil flr ran abs
max min normal besselj bessely erf erfc
arg1 ... arg9 @ $ + - / * ^ ** shift
| > < == >= <= != not # int sum of i'
These are mainly self-explanatory. The nonobvious ones are:
- { heav(arg1)} the step function, zero if { arg1<0} and 1 otherwise.
- { sign(arg)} which is the sign of the argument (zero has sign 0)
- { ran(arg)} produces a uniformly distributed random number
between 0 and { arg.}
- { besselj, bessely } take two arguments, n,x and return
respectively, J_n(x) and Y_n(x), the Bessel functions.
- { erf(x), erfc(x)} are the error function and the
complementary function.
- { normal(arg1,arg2)} produces a normally distributed random number
with mean { arg1} and variance { arg2}.
- { max(arg1,arg2)} produces the maximum of the two arguments
and { min} is
the minimum of them.
- { if(<exp1>)then(<exp2>)else(<exp3>)} evaluates { <exp1> }
If it is nonzero
it evaluates to { <exp2>} otherwise it is { <exp3>}. E.g. {
if(x>1)then(ln(x))else(x-1)} will lead to { ln(2)} if { x=2} and { -1
if x=0.}
- { delay(<var>,<exp>)} returns variable { <var>} delayed by the result of
evaluating { <exp>}. In order to use the delay you must inform
the program of the maximal possible delay so it can allocate storage.
(See the section on the NUMERICS menus.)
- { ceil(arg),flr(arg)} are the integer parts of{ <arg>} returning the
smallest integer greater than and the largest integer less than { <arg>}.
- { t } is the current time in the integration of the differential equation.
- { pi} is pi
- { arg1, ..., arg9} are the formal arguments for functions
- { int, #} concern Volterra equations.
- { shift(<var>,<exp>)} This operator evaluates the expression
{ <exp>} converts it to an integer and then uses this to indirectly
address a variable whose address is that of { <var>} plus the
integer value of the expression. This is a way to imitate arrays in
WINPP. For example if you defined the sequence of 5 variables, {
u0,u1,u2,u3,u4} one right after another, then { shift(u0,2)} would
return the value of { u2.}
- { sum(<ex1>,<ex2>)of(<ex3>)} is a way of summing up things.
The expressions { <ex1>,<ex1>} are evaluated and their integer
parts are used as the lower and upper limits of the sum. The index of
the sum is { i'} so that you cannot have double sums since there is
only one index. { <ex3>} is the expression to be summed and will
generally involve { i'.} For example { sum(1,10)of(i')} will
be evaluated to 55. Another example combines the sum with the shift
operator. { sum(0,4)of(shift(u0,i'))} will sum up { u0} and the
next four variables that were defined after it.
EXAMPLES
Here I present some example ODE files that you are welcome to try out.
EXAMPLE 1. DAMPED PENDULUM
This is a second order system rewritten as
a pair of first order equations. The friction is mu, the length of
the pendulum is l, mass is m, and gravity is
g.
# damped pendulum pend.ode
# equations:
dx/dt = xp
dxp/dt = (-mu*xp-m*g*sin(x))/(m*l)
# energy
pe=m*g*(1-cos(x))
ke=.5*m*l*xp^2
# auxiliary plottable quantities
aux P.E.=pe
aux K.E.=ke
aux T.E=pe+ke
# parameters
param m=10,mu=2,g=9.8,l=1
param scale=0.0083333
@ bounds=1000
# initial position
x(0)=2
# end of the equation file
done
COMMENTS:
- I have defined some extra quantities called pe and ke. These are
called "fixed variables" and are defined in terms of phase variables
and parameters etc. They are used internally but not accessible to the
user. Only phase variables (those that obey dynamics) and auxiliary
quantities can be plotted and viewed.
- I have defined some auxiliary variables that can be viewed by the
user.
- Run this and look at a 3D view with X along the X-axis, XP along the
Y-axis, and T along the Z-axis; use the Autofit option.
EXAMPLE 2. Stochastic equations
In this example a differential equation for voltage is coupled to a
single channel that stochastically jumps between two states, 0 and 1.
# hhnasing.ode
# single channel sodium channel with greatly enhanced single channel
# conductance.
dv/dt = -gl*(v-vl) -gna*m*(v-vna)+I
markov m 2
{0} {am(v)}
{bm(v)} {0}
#
init v=0,m=0
#
am(v)=PHI*.1e0*(.25e+02-v)/(EXP(.1e0*(.25e02-v))-.1e+01)
bm(v)=PHI*4.0e0*EXP(-v/1.8e01)
# parameters
par i=0,vna=115,vl=0,gna=.1,gl=.3,phi=1
#
done
COMMENTS:
- The random variable, m, has two states open (1) and closed (0). The
rates of switching between the states are defined by a 2x2 array (or if
there were 5 states, 5x5) of switching rates, a_ij. The probability of
jumping from state i to state j is a_ij*dt. In WINPP, the diagonal
terms are ignored since the row sums must always be 1. Here am(v) is
the voltage-dependent rate of jumping from state 0 to state 1.
- WINPP is case-insensitive so EXP and Exp and eXp are all the
exponential function.
- Floating point exponents are of the form #e+#, #e-#, #e#, where # is
any number.
- The functions am(v) and bm(v) are defined by a simple obvious
statement.
- Run this simulation and look at the voltage and the channel states
EXAMPLE 3: Delay-differential equations
In this example, a nonlinear scalar feedback system is defined.
# delayed feedback: delay.ode
# declare all the parameters, initializing the delay to 3
p tau=3,b=4.8,a=4,p=-.8
# define the nonlinearity
f(x)=1/(1+exp(-x))
# define the right-hand sides; delaying x by tau
dx/dt = -x + f(a*x-b*delay(x,tau)+p)
x(0)=1
@ total=40,xhi=40,ylo=0,yhi=1,delay=10
# done
done
COMMENTS:
- The "@" symbol is used to initialize internal plot and simuation
parameters so that theuser doesnt have to do it manually within the
program. Here I have set the total integrateion time to 40 and also
made the x and y axes optimal for the simulation. Finally, I have set a
parameter called "delay" which tells WINPP to allocate enough storage
so that it can access delays as far back as 10. I only needed 3 for
this simulation.
- Within a right-hand side the function delay(x,tau) produces,
x(t-tau).
- Run this. Then move the mouse into the plot window, hold down the
button and move around. Note how the (x,y) values are displayed. Find
the period of the oscillation. It is about 9.7
EXAMPLE 4. Integro-differential
equations
Here is a simple integro-diffential equation.
# tstvol2.ode
# Volterra example from Peter Linz's book
# solution is f1=1,f2=t but it is unstable
f1(0)=1
f1(t)=exp(-t)-t^2/2 + int{exp(-t)#f1*f1}+int{f2}
f2(t)=t-t^4/24+int{t#f2*f2/(1+f1*f1)}
@ yplot=f2,total=5,xhi=5,yhi=5
done
COMMENTS:
- The first integral operator convolves exp(-t) with f1(t)^2 and the
second one convolves t with f2(t)^2/(1+f1(t)^2)
- yplot=f2 tell WINPP to plot f2(t).
- Run this. Then click on Graphics Add/edit. In this dialog, click on
Add New. Make the Y axis of the new curve f1 and make the color Red.
Click on Done and then click on Graphics Redraw (Ctrl-R).
EXAMPLE 5. Differential-algebraic
equations
This is a simple DAE. We want to solve:
dx/dt=y
y+x=1
# dae_ex1.ode
x'=y
0= x+y-1
x(0)=0
solv y=1
aux yy=y
done
COMMENTS:
- The expression "0=" starting a line tells WINPP that the remaining
expression must be set to 0. The expression "solv y=1" tells WINPP that
it "y" is an algebraic variable and that an initial guess should be 1.
- The variable that is to be solved for, y, is hidden from the user
and in order to see it, the plottable auxiliary variable yy is defined.
Here is a more complex example in which the derivative is implicitly
defined.
dx/dt + exp(dx/dt)=-x
# dae_ex2.ode
x'=xp
0= xp+exp(xp)+x
x(0)=-3.7182
solv xp=1
aux xdot=xp
@ ylo=-4,yhi=0
done
COMMENTS:
- Note that the root solver solves for xp and this is just dx/dt.
In this last example of a DAE, we solve a relaxation oscillator by
making the tolerance for solving the root somewhat large and thus
letting numerical errors take the problem over singulities.
# dae_ex3.ode
w'=v_
0= v_*(1-v_*v_)-w
solv v_=1
aux v=v_
@ NEWT_ITER=1000,NEWT_TOL=1e-3,JAC_EPS=1e-5,METH=qualrk
done
COMMENT:
- Plot v and w in a phase-plane
EXAMPLE 6. Boundary-value problems.
Here we solve a nonlinear boundary value problem:
x''=x^3 ;x(0)=1 ; x(1)=0
We rewrite it as a pair of first order equations. The ODE file is
# bvp.ode
x'=xp
xp'=a*x^3
bndry x-1
bndry x'
init x=0
par a=1
@ dt=.01,total=1,xhi=1,yhi=1,ylo=0
done
COMMENTS:
- To run this, Click on Run BVP Show and repeated attempts at solving
using shooting will be shown. If a solution is found the intermediate
solutions will be erased and the final solution will be shown. The
initial value of xp is free to be chosen and will hopefully be found
to satisfy the end condition.
- Boundary conditions are put in with the statement "bndry". The program
attempts to make them zero. A primed variable means to evaluate the
variable at the end of the integration while unprimed means to evaluate
at the beginning.
- The total integration time is set here to 1 reflecting the length of
the domain.
Here is a very complicated boundary value problem.
Find a value of omega solving:
d (a'' - ak^2 + a'/t -a/t^2) + a(1-a^2)=0
d (k'+k/t+2ka'/a) + (omega + q a^2) =0
with
limit (t->0) a(t)/t exists
k(0)=0
k(1)=0
a'(1)=0
Note that 4 conditions must be satisfied (a(0)=0 is implied by the
first condition) but the system is only 3rd order. However omega is a free
parameter, so we write a differential equation for it of the form:
omega'=0
Now the system is 4 dimensional:
# gberg.ode
init a=0.0118 ap=1.18 k=0 omeg=-0.19
par d=0.177 q=0.5 sig=3
# the odes...
a'=ap
ap'=a*k*k-ap/t+a/(t*t)-a*(1-a*a)/d
k'=-k/t-2*k*ap/a-(omeg+q*a*a)/d
omeg'=0
# the boundary conditions
bndry a-t*ap
bndry k
bndry ap'
bndry k'
# set it up for the user
@ xhi=1 t0=.01,dt=.01,total=.99
done
COMMENTS:
- Since a(t)/t is undefined in the computer as t->0, we start t at a
small value and insist that a(t) = a'(0)t for t small. We integrate
for a total of .99 since we started at t=.01.
- It is written as a system of first order equations.
- This equation arises when one looks for spiral wave solutions to a
system of reaction-diffusion equations.
EXAMPLE 7. Equations with arrays
WINPP cannot really handle arrays of equations, but the parser will
expand an input file with an expression like
x[2..5]'=-x[j]+x[j-1]
into
x2'=-x2+x1
x3'=-x3+x2
x4'=-x4+x3
x5'=-x5+x4
This makes it easy to solve discretized PDEs. For example, the simple
Fisher equation:
u_t = u_xx + u(1-u)
discretized into 40 parts on the interval 0 < x > 40.
# discretization of Fisher's equation fisher.ode
# with no flux boundary conditions
par h=1
u0'=u0*(1-u0)+(u1-u0)/(h^2)
u[1..39]'=u[j]*(1-u[j]) + (u[j-1]-2*u[j]+u[j+1])/(h^2)
u40'=u40*(1-u40)+(u39-u40)/(h^2)
init u0=.1
@ total=40,dt=.2,meth=qualrk
@ xhi=40,yhi=1,ylo=0,yp=u20
done
COMMENTS:
- [j] refers to the current index in the array expansion. [j+n]
where n is an integer is an allowable syntax. Thus
u[4..8]' = u[j-3] will create a set of 5 equations, the first
being u4'=u1 .
- This is a crude discretization
- WINPP has some nice tricks for plotting arrays of equations.
Click on Graphics Array plot. When the window appears, click on Edit.
Click on the Last Column scroll button until u40 is highlighted. This
chooses u0..u40 as variables to look at. Fill in RowSkip =1, ZMin=0,
Zmax=1.01,and
then OK. After a few seconds you should see a nice color plot showing
the appearance of the traveling front. Color is magnitude, the
horizontal axis is column # and the vertical axis is time. Copy will
copy it to the clipboard. Print lets you make a postscript version of
it in color or greyscaled.
- You can look at the profile of u0..u40 at
a fixed point in time by clicking on the Numerics Transpose item.
Choose u0 as the first column, 41 as the number of columns, 50 for row
1. This will temporarily replace u0 by the row of values u0..u40 that
was 50 rows down in the BROWSER. Click Ctrl X and choose u0 to see the
"spatial" profile.
Here is a neural network model that illustrates the sum function. It is
a simplified version of a model proposed by Hansel and Sompolinsky for
orientation tuning:
du_i/dt = -u_j + F( sum_j J(i-j) u_j)
We choose J(x)=c+cos(beta*x). The equations in WINPP are:
# chain of 20 neurons wcring.ode
param a=.25,beta=.31415926,c=.2
u[0..19]'=-u[j]+f(a*sum(0,19)of((c+cos(beta*([j]-i')))*shift(u0,i')))
f(x)=tanh(x)
done
COMMENTS
Here the sum command is used along with the shift operator. i'
always is the name of the index for the sum operator. The variables
u0..u19 are defined successively so shift(u0,i') accesses the iith
variable defined after u0. That is shift(u0,3) would be u3 since u3 is
the third variable defined after u0. The index in an array definition
must always have brackets around it, eg [j]. So ([j]-i') will be expanded
to (0-i') for j=0, (1-i') for j=1, etc.
Initialize one of the u_j's to 1 and integrate the equation (Ctrl
G)
Click on Graphics Array and select u0..u19 with the remaining
defaults. Note how a local region grows and suppresses the remaining
region. Try some different parameters, for example c,making it closer
to 1.
In this last array example, I define a discrete convolution using the
special operator.
# neural network nnet2.ode
tabular wgt % 25 -12 12 1/25
f(u)=1/(1+exp(-beta*u))
special k=conv(even,51,12,wgt,u0)
u[0..50]'=-u[j]+f(a*k([j])-c*delay(u[j],tau)-thr)
par c=8,a=9,beta=10,thr=1,tau=4
init u0=1,u1=1
@ total=10,delay=10,xhi=10,ylo=0,yhi=1,yp=u20
done
COMMENTS:
The tabular command defines a lookup table that is
treated as a one variable function. It has 25 elements, is defined for
-12 > x < 12 and takes the value of 1/25 everywhere. If you want
to make a table out of a function, the argument in the function is
always t so that if instead of 1/25, you had
exp(-abs(t/4)) the function would be exponential. Tables can also
be read in as data files.
The special command defines another function of one
variable that is defined only on the integers 0 .. 50 since
the second argument is 51. This function convolves the function
wgt with the variable u0..u50 . The jth
value of k is
sum(-12,12)of(wgt(i')*shift(u0,i'+j)). If i'+j is less
than 0 or greater than 50, then the even directive implies
that the index is reflected. So the -2 becomes 2 and 56 becomes 44.
Integrate this and use the Array Plot to view
the traveling pulse.
EXAMPLE 8: Discontinuous systems
In this example there are three variables, (u,v,m). The last is the mass
and when the variable u falls below 0.2, the mass is divided by two.
This is a model of the cell cycle and represents a discontinuous
differential equation -- the right hand side of the equation contains
an effective delta function. Here is the ODE file:
# tyson.ode
init u=.0075,v=.48,m=1
p k1=.015,k4=200,k6=2,a=.0001,b=.005
u'= k4*(v-u)*(a+u^2) - k6*u
v'= k1*m - k6*u
m'= b*m
global -1 {u-.2} {m=.5*m}
@ total=1000,dt=.25,meth=qualrk,xhi=1000,ylo=0,yhi=2,yp=m
done
COMMENTS:
The global directive has three parts. The idea is to look
for a condition to occur and then operate on the variables.The first part
of the global is a
direction; -1 means crossing through 0 with negative slope while +1
means crossing through zero with positive slope. The second argument is
the condition that is to cross 0. The last is a brace delimited list of
conditions separated by semicolons. Thus in the present example if
u-.2 crosses 0 with negative slope, the mass is halved.
Integrate and fire models can be designed with this type of command.
WINPP Commands
Running a simulation
Click on Run and then one of the options. To
abort a simulation click on STOP . Initial data are specified by typing
values in the Initial Data boxes or via the menus.
- Go just uses the current initial data and runs for the amount of time
specified by the Numerics parameters
- Last uses the end result of the prior intergation as initial
conditions of the current simulation.
- Mouse lets you click on initial conditions in the graphics window.
NOTE that both axes must be phase-variables; i.e., they must satisfy
differential equations.
- Mice lets you click on lots of initial conditions repeatedly. Click
on STOP to finish.
- New lets you specify initial data for variables without clicking in
the initial data box which uses a font that is too small for my aging
eyes.
- More time lets you add additional time to the integration.
- Range integrate lets you repeat an integration while varying a
parameter or phase-variable.
Use old ICs means that the same initial
data (except the state variable if that is varied) for each integration.
Reset storage means that the data buffer is cleared after each
integration. Cycle color runs through 10 colors for the plot
curves.
NOTES: (i)The parameters are reset to their original values after
the task is completed; (ii) Only the last simulation remains in the
data buffer unless Reset storage is disabled; (iii) STOP stops a
single integration and ABORT ALL stops all simulations.
- Shoot uses special initial data that are generated when fixed
points are found. See the section on equilibria.
- BVP solves boundary value problems. "Show" shows intermediate attempts
at solution. The boundary conditions are given within the ODE file or
typed into the Boundary Condition window. If solutions fail to be found
you can try to adjust the BVP tolerances
- Equilibria
WINPP uses a Newton Solver to find fixed points. Once the fixed point
is found, a window appears
that contains information about the fixed point. The numerical value is
found by scrolling down the fixed point slider. The eigenvalues are
also shown. The title indicated stability. Finally a summary of the
eigenvalue information appears at the top; r+ is the number
or real positive eigenvalues, etc. the Grab button will load
the fixed point into the initial data buffer. Sometimes the Netwon
Sover fails to converge. You can change some parameters in the numerics section or perhaps get a better guess.
In addition to computing eigenvalues and equilibria WINPP can use this
information to compute 1-dimensional invariant manifolds. If there is a
single real positive eigenvalue and/or a single real negative
eigenvalue, WINPP will use the eigenvectors corresponding to these
eigenvalues to approximate the unstable/stable manifolds. WINPP will
ask if you want the Invariant manifolds to be computed when such a
fixed point is encountered. Answering yes will produce up to 4 new curves.
Since the integration is set forever for these curves, usually
something will go out of bounds or tend to an attractor. Click on
STOP to abort the calculation and move on to the next invariant
manifold. Once you have computed them, the initial conditions that led
to their calculation are stored in memory. There are up to 4 possible
branches (two for stable and two for unstable). You can access these
branches from the Run Shoot menu. Choose any or all of the
four branches. Remember to compute the stable manifolds,
integrate backwards in time!
The Equilibria menu has several options:
Plotting and browsing.
THE BROWSER
WINPP keeps the results of simulations in a data buffer. You can view
the numbers from the simulation by looking at the Browser. Use the
arrow keys or scroll bars to go through the numbers. There are several
commands:
- Find will look for the closest value of one of the browser
quantities to the number specified.
- Write will save the numbers in an ascii space delimited file.
- Load will load data into the browser for plotting
- Get will use the numbers in the top row of the browser as initial
data.
- Replace lets you replace a column with some formula.
- Unreplace unreplaces the last replace
- AddCol lets you add a formula column to the browser. All subsequent
integrations will compute the formula.
- DelCol lets you delete any added columns
GRAPHICS OPTIONS
WINPP lets you plot any quantity versus any other quantity, plot multiple
curves, save curves from previous runs, look at arrays, and run
animations.
Here are the commands:
- View This sets the view for the graphs. You can choose 2D
and 3D views. For 3D I always let WINPP define the limits of the view.
For 2D Xlo, XHi, YLo, YHi define the limits of the graph.
Type in the labels you want on the axes, determine whether or not you
want origin axes and determine their values. Click OK. For 3D views,
Xmin, etc define the scales for the actual variables. These
are plotted in the unit cube which is then rotated by the amount
specified by Phi,Theta and this is projected to the 2D plane
and fit in the box defined by XLo, etc . Experiment with this
and see what happens.
- Fit will maximize the screen area for the given plot
- Zoom in lets you magnify a region with the mouse.
- Zoom out projects the current view into the zoomed boz.
- New Window produces a new window with its own graph and
menus. Kill it in the usual fashion. Note the small check box next to
it. This means it is the active window and all new integrations are
drawn in it. Click on the check box in any defined window to make it
active.
- X vs time This lets you plot any variable as a function of
time. It is a shortcut to use instead of the View command.
- Keep curve WINPP clears resets the data buffer after
every integration, thus you cant compare the results of two
simulations.
This option freezes the curve in memory. You can name the
curves and give them different colors as well as delete them.
Autofreeze will keep the results of every integration until the maximum
number of curves is reached. Show key shows a key in the window
that has each saved curve's legend and the color of its line. Mouse
Position will prompt you for the position of the key after you
click Done .
- Add/edit lets you plot multiple variables in the same window
for a given integration.
Choose a color and the variables to plot; you can
change and delete these whenever you want. Line style draws lines
between points while Points just draws the points. Scroll through the
curve list to edit different curves and scroll through the colors to
change them. The difference between Add/Edit and Keep Curve
is that the curves in Add/Edit change with every
integration while Keep curves are frozen and attached only to
the window in which they first appeared.
- Redraw,Erase are pretty obvious.
- Toons See the section on animation .
- Array plot was already described
Animation
Many years ago, as a teenager, I used to make animated movies using
various objects like Kraft caramels (``Caramel Knowledge'') and
vegetables (``The Call of the Wild Vegetables''). After I got my
first computer, I wanted to develop a language to automate computer
animation. As usual, things like jobs, family, etc got in the way and
besides many far better programmers have created computer assisted
animation programs. Thus, I abandoned this idea until I recently was
simulating a simple toy as a project with an undergraduate. I thought
it would be really cool if there were a way to pipe the output of a
solution to the differential equation into some little cartoon of the
toy. This would certainly make the visualization of the object much
more intuitive. There were immediately many scientific reasons that
would make such visualization useful as well. Watching the gaits of an
animal or the waving of cilia or the synchronization of many
oscillators could be done much better with animation
than two and three dimensional plots of the state variables.
With this in mind, I have developed a simple scripting language that
allows the user to make little cartoons and show them in a dedicated
window. Here is an example:
- Run the numerical simulation pend.ode.
- Click on the Graphics Toon menu item from the
main WINPP window.
- Click on the File item in the animation window that pops
up and find pend2.ani
- If the file is OK, click on |> item and the animation will
begin.
The other buttons in the animation window are:
- | | Pause
- | |< < Reset to the beginning
- Up Advance one frame
- Down Go back one frame
- Copy Copy current frame to clipboard.
Here are some AVI movies I have made from sample animations
DASL: Dynamical Animation Scripting Language
In order to use the animation components of WINPP you must first
describe the animation that you want to do. This is done by creating a
script with a text editor that describes the animation. You describe
the coordinates and colors of a number of simple geometric
objects. The coordinates and colors of these objects can depend on the
values of the variables (both regular and fixed, but not auxiliary) at
a given time. The animation then runs through the output of the
numerical solution and draws the objects based on that output. I will
first list all the commands and then give some examples.
Basically,
there are two different types of objects: (i) transient and (ii)
permanent . Transient objects have coordinates that are recomputed at
every time as they are changing with the output of the
simulation. Permanent objects are computed once and are fixed through
the duration of the simulation. The objects themselves are simple
geometric figures and text that can be put together to form the
animation.
Each line in the script file consists of a command or object followed
by a list of coordinates that must be separated by semicolons. For
some objects, there are other descriptors which can be optional.
Since color is important in visualization, there are two ways a color
can be described. Either as a formula which is computed to yield a
number between 0 and 1 or as an actual color name started with the
dollar sign symbol, $. The .ani file consists of a lines of
commands and objects which are loaded into XPP and played back on the
animation window. Here are the commands:
- dimension xlo;ylo;xhi;yhi. This sets the
dimension of the animation window coordinates. At start up the
coordinates are (0,0) and (1,1). (0,0) is the bottom left and (1,1) is
the upper right.
- transient. Commands following are transient so that they
are recomputed at each time slice.
- permanent. Commands following are permanent so that they
are computeds and drawn only once.
- line x1;y1;x2;y2;color;thickness. This draws a
line from the first two coordinates to the second two. color
is optional and can be a named color (starting with $):
$WHITE, $RED, $REDORANGE, $ORANGE, $YELLOWORANGE,
$YELLOW, $YELLOWGREEN, $GREEN, $BLUEGREEN,
$BLUE,$PURPLE, $BLACK
or it can be a decimal number between 0 and 1. In the latter case, the
color is of the same spectrum as the array plot with red=0 and
violet=1. thick is optional and is an integer indicating the
thickness of the line.
- rline x1;y1;color;thickness. This draws a line
relative to the previous line. color,thick are optional.
- rect x1;y1;x2;y2;color;thickness. This draws a rectangle
with coordinates lower left coordinates (x1,y1) and upper right
(x2,y2). color,thick are optional.
- frect x1;y1;x2;y2;color.This is a filled
rectangle.
color is optional.
- circ x1;y1;rad;color;thickness.This draws a circle
centered at (x1,y1) with radius rad. color,thick are optional.
- fcirc x1;y1;rad;color. This is a filled circle
with optional color.
- ellip x1;y1;rx;ry;color;thickness. This draws
an ellipse centered at (x1,y1) with horizontal axis, rx and vertical
ry. color, thick are optional.
- fellip x1;y1;rx;ry;color.This draws a filled
ellipse with optional color.
- text x1;y1;s. This plots the text string s
at position (x1,y1) using the current text attributes.
- vtext x1;y1;s;z. This allows you to plot a
string s with a numerical value z which is
computed on the fly during the animation.
- settext size;font;color. This sets the text
size (0-4), font (roman,symbol) and named color. Color is not optional
here and must be included.
- end. This tells WINPP to close the animation file.
Animation Examples
Here is the animation file for the pendulum with lots of comments.
# fancy animation file for pendulum
# this stuff is always visible and doesnt change
PERMANENT
# set text to middle sized roman font in purple color
settext 3;rom;$PURPLE
# place the text at the top of the screen
text .25;.9;Pendulum
# Draw a thick blue line across the middle to hang the pendulum
line 0;.5;1;.5;$BLUE;4
# the rest of this stuff depends on time
TRANSIENT
# here is the pendulum's rod - its angle depends on x from the ode file
line .5;.5;.5+.4*sin(x);.5-.4*cos(x);$BLACK;3
# here is the bob centered at the end of the rod
# the color of the bob is proportional to the kinetic energy ke!
fcircle .5+.4*sin(x);.5-.4*cos(x);.05;1.-scale*ke
# now make some small text
settext 1;rom;$BLACK
# type out the current value of the time
vtext .05;.1;t=;t
# now lets use a symbol font to plot out the current angle
settext 1;sym;$BLACK
vtext .05;.05;q=;x
end
Array example
The next example is of a large dynamical system that represents a set
of coupled excitable cells. The ODE file is called wave.ode
.
Here it is
# wave.ode
# pulse wave with diffusional coupling
param a=.1,d=.2,eps=.05,gamma=0,i=0
vv0(0)=1
vv1(0)=1
f(v)=-v+heav(v-a)
#
vv0'=i+f(vv0)-w0+d*(vv1-vv0)
vv[1..19]'=i+f(vv[j])-w[j]+d*(vv[j-1]-2*vv[j]+vv[j+1])
vv20'=i+f(vv20)-w20+d*(vv19-vv20)
#
w[0..20]'=eps*(vv[j]-gamma*w[j])
@ meth=qualrk,dt=.25,total=150,xhi=150
done
Run this and then load the animation file wave.ani.
The animation just plots the voltages as colored beads on a string. Their
horizontal position is the index and vertical is their voltage. The color
of the beads id the degree of recovery, w.
Since I have run
the simulation, I know that the recovery variables are between 0 and 1
and that the voltages are between -1 and 1. Here is the one-line
animation file for this effect:
# wave.ani
# animated wave
fcircle .05+.04*[0..20];.5*(vv[j]+1);.02;1-w[j]
end
This file uses the ``array'' capabilities of the parser to expand
the single line into 21 lines from 0 to 20. I just draw a filled
circle at scaled vertical coordinate and horizontal coordinate with a
small radius and colored according to the recovery variable. The
horizontal coordinate is expanded to be .05 + .04*0 , .05 +
.04*1 , etc; the vertical is .5*(vv0+1), .5*(vv1+1), etc;
and the color is 1-w0}, 1-w1 , etc. Thus this is interpreted as
a 21 line animation file. Try it to see what it looks like. It is a
simple matter to add a vertical scale and time ticker.
There are several other examples included in the distribution. Try for
example, lorenz2.ode and the animation file
lorenz.ani .
Permanent movies
To make a permanent AVI file like I have, you
must obtain some imaging software. I use HyperCam a shareware program
available for $30. All I do is tell HyperCam to focus on the animation
window, turn on the record feature and start the animation. Then I save
it as an AVI file. Simple and cheap.
File Menu
The File menu controls IO stuff.
The commands are:
- Write set. This lets you save the current state of the
program into a file name. When you reload it with Read set you
can pick up where you left off.
- Read set. Resets WINPP to a previous state.
- Print produces black and white or color POSTSCRIPT
files. Windows Printer support is not yet available.
- Locbif opens up LOCBIF a local bifurcation
analyzer.
- Quit
Edit Menu
There are two items:
- Edit RHS This lets you edit the right hand sides of the
equation. Choose the equation to edit and type in the right-hand side.
Click on Reset to reset the highlighted equation to its
previous state.
- Copy copies the current graphics window into the
clipboard.
Numerics Menu
This menu has many numerically associated commands, most importantly
the ones that control the main integration routines.
The commands are:
- Int. Pars. define all the integration parameters.
The Pick method box lets you choose from many integrators.
- Discrete. Use this for difference equations like the
discrete second order logistic equation tstlbffp.ode
- Euler. First order Euler
- Heun. Second order Euler
- RK4. 4th order fixed step Runge-Kutta
- Adams. 4th order fixed step Adams-Bashforth predictor
corrector
- Gear. Variable order Gear for stiff systems
- BackEul. Backward differentiation first order Euler
- Stiff. Another stiff integrator
- QualRK Variable step Runge-Kutta
- Cvode. LSODE stiff integrator -- highly recommended!
- Volterra. Integrator for Volterra integral equations.
The Tolerances except for Bounds are all used by the
adaptive step integrators. Bounds is a global bounds and if
any variable exceeds this, the integration is halted.
- Total is the total amount of time to integrate.
- Transient is the amount of time to integrate before
plotting and saving the data.
- Delta T is the time step for fixed step integrators and
the output for adaptive integrators. If it is negative, then
integration is reversed in time.
- Nout is the frequency of output for fixed step
integrators and should be 1 for adaptive integrators.
- Tstart is the starting time.
- Max Delay tells WINPP how much to allocate for delay
equations. The delay in your problem should be less than this number.
- Control parameters control certain internal numerical
code.
- BVP controls control the solution of boundary value
problems: BVP Iter is the maximum number of iterations to try to
solve the problem; BVP Eps controls the numerical
differentiation for Newtons method as applied to BVPs; BVP Tol
controls the tolerated error in the boundaries.
- Newton control controls the analogous things for newtons
Method applied to fixed points and DAEs. Newt Dh is controls
numerical differentiation and Newt Eps controls the error
tolerated.
- Delay control controls the parameters used to test
stability for delay equations. I am not sure it works....
- Colorize lets you color the trajectories of an
integration according to the value of some quantity.
The choices are Turn off which disables it, Velocity
which takes the L1 norm of the vector field, or User which
uses the user quantity chosen with the scroll bar. Autoscale
will optimize the color scale.
- Poincare map lets you look at the
output at certain times
depending ohe type of map. Whenever the map crosses 0, the time and
values of all the variables are stored.
There are 3 types of maps:
- Section which stores values whenever one of the variables
or other computed quantities crosses a chosen number.
- Max/Min occurs when a variable or other quantity hits a
local max or min.
- Periodic events are somewhat different; eacg time a
variable crosses a section, the time difference from the last such
crossing is recorded in the Time column of the browser. This
lets you compute interspike intervals.
Direction indicates whetjer to stop at every crossing (
Both ) when the crossing occurs with Positive slope or
Negative slope. The Variable entry tells WINPP which
variable to test and the section number is the value at which the
crossing should occur. Choosing t as the variable is treated
differently; whenever mod(t,value)= 0 , the point is recorded.
Stop on section tells WINPP to stop the integration when the
crossing occurs.
- Ruelle plot lets you plot a lagged version of a variable
against itself or any other variable. The lag or shift is in
terms of number of time steps, rather than the actual time. For example
if you solve a system with a delay of 3 and the time step is 0.05, then
to see x(t) vs x(t-3), choose a shift of 3/0.05 = 60.
- Stochastics only seeds the random number in this version.
- Torus lets you treat chosen variables as though they lie
on a circle. Thus if you have the equation
x' = 1
then it will grow linearly in time if you just integrate it. But if it
is treated as a "torus" variable with period of 2 Pi, then it will be
reset to 0 each time it crosses 2 Pi. This is useful for phase
equations such as phase.ode. Try this and make
all the variables 2 Pi periodic. Note that there are 2 stable
phaselocked solutions.
- Transpose lets you transpose the data buffer and was
described above.
Averaging This allows you to compute the adjoint for a
periodic solution and to average it agains interactions to compute the
averaged equations. If this means nothing to you, skip it; only a few
people in the world will use this. Basically, one looks at:
X1'=F(X1) + eps G(X1,X2)
X2'=F(X2) + eps G(X2,X1)
where each system has a periodic solution, X0(t), when eps=0.
The averaging routine in WINPP compute the effect on the phase of the
oscillation when eps > 0 and small:
H(z) = (1/T) integral{X*(t) G(X0(t),X0(t+z)),t=0..T}
where X*(t) is the adjoint solution. H(z) is a
T-periodic function that yields information on the effects of
one oscillator on the other.
To properly use the items in this menu, you must first compute a
periodic orbit. Set the total integration time as close to the period
as possible and compute exactly one period.
- New Adjoint Will compute the adjoint, X*(t)
of the periodic
solution and normalize it so that
< X'(t),X*(t) > =1 .
The
adjoint satisfies
d/dt(X*(t)) = -A^T(t) X*(t) where
A(t) is the derivative of F(X) evaluate at the
periodic.
- Create H function will then do the required averaging
Input the influence of one oscillator on the other. The oscillator 1 in
the above equations has the name of the usual variables while the
"second" oscillator has the same name but with a prime. Thus for
example, if
G(X1,X2) = d (X2-X1)
then type in
d *(X1'-X1)
The H function, its odd part and its even part will be stored
in the second through fourth columns of the data browser. You can plot
them if you want to look at them.
- Adjoint,Hfunction,Orbit just exchange the various
function in and out of the data browser once they are computed.
PhasePlane Menu
This is a bunch of routines that are useful for the analysis of
two-dimensional vector fields. This is what the PP in WINPP
stands for. This works only if the view is two-dimensional and both
axes are differential equations variables. The equations are:
dx/dt=f(x,y)
dy/dt=g(x,y)
If the dimension of your system is higher than 2, the variables not
represented in the view are held constant.
- Nullclines draws the null curves, f(x,y)=0,g(x,y)=0
by breaking the plane into a mesh the size of which is defined in
the Grid parameters item.
- Direct fld draws a direction field of little lines,
indicating the flow.
- Flow Breaks the plane into a grid determined in the
Grid parameters and uses these points as initial condirions. It
integrates both backward and forward in time to produce a flow.
- redraw Nullclines redraws them once they have been
computed.
- Color phaseplane uses the color code you
defined in the colorize dialog and the same grid
as the
direction field to colorize the plane.
- Grid parameters sets the grid for the nullclines,
direction field, and flow.