README for WINPP

CONTENTS

  1. INTRODUCTION

    1. LOGFILES
    2. SUBWINDOWS
    3. QUICK GUIDE TO THE MOST USED FEATURES
  2. INPUT FILES

    1. INTEGRAL EQUATIONS
    2. OPTIONS
    3. RESERVED WORDS
  3. EXAMPLES

    1. Damped Pendulum
    2. Stochastic equations
    3. Delay-differential equations
    4. Integro-differential equations
    5. Differential-algebraic equations
    6. Boundary value problems
    7. Equations with arrays, PDEs etc
    8. Equations with delta function discontinuities
  4. WINPP Commands

    1. Running the simulation

      Click on Run Go, Run Last, or
      1. Using the mouse to set initial conditions (Run Mouse or Run Mice)
      2. Integrating over a range of parameters or ICs (Run Range)

    2. Finding Equilibria (Run Equilibria)
    3. Boundary value solving (Run BVP)
    4. Plotting and browsing

      1. Adjusting the view, setting 2D and 3D parameters (Graphics View)
      2. Printing (File Print)
      3. Plotting a variable versus time(Graphics XvsT)
      4. Multiple variables in the same plot and editing plot options (Graphics Add Curve)
      5. Create a new plot window (Graphics New Window)
      6. Keeping curves over different runs (Graphics Keep)

  • Plotting arrays (Graphics Array)
  • Transposing data (Numerics Transpose)
  • Animation (Graphics Toon)
    1. Example

  • File Menu

    1. Saving the state of WINPP (File Write)
    2. Recalling the state (File Read)
  • Changing right-hand side (Edit Edit RHS)
  • Copy to clipboard (Edit Copy) or Ctrl-C
  • Numerical parameters

    1. Time step, total time, method, max delay, output, transient,tolerances (Numerics Int) or (Ctrl-I)
    2. Delayed plots (Numerics Ruelle)
    3. Variables on a circle (Numerics Torus)
    4. Newton control, BVP control, DAE control etc (Numerics Control)
    5. Color coding trajectories (Numerics Colorize)
    6. Poincare maps, maxima, interspike intervals (Numerics Poincare)
    7. Transposing array data (Numerics Transpose)
    8. Averaging
  • Phase plane analysis

    1. Nullclines (Phaseplane Nullclines)
    2. Direction fields (Phaseplane Direct Fld)
    3. Flow diagrams (Phaseplane Flow)
    4. Colorized phase-plane (Phaseplane Color)
    5. Grid size for phase planes (Phaseplane Grid)


    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:

    1. ordinary differential equations (ODEs)
    2. delay differential equations (DDEs)
    3. differential-algebraic equations (DAEs)
    4. Volterra integro-differential equations (IDEs)
    5. discrete dynamical systems (MAPs)
    6. 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

    1. To integrate the equations, click on Run Go. Or just type Ctrl G.
    2. 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.
    3. 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.
    4. A quick way to plot a quantity versus time is to click Graphics XvsT or Ctrl-X.
    5. Ctrl-E erases the window and Ctrl-R redraws it.
    6. File Exit will get you out as will Ctrl-Q.
    7. Ctrl-L will use the endpoint of the most recent integration as the starting point of a new one.
    8. To change integration parameters, click on Numerics Int. Pars. or type Ctrl I.
    9. Click on Edit Copy to copy the picture to the clipboard and then paste it into your own documents or click Ctrl-C.
    10. 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.
    11. 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.

    The remaining options can be set from within the program. They are

    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:


    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:

    1. 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.
    2. I have defined some auxiliary variables that can be viewed by the user.
    3. 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:

    1. 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.
    2. WINPP is case-insensitive so EXP and Exp and eXp are all the exponential function.
    3. Floating point exponents are of the form #e+#, #e-#, #e#, where # is any number.
    4. The functions am(v) and bm(v) are defined by a simple obvious statement.
    5. 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:

    1. 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.
    2. Within a right-hand side the function delay(x,tau) produces, x(t-tau).
    3. 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:

    1. 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)
    2. yplot=f2 tell WINPP to plot f2(t).
    3. 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:

    1. 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.
    2. 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:

    1. 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:

    1. 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:

    1. 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.
    2. 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.
    3. 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:

    1. 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.
    2. It is written as a system of first order equations.
    3. 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:

    1. [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 .
    2. This is a crude discretization

    3. 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.

    4. 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.


    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:


    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:


    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:

    The other buttons in the animation window are:

    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:


    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:


    Edit Menu

    There are two items:


    Numerics Menu

    This menu has many numerically associated commands, most importantly the ones that control the main integration routines.

    The commands are:

    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.


    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.