ATH 437/537 非线性系统

对于经验丰富的用户来说，XPP/AUTO可以快速有效地分析常微分方程的非线性系统。然而，学习曲线可以起初相当陡峭，并且用户界面有时可能崩溃和/或行为不可预测，
所有这些都可以容易地挫败新用户。本文档和教程的目标是通过使用XPP/AUTO可以执行的一些基本功能希望能对其用法有所澄清。我将概述您开始使用XPP和使用它所需的一些关键点当然，尽管你可以用它做更多的事情，我不会涵盖。如果一开始是压倒性的，不要太担心；经过足够的练习键盘快捷键、最佳参数等都会刻在你的脑海中，最终
它们变得自然而然，甚至不需要再去想它们。请确保首先在下载（请按照特定于您的操作系统的下载说明进行操作）。如果您对本计划有任何其他问题，但此处未涉及，或者您想要更深入的解释，或者如果你想看更多的例子，请参阅在线文档MATH 437/537 – XPP/AUTO Tutorial

1 Introduction

To the experienced user, XPP/AUTO can be a quick and efficient way to analyze the basics ofa nonlinear system of ordinary differential equations (ODEs). However, the learning curve canbe quite steep at first and the user interface may sometimes crash and/or behave unpredictably,all of which can easily frustrate a new user. The goal of this document and tutorial is to walkthrough some of the essential functions that can be carried out using XPP/AUTO, and shouldhopefully offer some clarity into its usage.I will outline some of the key points that you need to get you started with XPP and to use itfor the course, although there are many more things that you can do with it that I won’t cover.Don’t worry too much if it’s understandably overwhelming at first; with enough practice, thekeyboard shortcuts, optimal parameters, etc. become engraved in your brain and eventuallythey become automatic without having to even think about them anymore.Make sure to first download XPP/AUTO at http://www.math.pitt.edu/~bard/xpp/xpp.html(follow the download instructions specific to your operating system).If you have any further questions about the program that are not covered here, or if youwould like more in-depth explanations, or if you want to see more examples, please referto the online documentation at http://www.math.pitt.edu/～bard/xpp/help/xpphelp.html2 Creating an ODE fileIn order to analyze the system of ODEs in question, you need to write a special file called anODE file (file extension: .ode). You can write one using a standard text file editor, making sureto save it as using a .ode file extension (I recommend using Sublime Text but use whichever texteditor you like).Say we want to analyze the following system of ODEs describing a simple SIR model(standardepidemiological model of susceptible, infected, and recovered individuals), with the R variable2
The system described by Eqs. (1)-(3) contains two differential equations (S and I with respect totime), one so-called “auxilliary” variable (R with respect to I), and 4 constant parameters (pi, a,b, and μ). A typical ODE file written to analyze this system in XPP/AUTO could look like this(see SIR.ode)

2.XPP/AUTO Tutorial

