% This is a LaTeX file. To use it you must first process it with % LaTeX. Then you can render it with whatever printer you want. If % you don't have TeX, then you should. It is free and anyone doing % mathematics should have it. \documentstyle{article} \newcommand{\beq}{\begin{equation}} \newcommand{\eeq}{\end{equation}} \newcommand{\beqa}{\begin{eqnarray}} \newcommand{\eeqa}{\end{eqnarray}} \newcommand{\bvb}{\begin{verbatim}} \newcommand{\evb}{\end{verbatim}} \newcommand{\beqann}{\begin{eqnarray*}} \newcommand{\eeqann}{\end{eqnarray*}} \newcommand{\nn}{\mbox{${\nonumber}$}} \newcommand{\labeq}[1]{\label{eq:#1}} \newcommand{\refeq}[1]{(\ref{eq:#1})} \newcommand{\tc}[1]{\addcontentsline{toc}{subsection}{#1}} \newcommand{\tcc}[1]{\addcontentsline{toc}{subsubsection}{#1}} \title{XPP1.6 -- the differential equations tool} \author{Bard Ermentrout} \begin{document} \maketitle \tableofcontents \newpage \section{Introduction} XPP is a tool for solving differential equations, difference equations, delay equations,functional equations, boundary value problems, and stochastic equations. It evolved from a chapter written by John Rinzel and myself on the qualitative theory of nerve membranes and eventually became a commercial product for MSDOS computers called PHASEPLANE. It is now available as a program running under X11 and UNIX. It is free of charge and the source code is freely distributed so if you don't like the way it works rewrite it. The code brings together a number of useful algorithms and is extremely portable. All the graphics and interface are written completely in Xlib which explains the somewhat idiosyncratic and primitive widgets. Indeed, I began this using the OpenLook widget set and then found out that no one but SUN used it so that most folks would not easily be able to use it. The SPARC that I purchased did not have the Athena widgets so I never bothered to learn them. After finally getting them, it was too late to learn something new so I have stuck with the primitive widgets that I wrote. Its ugly but it works. There are a number of other such programs available, but they all seem to require that your problems be compiled before using them. XPP does not; I have devised a simple and fairly fast formula compiler that is based on the idea of the inner interpretor used in the language FORTH (which remains my first love as far as language is concerned) Fear not, the differentiual equations and boundary conditions and other formulae are written in usual algebraic notation. XPP has been successfully compiled on a SPARC II under OpenLook, a SPARC 1.5 running generic X, a NeXT running X11R4, a DEC 5000, a PC using Linux, and an HP 730. I cannot vouch for other platforms particularly the rather picky IBM RS6000. Building XPP requires only the standard C compiler, Xlib, and the Curses library. The latter is only needed for reading in files and can be modified to work with stdlib with a minimum of work. XPP has the capabilities for handling up to 50 differential equations. There are solvers for delay and stiff differential equations a well as some code for boundary value problems. Difference equations are also handled. Up to 10 graphics windows can be visible at once and a variety of color combinations is supported. Primitive PostScript capabilities are supported. Post processing is easy and includes the ability to make histograms, FFTs and applying functions to columns of your data. (If not in this version, then in a later one...) Equilibria and linear stability as well as one-dimensional invariant sets can be computed. Nullclines and flow fields aid in the qualitative understanding of two-dimensional models. Poincare maps and equations on cylinders and tori are also supported. Some useful averaging theory tricks and various methods for dealing with coupled oscillators are included primarily because that is what I do for a living. I will assume that you are well versed in the theory of ordinary differential equations although you need not be to use the program. The basic unit for XPP is a single ASCII file that has the equations, parameters, variables, boundary conditions, and functions for your model.The methods of solving the equations, the graphics and postprocessing are all done within the program using the mouse and various menus and buttons. The impatient user should look at some sample {\tt *.ODE} files instead of actually reading the documentation. In additions, XPP uses a file called {\tt default.opt} that describes initialization and memory options for the program. \addcontentsline{toc}{subsection}{Interface} \begin{center} {\Large Notes on the Interface } \end{center} This interface is crude by most standards and it is rather ugly as well. I had originally thought to do the whole thing using the Openlook Widget Set but found that many of the systems to which I wanted to port XPP did not have the set. The present design is completely portable. The text editing has been improved over the original version of XPP. One can use the {\tt Home,End} keys to move the text cursor to the beginning and end of the line. The left and right arrow keys let you move back and forth. You can now insert or delet text at any point. I will eventually make it somewhat EMACS compatible. Almost every command has a keyboard shortcut. These are given below. \tc{Disclaimer} \begin{center} \Large Disclaimer \\ \end{center} XPP is distributed as is. The author makes no claims as to the performance of the program. Anyone is allowed to modify and distribute XPP as long as the original code is also made available. \tc{Acknowledgments} \begin{center} Acknowledgements \end{center} The porting of XPP to UNIX and X has been partially supported by NSF in the sense that they gave me the SPARC 2. MSRI in Berkeley gave me use of a lovely office facing the bay and a workstation on which the majority of the port was made in the Summer of 1992. The original version of PhasePlane benefitted by countless colleagues whose endless requests for features and whose constant use led to the version that you see now. Of these I would like to specifically cite John Rinzel and Artie Sherman at the NIH for being the main guinea pigs for the earliest versions. Please let me know of any bugs or other stuff that you'd like to see incorporated into XPP. I will usually fix them quickly. My EMAIL address is bard@mthbard.math.pitt.edu. \section{ODE Files} ODE files are ASCII readable files that the XPP parser reads to create machine useable code. They consist of combinations of lines each starting with a key letter. Lines cannot presently be continued onto the next line but the length of the lines is 256 characters. \bigskip \tc{Example} I will start with a very simple example: \beq x'=x(a+f(t)-x) \eeq the periodically varying logistic equation with $x(0)=.1,$ $f(t)=a\sin(\omega t),$ $a=1,b=.05,\omega=1.2$ The file is: \bvb 1 v x=.1 p a=1,b=.05,omega=1.2 u f 1 b*sin(omega*arg1) o x*(a+f(t)-x) d \end{verbatim} \medskip The first line says there is 1 variable. The name of the variable and a starting initial value is given next. Parameters are defined in the third line. A user-defined function is then created. The right-hand side of the equation is next and the last line tells XPP that the equation is done. This can be explored by typing: \bvb xpp \end{verbatim} at the terminal prompt. Here {\tt } is the name of the above file. \bigskip \tc{ODE File format} More generally, the ODE files have the following format : \bvb change variable =, ... ... fixed ,... ... markov ... parameter =, ... ... table ... table % ... user ... kernel ... ode ... integral ... ode (for fixed and auxilliary quantities) ... bdry ... r {t01} {t02} ... {t0k-1} {t10} ... ... {tk-1,0} ... # Comments done \end{verbatim} Because the program has a long history and there are many old versions around, the syntax above has been maintained even though greatly improved syntax is easy to create (see for example my code for PDEs called XTC.) Thus, there are some confusing restrictions that are necessary in order to maintain backward compatibility with the hundreds of ODE files that I have already written. The typical ODE file contains some or all of the above types of lines. {\em The first line of the file contain a positive number.} This tells XPP how much space to allocate for the value of your variables and auxilliary functions that you might like to maintain. This can be up to 50 columns. This number is the dimension of your system plus the number of extra quantities you record. (XPP cannot dynamically allocate more columns so that you should allocate extra columns if you want to use them. For each extra column, a line of the form \bvb o 0 \end{verbatim} should follow your main stuff before the last line. This fills the columns with zeros.) Each of the remaining lines contains one of the above key words at its beginning. Since only the first letter is used, you need only use the first letter. \medskip XPP uses a bunch of defaults when it is started up by looking for a file called ``default.opt.'' If it cannot find it, it uses internal options. Alternatively, you can tell XPP the name of the options file you want to use. The description of these files is below. The format for such a statement is: \bvb c \end{verbatim} which loads the options file specified in {\tt}. \medskip Variables are the quantities you wish to integrate in time. They (as can all names in XPP) can have up to 9 letters each. XPP is case insensitive. Any combination of letters and numbers is valid as is the underscore, ``\_''. The order of statements is not too important however you must not refer to names that have yet to be defined as XPP compiles on the run. Generally, one declares in this order: change options file, variables, fixed variables, parameters, kernels, tables, user functions, right-hand sides, boundary conditions. The right-hand sides are started with either ``o'' or with ``i'' depending on whether or not you are solving an differential equation or a Volterra equation. \medskip \tcc{Variables, Fixed, Markov} Variables can be declared one to a line or all on the same line separated by commas or spaces. Initial conditions are optional since they can be set within XPP; there can be {\em no spaces} separating the ``='' sign, the variable, and the number. If no initial data are given, the variable is set to zero. \medskip Fixed quantities are declared on one line or on several with no initialization; this is done later via the right-hand-side operation starting with ``o.'' \medskip Markov variables are finite state quantities that randomly flip from one state to another according to a given transition probability table that is given in the file. They are treated like variables in that they have initial conditions and are accessible to the user. The only requirement is that {\em they must be declared after all variables and fixed quantities are declared.} If you do not follow this convention then the program will definitely screw up. Each Markov variable has its own line. The name of the variable, the initial value, and the number of states are included on the line. The states are 0,1, and so on up to $k-1$ where $k$ is the number of states. \medskip \tcc{Parameters} Parameters are like variables and can be on one or several lines. The default for uninitialized parameters is 0. \medskip \tcc{Tables} The ``table'' declaration allows you to declare a function called {\tt } that reads in values from the file, {\tt } or as a function of one variable over some interval and then is used in your program as a function of 1 variable interpolated from the tabulated values. The values of this table are assumed to be equally spaced and the file is an ASCII file with the format: \bvb y1 y2 . . . yn \end{verbatim} Thus, $f(xlo)=y1$ and $f(xhi)=yn.$ These tables can also be read in from within XPP but a valid table must be given to start the program. This table can be arbitrarily long (as memory permits) and thus you can use experimental data as inputs to differential equations or even sketch curves and use that as your nonlinearity. Ypu can directly input tabulated functions as well to speed up computations with complicated functions. In this case, after the table name put a parenthesis symbol followed by the number of points, the minimum argument and the maximum argument and then the function. This should be written as a function of ``t''. Thau, the statement \begin{verbatim} table f % 501 -10 10 tanh(t) \end{verbatim} will produce a table of the hyperbolic tangent function from -10 to 10 consisting of 501 points. \medskip \tcc{User defined functions} User defined functions have the following form: \bvb user \end{verbatim} {\tt } is the name of the function. Next, the number of arguments is given. The maximum number of arguments is 9. Finally the expression for the function must be given. Line breaks are not allowed so the expression must be typed on one line. The arguments are \bvb arg1 arg2 ... arg9 \end{verbatim} and must be used as the formal arguments in a function. \medskip Integral equations are constructed in several different ways, the ``kernel'' declaration is one such way. This and the ``integral'' construct are described in a section devoted to Volterra equations. \tcc{Boundary conditions} Boundary conditions can be placed anywhere in the file or ignored altogether if you dont plan on solving boundary value problems. They have the form: \bvb b \end{verbatim} where ``b'' stands for boundary condition. The expression is one involving your variables and which will be set to zero. In order to distinguish left and right boundary conditions, the following notation is used. For the values of the variables at the left, use the symbol for the variable. For the values at the right end of the interval, use the symbol for the variable appended by a single quote, '. Thus to specify the boundary conditions $x(0)=1$, $y(1)=2$, the following is used: \bvb b x-1 b y'-2 \end{verbatim} Periodic boundary conditions would be written as: \bvb b x-x' b y-y' \end{verbatim} Note that it is not necessary to specify the BCs at this point. They can easily be specified within the program. Also note that the notation I have used allows the specification of mixed boundary conditions such as periodic BCs. (XPP has some additional special commands for periodic boundary conditions that enable the user to specify fewer equations than usual.) The number of boundary conditions must match the number of variables in your problem. \medskip \tcc{Right-hand sides} The most important and most confusing part of XPP is the declaration of right-hand sides. These are declared by starting the line with the letter ``o'' (historically for ODE) and following with some expression. (For Volterra equations, you will use the letter ``i''; see the description of this below in the section on Volterra models.) The order of these is very important. Basically the ``gross order'' is: RHS for variables, RHS for fixed quantities, RHS for auxilliary quantities. Let me first distinguish between these three. Variables refer to the quantities whose evolution determines the differential equation. Fixed quantities are temporary combinations of variables (for example a very complicated function of the variables) that can be used repeatedly in the right-hand sides determining the differential equation. These are like local quantities defined in a FORTRAN program. Auxilliary variables or quantities are additional combinations of the variables and fixed quantities whose values you want to keep during the simulation; things like energy or mean values. Within each of these three classes, the ordering is also important. The first RHS you define must be that of the first dependent variable. Continue defining each successive RHS until you have as many as there are variables. Then define the RHS for each of the fixed variables in the order in which they were declared. Finally, define the auxilliary variables which have no internal names but are called ``aux1,'' ``aux2,'' etc in the program. The examples shown below should clarify these rules a bit. The right-hand-sides of the ODEs can depend on any of the fixed variables as can the auxilliary quantities. A fixed variable can depend on any parameter or variable as well as a previously defined fixed variable. Thus, ordering is important. Just remember that all of these start with the command ``o'' and are in the order: \begin{verbatim} RHS1 RHS2 ... FIXED1 FIXED2 ... AUX1 AUX2 ... \end{verbatim} \medskip \tcc{Markov probabilities} The letter 'r' is used to tell XPP that you are ready to define a tRansition table for a Markov process whose name follows. (The letter 'r' has no meaning and was used since it is the first free letter I happened upon.) A finite state Markov variable with $k$ states must have associated with it a $k\times k$ matrix, $T_{ij}$ which contains the probability of going from state $i$ to state $j$ per unit of time. Thus the effective transition probability is {\tt DeltaT} times the one in the table; the large is {\tt DeltaT} the higher the probability. Since the probabilities must add to 1 in any row, the diagonal terms are automatic. If there are $k$ states to the variable, then there must be $k$ rows following the declaration of the transition matrix. Each row contains $k$ entries delimited by curly brackets ans separated by spaces; for example with $k=2$ the variable {\tt z} would have its matrix defined by': \begin{verbatim} r z {0} {P01} {P10} {0} \end{verbatim} where {\tt P01, P10} are any algebraic expression or number involving the parameters and variables of the XPP file. At each output time step, the probailities are computed, multiplied by the timestep, and a random number is chosen. If it in the appropriate range, then the transition will be made. Transitions of several such variables are made in parallel and then each is updates. The example below will help. \medskip The last line contains the single letter `d' telling the formula compiler that you are done. You should be aware of the following keywords that should not be used in your ODE files for anything other than their meaning here. \tc{Reserved words} \begin{center} \Large Reserved words \\ \normalsize \bvb 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 max min normal besselj bessely arg1 ... arg9 aux1 ... @ $ + - / * ^ ** \end{verbatim} \end{center} These are mainly self-explanatory. The nonobvious ones are: \begin{itemize} \item {\tt heav(arg1)} the step function, zero if {\tt arg1<0} and 1 otherwise. \item {\tt sign(arg)} which is the sign of the argument (zero has sign 0) \item {\tt ran(arg)} produces a uniformly distributed random number between 0 and {\tt arg.} \item {\tt besselj, bessely } take two arguments, $n,x$ and return respectively, $J_n(x)$ and $Y_n(x),$ the Bessel functions. \item {\tt normal(arg1,arg2)} produces a normally distributed random number with mean {\tt arg1} and variance {\tt arg2}. \item {\tt max(arg1,arg2)} produces the maximum of the two arguments and {\tt min} is the minimum of them. \item {\tt if()then()else()} evaluates {\tt } If it is nonzero it evaluates to {\tt } otherwise it is { \tt }. E.g. {\tt if(x>1)then(ln(x))else(x-1)} will lead to {\tt ln(2)} if {\tt x=2} and { \tt -1 if x=0.} \item {\tt delay(,)} returns variable {\tt } delayed by the result of evaluating {\tt }. 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.) \item {\tt ceil(arg),flr(arg)} are the integer parts of{\tt } returning the smallest integer greater than and the largest integer less than {\tt }. \item {\tt t } is the current time in the integration of the differential equation. \item {\tt pi} is $\pi.$ \item {\tt arg1, ..., arg9} are the formal arguments for functions and {\tt auxn} where {\tt n } is a positive integer are the names for auxilliary variables. \end{itemize} \section{ Examples} Nothing helps one understand how to use a program better than lots of examples. I offer here a few files which you can explore. \tc{Morris-Lecar} \begin{center} Morris-Lecar Equations \end{center} The Morris-Lecar equations arise as a simplification of a model for barnacle muscle oscillations. They have the form: \beqann C\frac{dV}{dt} &=& g_L(V_L-V)+g_Kw(V_k-V)+g_{Ca}m_\infty(V)(V_{Ca}-V)+I \\ \frac{dw}{dt} &=& \phi \lambda_w(V)(w_\infty(V)-w) \eeqann where \beqann m_\infty(V) &=& .5(1+\tanh((V-V_1)/V_2)) \\ w_\infty(V) &=& .5(1+\tanh((V-V_3)/V_4)) \\ \lambda_w &=& \cosh((V-V3)/(2V_4)) \eeqann The ODE file for these equations (dimensionless,so that $C=1$ and the conductances are relative to $g_{Ca}=1$ and the reversal potentials to $V_{Ca}=1$) are: \begin{verbatim} 2 # Declare the variables v v,w # Declare the parameters p gk=.5,gca=1,gk=2 p vk=-.7,vl=-.5,vca=1 p v1=.01,v2=.145,v3=.1,v4=.15 p i=.2,phi=.333 # Define some functions u minf 1 .5*(1+tanh((arg1-v1)/v2)) u winf 1 .5*(1+tanh((arg1-v3)/v4)) u lamw 1 cosh((v-v3)/(2*v4)) # define the right-hand sides ov gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v)+i ow phi*lamw(v)*(winf(V)-w) # Now we are done d \end{verbatim} This is a very straightforward model and so the equation file is pretty simple. Suppose you want to keep track of the calcium current as an auxilliary variable. Then the file looks like: \begin{verbatim} 3 # Want an extra quantity, ICa, so allocate 3 slots, v,w, and ICa # # Declare the variables v v,w # Declare the parameters p gk=.5,gca=1,gk=2 p vk=-.7,vl=-.5,vca=1 p v1=.01,v2=.145,v3=.1,v4=.15 p i=.2,phi=.333 # Define some functions u minf 1 .5*(1+tanh((arg1-v1)/v2)) u winf 1 .5*(1+tanh((arg1-v3)/v4)) u lamw 1 cosh((v-v3)/(2*v4)) # define the right-hand sides ov gl*(vl-v)+gk*w*(vk-v)+gca*minf(v)*(vca-v)+i ow phi*lamw(v)*(winf(V)-w) # # Now we simply keep ICa: o -gca*minf(v)*(vca-v) # Now we are done d \end{verbatim} This can be accessed in XPP by requesting ``AUX1'' as it is the first auxilliary variable. Note that we are wasting computational time since we compute $I_{Ca}$ twice; once as a contribution to the potential change and once as an auxilliary variable. In a FORTRAN or C program, one would compute it as a local variable and use it in both instances. This is the purpose of fixed variables. (For the present problem, the computational overhead is trivial, but for coupled arrays and other things, this can be quite substantial.) The last program uses a fixed variable to reduce the computation: \begin{verbatim} 3 # Want an extra quantity, ICa, so allocate 3 slots, v,w, and ICa # # Declare the variables v v,w # # Declare a fixed quantity ICA, the calcium current # f ica # Declare the parameters p gk=.5,gca=1,gk=2 p vk=-.7,vl=-.5,vca=1 p v1=.01,v2=.145,v3=.1,v4=.15 p i=.2,phi=.333 # Define some functions u minf 1 .5*(1+tanh((arg1-v1)/v2)) u winf 1 .5*(1+tanh((arg1-v3)/v4)) u lamw 1 cosh((v-v3)/(2*v4)) # define the right-hand sides # Note that we include ICa as a computed quantity # ov gl*(vl-v)+gk*w*(vk-v)+i-ica ow phi*lamw(v)*(winf(V)-w) # # Now we need to define ICa. Remember first RHS of ODEs, then # fixed variables, then auxiliary # oica -gca*minf(v)*(vca-v) # Now we simply keep ICa: o ica # Now we are done d \end{verbatim} The calcium current is computed once instead of twice. This is why ``fixed'' variables are useful. \tc{Cable model} \begin{center} A linear cable equation with boundary conditions \end{center} In studying a dendrite, one is often interested in the steady-state voltage distribution. Consider a case where the dendrite is held at $V=V_0$ at $x=0$ and has a leaky boundary condition at $x=1.$ Then the equations are: \beqann \lambda^2\frac{d^2V}{dx^2} &=& V(x) \\ V(0)&=& V0 \\ a\frac{dV(1)}{dx} + b V(1) &=& 0 \eeqann When $a=0,b\ne0$ the voltage at $x=1$ is held at 0. When $a\ne0,b=0$ there is no leak from the cable and the conditions are for sealed end. Since this is a second order equastion and XPP can only handle first order, we write it as a system of two first order equations in the file which is: \begin{verbatim} 2 # define the variable, v and its derivative vx v v=1,vx=0 # define the 4 parameters, a,b,v0,lambda^2 p a=1,b=0,v0=1,lam2=1 # now do the right-hand sides o vx o v/lam2 # # and finally, boundary conditions # First we want V(0)-V0=0 b v-v0 # # We also want aV(1)+bVX(1)=0 b a*v'+b*vx' # Note that the primes tell XPP to evaluate at the right endpoint d \end{verbatim} Note that I have initialized $V$ to be 1 which is the appropriate value to agree with the first boundary condition. \tc{Delayed feedback}\begin{center} The delayed inhibitory feedback net \end{center} A simple way to get oscillatory behavior is to introduce delayed inhibition into a neural network. The equation is: \[ \frac{dx(t)}{dt} = -x(t) + f(ax(t)-bx(t-\tau)+P) \] where $f(u)=1/(1+\exp(-u))$ and $a,b,\tau$ are nonnegative parameters. Here $p$ is an input. The XPP file is: \begin{verbatim} 1 # declare the scalar x v x=1 # declare all the parameters, initializing the delay to 3 p tau=3,b=4.8,a=4,p=-.8 # define the nonlinearity u f 1 1/(1+exp(-arg1)) # define the right-hand sides; delaying x by tau o -x + f(a*x-b*delay(x,tau)+p) # done d \end{verbatim} Try this example after first going into the numerics menu and changing the maximal delay from 0 to say, 10. Glass and Mackey describe a few other delay systems some of which have extremely complex behavior. Try them out. \tc{Random mutation}\begin{center} A population problem with random mutation \end{center} This example illustrates the use of Markov variables and is due to Tom Kepler. There are two variables, $x_1,x_2$ and a Markov state variable, $z.$ The equations are: \beqann x_1' &=& x_1(1-x_1-x_2) \\ x_2' &=& z(ax_2(1-x_1-x_2)+\epsilon x_1) \eeqann and $z$ switches from 0 to 1 proportionally to $x_1.$ The transition matrix is 0 everywhere except in the 0 to 1 switch where it is $\epsilon x_1.$ So z=1 is an absorbing state. Initial conditions should be $x_1(0)=1.e-4$ or so, $x_2(0)$ the same (this is just for convenience; we really want $x_2(0)=0$ and then have it jump discontinuously to $1.e-4$ or so when $z$ makes its transition, but this shouldn't matter that much) This models the population dynamics of two populations $x_1,x_2$ in competition with each other. $x_2$, initially absent, is a mutant of $x_1.$ The mutation rate $\epsilon$ should be smaller than one. The relative advantage, $a$, should be larger than one. One expects that $x_1$ grows for a while, eventually $z$ makes its transition, $x_2$ begins to grow and eventually overtakes $x_1.$ The equation file for this is \begin{verbatim} 3 v x1=1.e-4,x2=1.e-4 m z 0 2 p eps=.1,a=1 o x1*(1-x1-x2) o z*(a*x2*(1-x1-x2)+eps*x1) r z {0} {eps*x1} {0} {0} d \end{verbatim} The contiunuous variables are declared first, then the Markov variable with 2 states. The parameters are defined. The two differential equation are defined, and finally the transition matrix for $z.$ \medskip XPP comes with some other examples that I urge you to look at. A Volterra example is shown below. \section{Using xpp} XPP runs from the command line by typing \bvb xpp \end{verbatim} where {\tt } is an ODE file as described above. This argument is optional. If you invoke XPP without the argument, you are asked whether you want to read or create a file. Choosing {\tt (r)} to read the file prompts you for an ODE file name. You can get a directory by pressing {\tt }. If you choose to create a new file, you must write a set of commands such as above and then save the file. It will bounce you out if you make any mistakes. It is not the recommended way to do it -- use an editor and read in the file. XPP will look at the extension of a file and if it is {\tt .dif } the program will assume that the model is a difference equation and will treat it as a map. The solution, equilibria, stability and listing of difference equations is treated differently from differential equations. The default plot style is also plotting of points instead of lines as for the continuous time case. \tc{Windows} \tcc{Main Window} Once you get a file loaded and it is valid, the program and all of the windows will pop up. There is the main window with the graphics and menus. This has 4 regions on it: the graphics window, the menu window, the command window across the top and the message window at the bottom. Initial conditions and parameters etc are input here as well as text labels. The message window gives the position of the pointer. It used to give messages but this action has been supplanted by the message widget and the error widget. If you click the second or third mouse button in the graphics window, the coordinates of the pointer are shown in the message window. The left most region is the menu. Pressing the mouse on a command is the same as typing the hot letter for the command and will execute it. The pointer is a hand in the menu window. \tcc{Initial data,delays,parameters} In addition to the main window, there is the initial data window which contains the names of the variables and the current initial data. The delay window is similar containing the initial data for delay equations. The parameter window contains the list of all parameters and their current values. The boundary condition window is similar with the list of boundary conditions given. These 4 windows constitute the main user editable windows. All of the quantities and formulae contained within these can be edited within the program. Clicking on any item in these windows will place the item on the command line where you can edit it. Pressing {\tt Back Space} moves back one character. Pressing {\tt Esc} exits without changing while pressing{\tt Return} accepts the current value. Clicking {\tt DONE} in the relevant window accepts the present value and exits. If you have chosen ICs, BCs, or Delay ICs, then after accepting a value, the program prompts you for the next one until until you have finished or exited the process. If you are changing parameters, then after accepting a parameter, you will be prompted for the name of anew prameter. Either type it in or click on the name in the parameter window. Click on {\tt DONE}, type {\tt Esc} or type {\tt Return} twice to exit. \tcc{Equation window} Another window contains the listing of the differential equations. If the file is a difference equation, the listing shows this. These equations can be scrolled up and down by pressing the mouse on the {\tt } buttons or using the cursor keys. Resizing lets you see more and closing iconifies the window. Once you are in XPP you cannot alter the differential equations. Thus , to solve a different problem you must exit and start again. I hope to change this in the future. All of these windows can be closed with no effect on the performance of the program. \tcc{Browse window} The other major window contains listings of the numerical values of the solutions to the differential equations. This browser window allows you to postprocess data and save and load numbers to the disk. Pressing the cursor keys, the {\tt HOME, END, PgUp,PgDn} keys scrolls you through tthe data either up and down or left and right. Resizing the window lets you see more data. There are many other keys and you can click on the buttons or use the keyboard if the mouse pointer is in the Browser. You can load or save data files by choosing the {\tt Load} and {\tt Sav} e keys. You will be asked for a filename which you should enter. Choosing {\tt Get} will load the top row of numbers into the initial data. {\tt Find} will find the data row for which the specified variable is closest to the specified value. {\tt Replace} allows you to replace a column by some formula involving the column itself and all other columns. You are asked for the column to replace. Type in one of the names across the top row. Then you are asked for the formula. Type it in and the entire column will be replaced. For example if the names of your variables are {\tt X} and { \tt Y} and you choose {\tt X } as the column to replace and type \bvb X^2 + T*Y \end{verbatim} then the column {\tt X } will be replaced by the function at each row. If the first character in the formula is the {\tt \& } symbol, an integral of the formula will be placed in the specified column. If the formula is anything of the form : \bvb @ \end{verbatim} the numerical derivative of the {\tt } will be computed ad placed in the column. Choosing {\tt Unreplace} will restore the most recent column. Each time you integrate the equations or solve a boundary value problem or look at a range of equilibria, the numbers are replaced by the most recent. You must save the data to a file to use it again. Pressing {\tt First} and {\tt last} will mark the data so that you can save a portion or restore a portion. {\tt Home} takes you to the beginning and {\tt End} to the end of the list. Other windows are opened on occasion and can remain open throughout the session. The usual methods for resizing windows works on the graphics and listing and browser windows. When there are multiple graphics windows, only one is active at a given time and only one can be drawn to. The window is made active by clicking the left button when inside the window. A little square appears in the upper left corner of the window to indicae its activity. The variables on the axes are given in the title of the window and are Z-Axis vs Y-Axis vs X-Axis. One other window appears permanently, when neeeded, the Equilibrium window that gives values and stability of the most recently computed equilibrium. The main commands are on the first menu and the numerics menu. Communication to the program is via the keyboard and the mouse and through parameter boxes that are displayed. Almost all commands have keyboard shortcuts so that if you know the original DOS version, you can use the X version with little help. When you are prompted for floating point numbers at the command line, then you can insert either numbers or mathematical expressions provided the first character is the {\tt \%} symbol. Thus if you want to input {\tt sin(1.5)} as a floating point number, you would write \bvb %sin(1.5) \end{verbatim} at the command line. \tc{Changing ICs and parameters} \begin{center} \Large Changing initial data and parameters \\ \end{center} Click on the variable you want to change in the IC window. The current value and the name will appear in the command window. Type {\tt Delete} to completely start over or {\tt Back Space} to move back one character at a time. Press {\tt Esc} to abort. Press {\tt Return} and the value will be accepted and the next IC will come up. If you are done, press {\tt } in the IC window. Otherwise XPP will cycle through all of the variables. If you start the value of the initial condition with a {\tt \%} symbol, then you can put a formula in which will be evaluated. For example {\tt \%sqrt(pi) } will be evaluated to 1.772... . The same technique works for delay equation data except that you put in formulae that depend on time and this will create a list of initial data for $t } in the parameter window. You can press the letter {\tt p} to change parameters or simply click on the parameter command in the main menu. Typing {\tt default } will load all the parameters with their values when the program first fires up. \section{The main commands} All commands can be invoked by typing the hot key for that command (capitalized on the menu) or clicking on the menu with the mouse. Usually most commands can be aborted by pressing the {\tt Esc} key. Once one of these is chosen, the program begins to calculate and draw the trajectories. If you want to stop prematurely, press the {\tt Esc} key and the integration will stop. \begin{description} \tc{ICs}\item[(I)nitial conds] This invokes a short list of options for integrating the differential equations. The choices are: \begin{description} \item[(G)o] which uses the initial data in the IC window and the current numerics parameters to solve the equation. The output is drawn in the current selected graphics window and the data are saved for later use. The solution continues until either the user aborts by pressing {\tt Esc}, the integration is complete, or storage runs out. \item[(N)ew] This prompts you at the command line for each initial condition. Press {\tt Return} to accept the value presented. \item[(L)ast] This uses the end result of the most recent integration as the starting point of the curret integration. \item[(S)hift] This is like Last except that the stating time is shifted to the current integration time. This is irrelevant for autonomous systems but is useful for nonautonomous ODEs. \item[(O)ld] This uses the most recent initial data as the current initial data. It is essentially the same as Go. \item[(R)ange] This lets you integrate multiple times with the results shown in the graphics window. Pressing this option produces a new window with several boxes to fill in. First choose the quantity you want to range over. It can be a parameter or a variable. The integrator will be called and this quantity will be changed at the beginning of each integration. Then choose the starting and ending value and the number of steps. The option Reset storage only stores the last integration. If you choose not to reset, each integraion is appended to storage. Most likely, storage will be exceeded and the integration will overwrite or stop. The option to use last initial conditions will automatically use the final result of the previous integration as initial dat for the next integration. Otherwise, the current ICs will be used at each step (except of course for the variable through wich you are ranging.) When you are happy with the parameters, simply press the OK button. Otherwise, press the Cancel button to abort. Assuming that you have accepted, the program will compute the trajectories and plot them storing none of them or all of them. If you press {\tt Esc} it will abort the current trajectory and move on to the next. Pressing the {\tt / } key will abort the whole process. \item[(M)ouse] allows you to specify the values with the mouse. Click at the desired spot. \end{description} \tc{Continue}\item[(C)ontinue] This allows you to continue integrating appending the data to the current curve. Type in the new ending time. \tc{Nullclines}\item[(N)ullclines] This option allows you to draw the nullclines of systems. They are most useful for two-dimensional models, but XPP lets you draw them for any model. The constraints are the same as in the direction fields option above. The menu has 4 items. \begin{description} \item[(N)ew] draws a new set of nullclines. \item[(R)estore] restores the most recently computed set. \item[(A)uto] turns on a flag that makes XPP redraw them every time it is necessary because some other window obscured them. \item[(M)anual] turns this flag off so that you must restore them manually. The X-axis nullcline is blue and the Y-axis nullcline is red. \end{description} \tc{Direction Field/Flow}\item[(D)irection Field/Flow] This option is best used for two-dimensional systems however, it can be applied to any system. The current graphics view must be a two-d plot in which both variables are different and neithet is the time variable, T. There are two items. \begin{description} \item[(D)irection fields] Choosing the direcion field option will prompt you for a grid size. The two dimensional plane is broken into a grid of the size specified and lines are drwan at each point specifying the direction of the flow at that point. If the system is more than two-dimensional, the other variables will be held at the values in the initial conditions window. \item[(F)low] Choosing the flow, you will be prompted for a grid size and trajectories started at each point on the grid will be integrated according to the numerical parameters. Any given trajectory can be aborted by pressing {\tt Esc} and the whole process stopped by pressing {\tt /}. The remaining variables if in more than two-dimensions are initialized with the values in the IC window. \end{description} \tc{Window}\item[(W)indow] allows you to rewindow the current graph. Pressing this presents another menu with the choices: \begin{description} \item[(W)indow] A parameter box pops up prompting you for the values. Press OK or CANCEL when done. \item[(Z)oom in] Use the mouse to expand a region by clicking, dragging and releasing. The view in the rectangle will be expanded to the whole window. \item[Zoom (O)ut] As above but the whole window will be shrunk into the rectangle. \item[(F)it] The most common command will automatically fit the window so the entire curve is contained within it. For three-D stuff the window data will be scaled to fit into a cube and the cube scaled to fit in the window. Use this often. \end{description} \tc{Phase Space}\item[ph(A)se space] XPP allows for periodic domains so that you can solve equations on a torus or cylinder. You will be prompted to make (A)ll variables periodic, (N)o varibales periodic or (C)hoose the ones you want. You will be asked for the period which is the same for all periodic variables (if they must be different, rescale them) Choose them by name or by clicking the appropriate names from the list presented to you. An {\tt X} will appear next to the selected ones. Clicking toggles the {\tt X}. Type {\tt Esc} when done or CANCEL or DONE. XPP mods your variables by this period and is smart enough when plotting to not join the two ends. \tc{Kinescope}\item[(K)inescope] This allows you to capture the bitmap of the active window and play it back. It works poorly in X, you cannot save the bitmaps, nor will it run properly across the network. This is not my fault, it is idiosyncratic with X. It is slow also. My PC makes much smoother animations. In any case, another menu pops up with the choices: \begin{description} \item[(C)apture] which takes a snapshot of the currently active window \item[(K)ill] which deletes all the snapshots \item[(A)utoplay] which continuously plays back snapshots until you press the middle mouse button \item[(P)layback] which cycles thru the pictures each time you click the left mouse button and stops if you click the middle. It will beep when done. \end{description} \tc{Graphics}\item[(G)raphic stuff] This induces a popup menu with several choices. \begin{description} \item[(A)dd curve] This lets you add another curve to the picture. A parameter box will appear aking you for the variables on each axis, a color, and line type (only 1 - solid and 0 - dotted are allowed currently) All subsequent integrations and restorations will include the new graph. Up to 10 per window are allowed. \item[(D)elete last] Will remove most recent curve. \item[(A)ll but first] Deletes all curves but first. \item[(E)dit curve] You will be asked for th curve to edit. The first is 0, the second 1, etc. You will get a parameter box like the add curve option. \item[(P)ostscript] This will ask you for a file name and write a postscript representation of the current window. Nullclines, text, and all graphs will be plotted. Color is not supported. \item[(F)reeze] This will create a permanent curve in the window. Usualy, when you reintegrate the equations or load in some new data, the current curve will be replace by the new data. Freeze prevents this. Up to 10 curves can be frozen. It is not implemented yet. Another menu will pop up with the options: \begin{description} \item[(F)reeze current] This freezes the current curve 0. \item[delete (A)ll] This gets rid of all curves. \item[delete (L)ast] This gets rid of the most recently frozen curve. \end{description} \end{description} \item[n(U)merics] This is so important that a section is devoted to it. See below. \tc{File}\item[(F)ile] This brings up a menu with several options. Type {\tt Esc} to abort. \begin{description} \item[(Q)uit] This exits XPP first asking if you are sure. \item[(P)rt info] This prints info about the current simulation to the console. The equations, the numerical parameters, the initial data, boundary conditions and all parameters are printed. \item[(W)rite set] This creates a file with all of the info about the current numerics, etc as well as all of the currently highlighted graphics window. It is not meant to be readable by the user. It in some sense saves the current state of XPP and can be read in later. \item[(R)ead set] This reads a set that you have previously written. The files are very tightly connected to the current ODE file so you should not load a saved file from one equation for a different problem. \item[(S)ave info] This is like {\tt (P)rt info} but saves the info to a file. It is human readable. \item[(C)alculator] This pops up a little window. Type formulae in the command line involving your variables and the results are displayed in the popup. Click on Quit or type {\tt Esc} to exit. \end{description} \tc{Parameters}\item[(P)arameters] This is just like clicking on the parameter window. Type the name of a parameter to change or click on its name in the parameter window. Type {\tt default} to get back the values when you started XPP. \tc{Erase}\item[(E)rase] erases the contents of the ative window, redraws the axes, deletes all text in the window, and sets the redraw flag to Manual. \tc{Multiple windows}\item[(H)alfwindow] This name is a holdover from the DOS version. The option allows you to create and destroy graphics windows. There are several choices. \begin{description} \item[(C)reate] makes a copy of the currently active window and makes itself active. You can change the graphs in this window without affecting the other windows. \item[(K)ill all] Removes all but the main graphics window. \item[(D)estroy] This destroys the currently active window. The main window cannot be destroyed. \item[(A)uto] turns on a flag so that the window will automatically be redrawn when needed. \item[(M)anual] turns off the flag and the user must restore the picture manually. \item[(B)ottom] puts the active window on the bottom. \end{description} \tc{Adding text}\item[(T)ext] allows you to write text to the display. Type in the text and move the mouse to the position and click the button. The text is placed in the window. The text will only be drawn if you have clicked in the active window. Text is remembered until you erase the screen. If you resize or rewindow, it is placed in the appropriate position. \tc{Equilibria}\item[(S)ing pts] This allows you to calculate equilibria for a system. There are three options. \begin{description} \item[(G)o] begins the calculation using the values in the initial data box as a first guess. Newton's method is applied. If a value is found XPP tries to find the eigenvalues and asks you if you want them printed out. If so, they are written to the console. Then if there is a single real positive or real negative eigenvalue, the program asks you if you want the unstable or stable manifolds to be plotted. Answer yes if so and they will be approximated. The calculation will continue until either a variable goes out of bounds or you press {\tt Esc}. If {\tt Esc} is pressed, the other branch is computed. The program continues to find any other invariant sets until it has gotten them all. These are not stored. The program will continue integrating beyond the current numerical parameters until you press {\tt Esc}. Once an equilibrium is computed a window appears with info on the value of the point and its stability. The top of the window tells you the number of complex eigenvalues with positive,{\tt (c+)}, negative {\tt (c-)}, zero {\tt (im)} real parts and the number of real positive {\tt (r+)} and real negative { \tt (r-)} eigenvalues. If the equation is a difference equation, then the symbols correspond the numbers of real or complex eigenvalues outside {\tt (+)} the unit circle or inside { \tt (-)} This window remains and can be iconified. \item[(M)ouse] This is as above but you can specify the initial guess by clicking the mouse. Only the two variables in the two-D window will reflect the mouse values. This is most useful for 2D systems. \item[(R)ange] This allows you to find a set of equilibria over a range of parameters. A parameter box will prompt you for the parameter, starting and ending values, number of steps. Additionally, two other items are requested. Column for stability will record the stability of the equilibrium point in the specified column. If you elect to shoot at each, the invariant manifolds will be drawn for each equilibrium computed. The stability can be read as a decimal number of the form {\tt u.s } where {\tt s} is the number of stable and {\tt u} the number of unstable eigenvalues. So {\tt 2.03} means 3 eigenvalues with negative real parts (or in the unit circle) and 2 with positive real parts (outside the unit circle.) The result of a range calculation is saved in the data array and replaces what ever was there. The value of the parameter is in the time column, the equilibria in the remaining columns and the stability info in whatever column you have specified. {\tt Esc} aborts one step and {\tt /} aborts the whole procedure. \end{description} \tc{Axes}\item[(V)iew axes] This selects one of 9 different types of graphs. The first 5 are two-d and the last 4,three D. If you select a two-D curve, you will be asked for limits as in the window command as well as the variables to place on the axes. Three D is more complicated. You will be asked for the 3 variables for the 3 axes, their max and min values and 4 more numbers, {\tt XLO}, etc. XPP first scales the data to fit into a cube with corners (-1,-1,-1) and (1,1,1). Rotation of this cube is performed and then projected into the two-D window. {\tt XLO}, etc define the scales of this projection and are thus unrelated to the values of your data. Use the {\tt (F)it } option if you don't know whats going on. Numbers appear on the 2D axes but not the 3D, however the axes labels are drawn. \tc{Time plots}\item[(X)i vs t] This chooses a certain 2D view and prompts you for the variable name. The window is automatically fitted and the data plotted. It is a shortcut to choosing a view and windowing it. \tc{Redrawing}\item[(R)estore] redraws the most recent data in the browser in accordance with the graphics parameters ofthe active window. \tc{3D Parameters}\item[(3)d params] This lets you choose rotations of the axes and persepctive planes. Play with this to see. You must be have a 3D view in the active graph to use this. \tc{Boundary values}\item[(B)ndry val] This solves boundary value problems by numerical shooting. There are 4 choices. \begin{description} \item[(S)how] This shows the successive results of the shooting and erases the screen at the end and redraws the last solution. The program uses the currently selected numerical integation method, the current starting point, {\tt T0} as the left end time and {\tt T0+TEND} as the right end. Thus, if the interval of interest is {\tt (2.5,6)} then set {\tt T0=2.5} and {\tt TEND=3.5} in the numerics menu. \item[(N)o show] This is as above but will not show succesive solutions. \item[(R)ange] This allows you to range over a parameter keeping starting or ending values of each of the variables. A window will appear asking you for the parameter, the start, end, and steps. You will also be asked if you want to cycle color which means that the results of each successful solution to the BVP will appear in different colors. Finally the box labeled{\tt side} tells the program whether to save the initial {\tt (0)} or final{\tt (1)} values of the solution. As the program progresses, you will see the current parameter in the info window under the main screen. You can abort the current step by pressing {\tt Esc} and the whole process by pressing {\tt /}. \item[(P)eriodic] Periodic boundary conditions can be solved thru the usual methods, but one then must write an addition equation for the frequency parameter. This option eliminates that need so that a 2-D autonomous system need not be suspended into a 3D one. You will be asked for the name of the adjustable parameter for frequency. You will also be asked for the section variable and section. This is an additional condition that must be satisfied, namely, $x(0)=x_0$ where $x$ is the section variable and $x_0$ is the section. Type {\tt yes} if you want the progress shown. \end{description} \end{description} \section{Numerical parameters} When you click on the {\tt nUmerics} command in the main menu, a new list appears. This is the numerics menu and allows you to set all of the numerical parameters as well as some post-processing. Some of these may not yet be implemented. Press {\tt Esc} or click on the {\tt exit} to get the main menu back. The items on this menu are: \begin{description} \tc{Time control}\item[(T)otal] This is the amount of time to integrate. It is called {\tt TEND} in the documentation. If it is negative then it will be made positive and no data will be stored. Thus you can integrate for very long periods of time without being told that the storage is full. \item[(S)tart time] This is the initial time {\tt T0} For autonomous systems it is usually irrelevant. \item[tRansient] The program will integrate silently for this amount of time before plotting output. It is used to get rid of transients. \item[(D)t] This is the step size used by the fixed step integrators and the output step for Gear's algorithm. It is positive or negative depending on the direction of integration. \tc{Iteration control} \item[n(C)line ctrl] will prompt you for the grid size for comptuing nullclines. \item[s(I)ng pt ctrl] This prompts you for errors and epsilons for eigenvalues and equilibria calculations as well as the maximum iterates. \item[n(O)ut] This sets the number of integration steps to perform between output to storage. Thus, if you output every 10 steps with a {\tt Dt} of .05, XPP will yield output at times that are $10*.05=.5$ timesteps apart. The advantage of this is that lengthier records of data can be made without losing accuracy of the integrator. \item[(B)ounds] This ses a global bound on the integrator. If any variable exceeds this value in magnitude, a message appears and the integration stops. \tc{Numerical method}\item[(M)ethod] allows you to choose the methos of integration. Only Gear is adaptive and is the best to use for stiff problems. If you choose Gear, you will be asked for error tolerance, minimum step, and maximum step size. The discrete method should be used for difference equations. If the method is GEAR, the output {\tt NOUT} is set to 1. Also Gear generates several possible error messages that you may have to respond to. These suggest how to fix the error. Backward Euler prompts you for a tolerance and the maximum number of iterates for each step. The Volterra method is described below; it prompts for a tolerance, maximum iterates, and a ``memory size.'' You will also be asked if you want the convolution kernels to be re-evaluated after any parameter is changed in either the range integration or through a manual change in parameters. If this flag is 1 then the kernels will be re-computed. The default is to not recompute them. \item[d(E)lay] This sets the upper bound for the maximum delay. If in your delay equations, the delay exceeds this, a message will appear and the integration will stop. Each time this is changed all previous delay data is destroyed and you must begin your integration anew. Thus, it should be the first thing you set when solving a delay equation. Since the storage depends on the size of {\tt Dt} when you change this, then the delay storage will also be destroyed. Delay equations require data for $t.) To activate the DB, bring the window to the top and put the pointer inside the box. Across the top is a menu of commands and then there follows a list of titles for the variables, time and the auxiliary functions. Using the arrow keys, the page keys or clicking on the appropriate commands allows you to scroll through the data. By resizing the DB you can get more or less data. Clicking on Left shifts the data window to the left and Right moves it to the right. The time column always stays fixed. If you do parametric or range calculations, the range variable is kept in the time column. Home takes you to the top of the data and End to the last row. The remaining commands will be described separately. The keyboard shortcuts to invoke them are in parentheses. \begin{description} \item[(F)ind] This pops up a window and asks for a variable and a value. It then looks through the data until it comes to the closest value to the specified that the variable takes. It only moves down the file so that you can find successive values by starting at the top. \item[(G)et] This makes the top row of data the initial coonditions for a new run. \item[Re(p)lace] This pops up a window asking you for a column to replace. Then it prompts you for a formula. Suppose you want to replace {\tt AUX1} with {\tt x+y-t} where {\tt x,y} are two variables. Then type this in when prompted and the column that held {\tt AUX1} will be replaced by the values in these columns. Any valid XPP function or user function can be used. Two special symbols can also be used when applied to single variables: \tc{Numeric integration/differentiation}\begin{description} \item[@VARIABLE] replaces the column with the numerical derivative of the variable. You cannot use this within a formula, but once the column is replaced, it can be treated as any other column. \item[\&VARIABLE] replaces the colmun with the numerical integral. Note that successive applications of the derivative and then the integral will result in the original plus a constant. \end{description} \item[(U)nreplace] This undoes the most recent replacement. \item[Fir(s)t] This marks the top row in the DB (nothing is shown)as the start or first row for saving or restoring. \item[Last (e)] This marks the top row as the last or end row for saving or restoring. \item[(R)estore] replots the data marked by First and Last. The default is the entire data set \tc{Writing output}\item[(W)rite] prompts you for a filename and writes the marked data to a file. The files are of the form: \begin{verbatim} t0 x1(t0) ... xn(t0) t1 x1(t0) ... xn(t0) . . . tf x1(tf) ... xn(tf) \end{verbatim} and reflect the current contents of the DB. They are ASCII readable. \item[(L)oad] Will load in as much of a similarly formated data file as possible for graphing. \end{description} \section{C Files} For some types of problems that have lengthy and complicated right-hand sides, you would like to use directly compiled C code rather than using the built in interpreter. I have found that one generally does not get as much speed up as merits the extra work of recompilation for each problem. Nevertheless, I have included an option for doing that in XPP. If you want to do this for a particular ODE file, simply append the call to XPP with the option {\tt -m} which means to make a C file. XPP goes on as before but when you exit the program, there will be two files called {\tt my\_cfile.c} and {\tt my\_hfile.h.} These files should help you in your creation of your own file. Basically, you should replace the call to {\tt my\_rhs()} in {\tt my\_rhs.c} with the contents of {\tt my\_cfile.c.} At this time, this is incompatible with functional equations, so dont use it if you are solving them. Nor is it compatible with tabulated stuff. For example, the Morris-Lecar model \begin{verbatim} 2 variables v,w fixed m params v1=-.01,v2=0.15,v3=0.1,v4=0.145,gna=1.33,phi=.333 params vk=-.7,vl=-.5,iapp=.075,gk=2.0,gl=.5,om=1 user minf 1 .5*(1+tanh((arg1-v1)/v2)) user ninf 1 .5*(1+tanh((arg1-v3)/v4)) user lamn 1 phi*cosh((arg1-v3)/(2*v4)) odev (iapp+gl*(vl-v)+gk*w*(vk-v)+gna*m*(1-v))*om odew (lamn(v)*(ninf(v)-w))*om odem minf(v) b v-v' b w-w' done \end{verbatim} produces the header file: \begin{verbatim} #define v y__y[0] #define w y__y[1] double m; #define v1 constants[1] #define v2 constants[2] #define v3 constants[3] #define v4 constants[4] #define gna constants[5] #define phi constants[6] #define vk constants[7] #define vl constants[8] #define iapp constants[9] #define gk constants[10] #define gl constants[11] #define om constants[12] double minf(); double ninf(); double lamn(); \end{verbatim} and the C file \begin{verbatim} #include "my_hfile.h" extern double constants[]; double minf(arg1) double arg1; { return(.5*(1+tanh((arg1-v1)/v2)) ); } double ninf(arg1) double arg1; { return(.5*(1+tanh((arg1-v3)/v4)) ); } double lamn(arg1) double arg1; { return(phi*cosh((arg1-v3)/(2*v4)) ); } my_rhs(t,y__y,ydot,neq) double t,*y__y,*ydot; int neq; { do_fix(); ydot[0]= (iapp+gl*(vl-v)+gk*w*(vk-v)+gna*m*(1-v))*om; ydot[1]=(lamn(v)*(ninf(v)-w))*om; } do_fix() { m=minf(v); } \end{verbatim} Copy the file, {\tt my\_rhs.c} to some backup. Delete the function called {\tt my\_rhs(...)} and substitute the contents of {\tt my\_cfile.c} putting the top declarations of {\tt my\_cfile.c} at the beginning of {\tt my\_rhs.c}. Then, remake XPP. Invoke it as before by typing {\tt xpp} followed by the file name. \subsection{Warnings} This tool is only a guide to the production of C code. There are several things to watch out for: \begin{enumerate} \item C is case sensitive and XPP is not so make sure you are consistent. \item If your parameters are variables are the same as any of the local variables or key words in the C file {\tt my\_rhs.c} then there will be compilation errors. For example, don't use {\tt i.} When in doubt, capitalize all your parameters and variables in the ODE file keeping the system functions such as {\tt exp} in lower case. \item Some XPP functions like {\tt heav} and {\tt max} are not known to C so you must define them. I will probably fix this later. \item The expression {\tt y**x} or {\tt y\^x} makes no sense in C and should be replaced by {\tt pow(y,x)}. \end{enumerate} \section{The options file} You can have many options files. They are useful for initializing XPP. They have the following format: \bvb 9x15 BIG_FONT_NAME <-- menu fonts fixed SMALL_FONT_NAME <-- IC, browser, etc fonts 0 BACKGROUND (1=white,0=black) 0 IXPLT <-- X-axis variable (0=time) 1 IYPLT <-- Y-axis variable 1 IZPLT <-- Z-axis variable 3 AXES <-- type of axis 1 NJMP <-- nOutput 40 NMESH <-- Nullcline mesh 4 METHOD <-- Integration method 1 TIMEPLOT <--- set to zero if one axis is not time 8000 MAXSTOR <--- maximum rows stored 20.0 TEND <-- total integration time .05 DT <--- time step 0.0 T0 <-- start time 0.0 TRANS <-- transient 100. BOUND <--- bounds .0001 HMIN <-- min step for GEAR 1.0 HMAX <-- max `` `` `` .00001 TOLER <-- tolerance for GEAR 0.0 DELAY <-- maximal delay 0.0 XLO <--- 2D window sizes 20.0 XHI -2.0 YLO 2.0 YHI \end{verbatim} Within an ODE file, you can call up a different options file by typing \bvb c \end{verbatim} where {\tt } is the name of an options file. \section{Functional equations} In the course of some research problems, I was pursuing, I ran into some Volterra equations of the form: \[ u(t)=f(t)+\int_0^t K(t,s,u(s))ds \] that I could not convert to ODES. (If $K$ is a convolution with a sum of powers and exponentials, then it can be converted to an ODE. Since XPP is much more efficient with ODEs and has been thoroughly debugged with respect to them, you should always attempt this conversion first.) Thus, I have added the capability to solve equations with this type of term in them. This has necessitated the addition of two new commands for the ODE file and a solver for such problems. The solver is described below in the numerical section. Since the equation above requires ``memory'' all the way back to $t=0$ and one often is interested in long time behavior, XPP truncates this to: \[ u(t)=f(t)+\int_{\max(t-T,0)}^t K(t,s,u(s))ds \] where {\tt T=DT*MaxPoints}, the latter being a parameter that you set in the numerics menu. The default is 4000. Thus, one assumes that the kernel function decays for large $t$ and so the tail will be small. As {\tt MaxPoints} is a parameter, you can always make it larger at the price of taking longer to evaluate right-hand sides. Let us consider the following equation: \[ u(t) = \sin(t)+\frac{1}{2}(\cos(t)-\exp(t)-t\exp(t))+\int_0^t (t-s)\exp(-(t-s))u(s)ds \] whose solution is $u(t)=\sin(t).$ The following ODE file will create this model and also add an auxiliary variable with the solution for purposes of comparison. \begin{verbatim} 2 v u k w 0.0 (t-t')*exp(t'-t)*u i sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+w o sin(t) d \end{verbatim} There are two new ``commands'': The ``integral'' operator {\tt k} declares that the symbol {\tt w} will denote an integral with respect to $t'$ of the expression. The number directly following the name ``w'' denotes the amount of a weak singularity; thus: \begin{verbatim} k w K(t,t',u) \end{verbatim} will attempt to evaluate: \[ \int_0^t (t-t')^{-mu} K(t,t',u(t')) dt' \] \tc{Convolution}If $mu=0$ then no special integration methods are used; otherwise, a product integration technique is utilized. If your problem can be cast as a convolution problem, considerable speedup can be obtained since lookup tables are created. For example, the present example is in fact a convolution problem, so that instead of the full declaration, one could instead write: \begin{verbatim} k w 0.0 t*exp(-t)#u \end{verbatim} which convolves the first expression (of $t$ {\bf only}) with $u.$ For this example, using a time step of .05 and integrating to 40, there is a 3-fold speed-up using the convolution. For more complicated kernels, it will be more. The ``k'' declaration is unnecesary and the above file can be written more simply as: \begin{verbatim} 2 v u i sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{t*exp(-t)#u} o sin(t) d \end{verbatim} Here the $mu$ is set to 0 by default. If one wants to solve, say, \[ u(t) = exp(-t) + \int^t_0 (t-t')^{-mu} K(t,t',u(t'))dt' \] the file is: \begin{verbatim} 2 v u i exp(-t) + int[mu]{K(t,t',u} o sin(t) d \end{verbatim} Note that the ``mu'' must be a number. The expression \[ {\tt i } \] tells XPP to interpret the right-hand-side to mean {\tt u(t) = } instead of $u'(t)={\tt }$ which is implied if the first letter is {\tt o} instead of {\tt i.} \tc{Auto evaluation}\begin{center} {\bf Warning} \end{center} If you have parameters in your definition of the kernel and you change them, then you will have to go into the numerics menu and recompute the kernels by calling the Method command which automatically recomputes the kernels. Alternatively, turn on the AutoEval flag and it will be done automatically. \tc{Example}I close this section with an example of a pair of coupled oscillators in an infinite bath which has diffusion and passive decay. The equations are: \beqann u(t) &=&\int_0^t k(t-s)F(u(s),v(s))ds + \int_0^t k_d(t-s)F(u_1(s),v_1(s))ds \\ v(t) &=&\int_0^t k(t-s)G(u(s),v(s))ds + \int_0^t k_d(t-s)G(u_1(s),v_1(s))ds \\ u_1(t) &=&\int_0^t k(t-s)F(u_1(s),v_1(s))ds + \int_0^t k_d(t-s)F(u(s),v(s))ds \\ v_1(t) &=&\int_0^t k(t-s)G(u_1(s),v_1(s))ds + \int_0^t k_d(t-s)G(u(s),v(s))ds \eeqann where \beqann k(t) &=& \exp(-t)/\sqrt(\pi t) \\ k_d(t) &=& \exp(-t)\exp(-d/t)/\sqrt(\pi t) \eeqann and \beqann F(u,v) &=& \lambda u -v - (u+qv)(u^2+v^2) \\ G(u,v) &=& \lambda v + u - (v-qu)(u^2+v^2). \eeqann Note that the kernel $k$ is weakly singular and thus $\mu=.5.$ There is a singularity at $t=0$ for the diffusive kernel; this can be rectified by adding a small amount to the denominator. The XPP file is as follows \begin{verbatim} 4 # the four variables: v u,v,u1,v1 # the functions F(u,v), G(u1,v1) f f,g,f1,g1 p lam=1.5,q=0.8,d=1,u0=1,u10=.95 # 1/sqrt(pi)= p spi=.5641896 # the integral equations; since (0,0,0,0) is a rest point, I # add a small quickly decaying transient i u0*exp(-5*t)+spi*(int[.5]{exp(-t)#f}+int[.5]{exp(-t-d/(t+.0001))#f1}) i spi*(int[.5]{exp(-t)#g}+int[.5]{exp(-t-d/(t+.0001))#g1}) i u10*exp(-5*t)+spi*(int[.5]{exp(-t)#f1}+int[.5]{exp(-t-d/(t+.0001))#f}) i spi*(int[.5]{exp(-t)#g1}+int[.5]{exp(-t-d/(t+.0001))#g}) # the four functions f,g,f1,g1 o lam*u-v-(u*u+v*v)*(u+q*v) o lam*v+u-(u*u+v*v)*(v-q*u) o lam*u1-v1-(u1*u1+v1*v1)*(u1+q*v1) o lam*v1+u1-(u1*u1+v1*v1)*(v1-q*u1) d \end{verbatim} Try it. \section{Good luck} Look at the example ODE files to see how to write your own. Some things I will do shortly: \begin{enumerate} \item Better graphics including HPLG plotter support \item Interface to AUTO the bifurcation program \item PDE version of XPP. (This is now working and handles 1-dimensional integral and partial differential equations) \item Better delay integrator. \item As always make it faster. \item FFTs anyone? \end{enumerate} \tc{Numerical methods}{\bf Some comments on the numerical methods} Most are standard. The BVP solve works by shooting and using Newton's method. All Jacobi matrices are computed numerically. The nullclines are found by dividing the window into a grid and evaluating the vector field at each point. Zero contours are found and plotted. Equilibria are found with Newton's method and eigenvalues are found by the QR algorithm. Invariant sets are found by using initial data along an eigenvector for the corresponding eigenvalue. The eigenvector is computed by inverse iteration. The solvers are standard. The Gear algorithm is out of Gear's text on numerical methods. Delay equations are solved by storing previous data and linearly interpolating from this data. Adjoints are computed for stable periodic orbits by integrating the negative transpose of the numerically computed variational equation backwards in time using backward Euler. The Volterra solver uses essentially an integrator (second order) based on the implicit product scheme described in Peter Linz's book on Volterra equations (SIAM,1985). For ODEs implicit schemes take considerably more time than explicit ones, but since most of the compute time for Volterra equations is in approximating the integral, this time penalty is minimal. Performace is gained primarily by taking advantage of convolution type equations. Normally distributed noise is computed by the Box-Muller transformation of uniform noise. \end{document}