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:
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
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.
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.
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.
# 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
The general integral equation:
becomes
u = f(t) + int{K(t,t',u)}
The convolution equation:
would be written as:
v(t) = exp(-t) + int{exp(-t^2)#v}If one wants to solve, say,
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}
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
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:
Here I present some example ODE files that you are welcome to try out.
# 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
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
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
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
This is a simple DAE. We want to solve:
y+x=1
# dae_ex1.ode x'=y 0= x+y-1 x(0)=0 solv y=1 aux yy=y done
Here is a more complex example in which the derivative is implicitly defined.
# 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
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
Here we solve a nonlinear boundary value problem:
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
Here is a very complicated boundary value problem.
Find a value of omega solving:
d (k'+k/t+2ka'/a) + (omega + q a^2) =0
k(0)=0
k(1)=0
a'(1)=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
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+x4This makes it easy to solve discretized PDEs. For example, the simple Fisher equation:
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
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.
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:
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
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
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
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.
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.
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:
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:
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:
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 .
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.
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
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:
$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.
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
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 doneRun 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] endThis 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 .
The File menu controls IO stuff.
The commands are:
There are two items:
This menu has many numerically associated commands, most importantly the ones that control the main integration routines.
The commands are:
The Pick method box lets you choose from many integrators.
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.
There are 3 types of maps:
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.
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:
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.
The adjoint satisfies
where A(t) is the derivative of F(X) evaluate at the periodic.
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
then type in
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.
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:
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.