? Equations: this is where you put your differential equations. You can write these in dif-ferent ways, either as S’ = ... (as shown above), or more explicitly as dS/dt = ...; bothwork exactly the same.? Functions: these are the auxilliary functions, i.e. the quantities which are functions of oneor more of your time-dependent variables (in this case R is a function of I(t) as in Eq. (3)).Including the aux prefix allows you to still plot R in XPP (but not in AUTO). WARNING:Do not forget to write R(I) (instead of just R); XPP will still compile your file, and it willseem like it works, but sometimes it won’t! (learned the hard way)? Parameters: these are the constant parameters of your model (4 of them in this case). Youcan start the line with either p or par or even parameter, and can put multiple parameterson one line separated by commas as I did on line 10.? Numerics: These let you modify the internal settings within XPP/AUTO that have to dowith various things such as how you want XPP to display time series/bifurcation plots(which variables to plot, the range, etc.), and also numerical parameters like min/maxstep size (dsmin, dsmax), upper bound, etc. We’ll go over these in more detail later. Youtechnically don’t need to set any of them beforehand since they can all be changed withinthe program, but often it is helpful to do so since you save yourself the trouble of resettingthem manually each time you open a file. IMPORTANT: exceptionally in thenumericssettings, you cannot add spaces before or after the commas and equal symbols.? Initial Conditions: Set the value of your time-dependent variables at the first time point,i.e. at t = 0. If you leave these out, they default to 0. You could also write init S=0.5, e.g..Note that variable and parameter names are NOT case sensitive, so for example you cannotname one variable as R and another as r because XPP will register these as the same variable/-parameter. Consequently, you should not name a variable T, because T (i.e., t) is reserved forthe independent variable (time). Other keywords (like ’pi’) are also reserved for other purposes.There are a lot more things you can write into .ode files that XPP can handle (e.g. step file will fail to open because of syntax or other errors in your .ode file.If this is the case, make sure you read the output text on the XPP Shortcut/Terminal window,there should be an explanation of the error at the end of this text.3.1 Running XPP/AUTO on WindowsIf you are running a Windows machine, then you first need to run Xming before opening your.ode file (attempting to run the ODE file before doing so will result in an error). Then, in orderto start XPP/AUTO, you can click and drag the .ode file, e.g. from your File Explorer, anddrop it into the XPP Desktop shortcut (instructions for setting up example, to open the ODE file called SIR.ode, first ensure that Xming is running, thendrag the file SIR.ode and drop it into the XPP Shortcut on the Desktop.3.2 Running XPP/AUTO on Mac or LinuxTo run the program if using a Linux machine or a more recent version of Mac OS, you shouldgo through the system terminal. Notably, for Mac users, the ability to drag and drop ODEfiles into the Desktop shortcut is no longer available (because Mac stopped supporting 32-bitprograms for some reason). If you can get it to work, great, but if not, follow the instructionsbelow.First, open a terminal window (for Mac users that may not be familiar with terminal, do a Spot-light search and type in ’Terminal’ and it will come up). Then, you’ll want to navigate to thedirectory containing the .ode file, using the cd command (‘change directory’). For example, onmy Mac I have the two ODE files, SIR.ode and dengue.ode in a subfolder called The first line of code brings me into the directory containing the file of interest; this line of codewill vary from person to person depending where you store the file called SIR.ode. The sec-ond line varies for Mac and Linux users; for Linux omit /Applications/ and just write xppautSIR.ode.Of course, later on when you want to open different .ode files, make sure you change the nameof the file from SIR.ode to whatever you named the file in question.A few notes for those of you who might be completely unfamiliar with Shell commands:? You can always use the ls command in the terminal window; typing “ls” will display allfiles and folders within the current directory, which can be useful if for example you forgetthe file or folder names. You can use the Tab key to auto-complete in some cases.? If you want to list only files of a given type (e.g. only the .ode files in the directory) youcan type ls *.ode.? If you mistakenly went into the wrong directory and want to go back to the previous one,type “cd ..”
Figure 1: The XPP window that should appear when the program successfully reads SIR.ode? Param: View and modify, or reset, model parameters. Tip: if the parameter list is long andcut off, press the full screen button to show all parameters at the same time.? Eqns: View model equations (ODEs). Do not attempt to modify the model equationswithin XPP, this has never worked for me. If you need to change it, close the program,modify the source file (.ode), save, and recompile in XPP.? Data: View the data of each variable at each time point of the integration (note: this will beblank if you haven’t run the simulation yet).On the left hand side, you have all the options for running simulations, plotting nullclines, mod-ifying numerics (e.g. step size, integration method, etc.), modifying view options (e.g. whichvariables to plot, time series vs. phase portrait), adding curves, solving for fixed points, plottingmanifolds, and much more. KEYBOARD SHORTCUTS: Note that each button in this columnhas one uppercase letter, e.g. (I)nitialconds or ph(A)sespace. You can navigate these differentoptions and tools using your cursor, but for the most part and as you get used to the program,you learn to use these uppercase letters as keyboard shortcuts to perform tasks very quickly andefficiently. We’ll now go through some of the most important of these functions.4.1 Time series solutionsNotice at the top of the program above the empty plot, the title reads ”S vs T” meaning variableS with respect to time. This is because under the ’Numerics’ section in the ODE file, on line 15,7XPP/AUTO Tutorial
we set XP1=t (i.e. plot time on the x-axis) and YP1=S (plot S on the y-axis). If these are left out,then XPP will simply plot the first variable in the equations list w.r.t. time.1. To runthesimulation/integrate the solution using the default initialconditions (ICs) andparameters set in the ODE file, press(I)nitialconds → (G)o, or simply type I then G onyour keyboard. You should get the following trace to appear (Fig. 2).Figure 2: By default, XPP is plotting YP1 vs. XP1, which we set to S and t, respectively, in the ODE file.2. To view the trace for I(t) on the same graph, press (G)raphic stuff → (A)dd curve, thenset the Y-axis variable to I, and change the colour by typing in an integer from 1 to 10(0 is black, you can play around with the colours to see what you get for each number).Choosing 1 yields red as in Fig. 3.Figure 3: Notice the new red trace representing As you see in Fig. 3, S(t) (black) and I(t) (red) converge to a stable steady-state. A veryuseful feature, and one way to see the numerical values of S and I at steady state, is to press(I)nitialconds→ (L)ast (do this a few times). This makes the ICs of the new simulation thelast data point of the previous simulation, in this case near the steady state, so each timeyou perform (I)→ (L) you get closer and closer to the steady state (Fig.4).Figure 4: Performing the (I) → (L) operation 3 times, solutions converge closer and closer to the truesteady state value at S = 2 and I = 0.5 for the default parameter values. Window at top-left is the ’ICs’window. The sliders that you see at the bottom of the XPP window allow you to vary up to 3 pa-rameters at a time on a sliding scale, between a lower and upper range that you set. Forexample, try setting a slider for the parameter a between 0 and 1 (Fig. 5), and dragging thesmall black square across that range. You should see S(t) and I(t) vary as you vary a.Figure 5: Slider window to quickly scan over parameter value To delete any/all the extra curves and prevent them from being plotted in the future (e.g.the red trace in previous Figs), press/type (G)raphic stuff → (D)elete last, or (G)raphicstuff→ (R)emove all (to clear everything but the original trace).? Running new simulations simply plots new traces on top of the old ones. To tidy up,simply press/type (E)rase and they should be cleared off the plot.4.2 Phase portraitUp until now, we’ve only been plotting solutions with respect to time. Often, however(andespecially with 2D systems like the one we effectively have), it is useful to plot the variables withrespect to each other. Before proceeding, make sure you set parameters back to their defaults(click ’Param’ at the top, and ’Default’ in the pop-up window).1. To change the axes, press (V)iewaxes→ (2)D. A window will appear allowing you to setthe axes, ranges, and labels; set them as shown below in Fig. 6.Figure 6: Changing what’s being plotted on the x- and y-axes.2. In this view, we can do a few cool things. For example, we can get a good idea of the vectorfield by pressing (I)nitialconds→ m(I)ce. This allows you to click anywhere on the graphto plot many traces from different initial conditions, and eventually a vector field revealsitself as in Fig. 7 and we can see all solutions converging to a single point in the graph.When you’re satisfied, press the ‘Escape’ key on your
Figure 7: Solution flows in phase space reveals behaviour of system (evolution to stable focus).3. Let’s plot nullclines! First, we need to tweak some numerics to get smoother, more accurateestimates from XPP as to where the nullclines actually are. Go into the n(U)merics menu,and you will see a bunch of options on the left. In this menu, you can change things likeintegration method (M), step size ?t in the integration (D), total integration time (T), upperallowable bound for variables to assume (B), and more.4. Right now, we are interested in nullcline numerics; press (N)cline ctrl and increase “nclinemesh’ at the top from 40 to at least 400 (the larger this value, the more mesh points usedin computing the nullclines and the smoother and more accurate they will be). Then press’Escape’ to go back to the main menu5. To plot the nullclines, press (N)ullclines → (N)ew. They should appear as in Fig. 8, withx-axis nullclines shown in orange and y-axis nullclines in green.Figure 8: One S (horizontal)-nullcline plotted in orange, two I (vertical)-nullclines in green.11Determining equilibria and plotting manifoldsThere are a few ways to find steady states. In anticipation of performing continuation analysiswith AUTO (i.e. plotting bifurcation diagrams), it is necessary that the data stored in the ICscolumn contains the values of the variables at, or near enough to, a fixed point; failure to doso will result in convergence error once you get to the AUTO stage. Below are a few ways toaccomplish this We’ve already seen the first way to reach a stable fixed point, which is to run (I)nitialconds→ (L)ast a few times until the ICs converge to the fixed point (refer to Fig. 4)2. A second way is to use the built-in tool that computes fixed points. To do so, click/type(S)ing pts→ (G)o. It will search for a numerical fixed point in the neighbourhood of thevalues stored in the ICs data. Then it will ask if you want to print eigenvalues; if you click‘Yes’ they will appear in the Terminal window. Important: make sure to click “Import”in the ‘Equilibria’ window (I like to click it twice); this write the equilibrium point to theICs data, effectively replacing the need to perform (I)→ (L) as many times as youwouldotherwise.Figure 9: Using (S) → (G) to find the stable equilibrium point. Eigenvalues are shown in the terminalwindow; both are negative and complex indicating that the fixed point at S = 2, I = 0.5 is a stable focus.Clicking “Import” twice in the ‘Equilibria’ window (top left) sets the ICs to these fixed-point values.3. Another way to converge to a fixed point involves using the phase portrait configurationwith the nullclines plotted, as in Fig. 8. The two intersections of green and orange nullclinesrepresent the two fixed points of the system. By clicking on (S)ing pts→ (M)ouse, you use12
XPP/AUTO Tutorialyour mouse to click on/near the left-most intersection (the stable fixed point). The same‘Equilibria’ window will pop up as before, as well as the option to display eigenvalues(Fig. 10). Again, make sure to press ‘Import’ before continuing toward AUTO.Figure 10: Using (S)→ (M) to click on and converge to the stable equilibrium point. A circle appears overthe stable fixed point when in phase portrait mode (hard to see here but it’s there).Plotting manifolds of a saddle point can be done by following the same instructions as in step3 above, only this time click on the saddle point, and make sure to press ‘Yes’ (or type Y) whenprompted to “Draw Invariant Sets?”; the result is shown in Fig. 11.Figure 11: Using (S) → (M) to click on and converge to the saddle point (eigenvalues: ?1 and 1). XPPmarks the saddle point with a triangle. Stable manifolds shown in cyan, unstable manifolds in ? have all parameters set to their default values (Param window: click “Default,” then click“Close” to close the window (note: pressing your system’s built-in close button won’t workin XPP/AUTO, you have to use the buttons marked “Close” within the program).? have imported the stable fixed-point values of S = 2 and I = 0.5 (when a = 1); followinstructions in the previous section.Side note: AUTO is actually its own program, but the creator of XPP/AUTO, Bard Ermentrout,embedded AUTO into XPP in order to use time series simulations and bifurcation analyses inconjunction with each other.5.1 Single-parameter steady-state bifurcationsLet’s create a bifurcation diagram with a ∈ [0, 1] as our bifurcation parameter/range. To start:1. Launch the AUTO window: press (F)ile → (A)uto; the AUTO window should pop up(Fig 12). Set which parameter you want as your bifurcation parameter in the (P)arameter window.Make sure to set Par1 to the parameter you want to vary, in our case to a (see Fig. 13).Note that, if you intend to do a two-parameter bifurcation, this is the moment for youto set Par2 to be your second parameter of interest; if you run the continuation beforechecking for this and realize after, you’ll have to restart and no one wants that (more ontwo-parameter bifurcations in the last subsection).Figure 13: Setting Par1 as your primary bifurcation parameter, and Par2 as your second parameter (fortwo-parameter bifurcations).3. Open the (N)umerics window, and set the values such that they match those shown inFig. 14, click ‘Ok’ or press Tab key when you’ve finished.Figure 14: Numerics window. Do not change the values for those parameters boxed in red.There are many AUTO numerics to set, and each affects the output of your bifurcation diagramin different ways. Below (Table 1) is a list of what these actually are, as welldocumentation).4. Next, set the diagram axes for the bifurcation diagram. Press/type (A)xes → h(I)-lo andset the values in the ‘AutoPlot’ window to those shown in Fig. 15, then press ‘Ok’/Tab key.Figure 15: Setting the axes and their limits.5. We’re ready! To run the bifurcation, press (R)un → (S)teady state. You should get thebifurcation diagram shown in Fig. 16.Figure 16: AUTO-generated bifurcation diagram, which reveals a transcritical bifurcation when a ≈ 0.67.Stable stead-state branches are plotted in red, unstable steady state branches in black.Pay attention to the output text in the Terminal window (Fig. 17). AUTO produces a label forpoints along the computation every Npr points, whenever a special point (e.g. bifurcation point)is found, and when the endpoints are attained (e.g.
Param. Description Value?Ntst Number of mesh intervals for discretization of periodic or-bits; increase if periodic solution
Nmax The maximum number of points that AUTO will computebefore terminating the computation; increase if you find thatthe computation terminates prematurelyDepends, ≥ 500NPr Every Npr steps, AUTO spits out information about whereit’s at; it doesn’t really matter what value you choose, but toosmall a number will make the diagram, and Terminal output,appear cluttered with tags.50-500 or moreDs The initial step size and direction in the bifurcation parame-ter that AUTO will take; this must be made negative if in-tend to move from right to left; try decreasing its absolutevalue if you get an MX error in first stepDepends, (+) orDsmin Smallest allowable step that AUTO takes (step size is adap-tive); try decreasing if AUTO throws an MX error, or if it ap-pears to be missing bifurcation pointsDepends
EPSL Tolerance criterion (convergence of equation parameters); noneed to worry about what it is but keep it small10?7 to 10?6Dsmax Largest allowable step that AUTO takes; increase if going tooslow, but be careful (too large and it can miss bfn points)DependsPar Min Minimum allowable value for bfn parameter set as Par1 DependsPar Max Maximum allowable value for bfn parameter set as Par1 DependsNorm Min Minimum allowable value for the solution’s L2-norm 0Norm Max Maximumvariables in the ODE are large, though itmight take longer if too big (maybe 10 to 100 times the ex-pected L2-norm, if known)DependsEPSU Tolerance criterion (convergence of solution components),keep it small and same as EPSL10?7 to 10?EPSS Tolerance criterion (arclength convergence); keep it small,and about 100-1000 times the value set for EPSU and
Table 1: Parameters within AUTO’s numerics window. IMPORTANT: Values/ranges listed in the right-most column are only my suggestions, please take with a grain of salt and use your judgement
WX：codehelp