Computer assignments from Borelli & Coleman

Using XPP with the book problems

Page 67 - level curves of a function

To plot level curves of a function F(x,y) you plot trajectories for the differential equation

x'= Fy(x,y)
y'=-Fx(x,y)

where Fx(x,y) is the partial derivative of F with respect to x etc. Consider the function x^4+2xy^3-y^2-x^2 . Then the associated ODE is x'= 6xy^2-2y
y'=-(4x^3+2y^3-2*x)

Here is the XPP file bc67.ode

# generic second order equation
f(x,y)=6*x*y^2-2*y
g(x,y)=-(4*x^3+2*y^3-2*x)
#
x'=f(x,y)
y'=g(x,y)
@ xlo=-1.25,xhi=1.25,ylo=-1.25,yhi=1.25
@ xp=x,yp=y
done

Lines beginning with @ are control lines that change the defaults in XPP so you don't have to do it within the program. Thus, I have told XPP the window to use and to put x on the x-xaxis and y on the y-axis; a parametric plot.

Here are several ways in which you can look at the behavior.

Sketch the direction fields with the command Dir.Field/Flow Scaled Dir.Field Accept the default of 10.
Use the Graphics Freeze On-Freeze to permanently store curves and then use the Initialconds mIce command to draw lots of trajectories.
Print you picture by using Graphics Postscript
A quick way to explore the whole phase space is to use the Dir.Field/Flow Flow command with a grid of 5. The trajectories are not saved in this case, so you cannot get a permanent copy. (Unless you capture the window using xv or in Windows with the Alt-PrntScrn )

First order equations in the y-t plane

XPP will not let you draw direction fields and nullclines of scalar time-dependent differential equations. But you can fool it into thinking that it is an autonomous two-dimensional problem:

s'=1
y'= f(s,y)

Note that s is just a time-like variable perhaps shifted. Here is the ODE using the example from Figure 2.2.1 gen1d.ode

# generic 1-d problem with a few parameters
f(t,y)=-t*t/((y+2)*(y-3))
s'=dir
y'=f(s,y)*dir
par dir=1,a=0,b=0,c=0,d=0
@ xlo=-5.1,ylo=-5.1,yhi=5.1,xhi=5.1
@ xp=s,yp=y
@ meth=qualrk,tol=1e-7
done

A number of comments are warranted.

I have included a parameter called dir When this is -1, then time goes backwards and when it is positive 1, time goes forwards. This saves you from having to change the sign of Dt in the Numerics menu.
I have added a bunch of free parameters so that you can study sensitivity to parameters.
The window of the plot is a bit weird to avoid hitting the singularities when direction fields are computed.
I am using a fairly fancy integrator instead of the default.

To get figure 2.21, do the following:

Integrate the equations with Initialconds Go. (Click OK when it says Step size too small -- the problem can't be --> -- solved beyond this point!)
Freeze the curve Graphics Freeze Freeze and accept the defaults.
Change the parameter dir to -1. And integrate again. Cahnge it back to 1 again.
Draw some direction fields Dir.Field/Flow Scaled-Dir.Field and use a grid of 15.

To get Figure 2.2.2 and 2.2.3, use the above file, gen1d.ode but

Edit the function f(t,y) using the File Edit Functions command.
Change the method of integration to Runge-Kutta by clicking Numerics Method Runge-Kutta and then clicking Esc-exit to go back to the main menu.
Draw scaled direvtion fields.
Click Graphics Freeze On-Freeze .
In the initial conditions window, set the initial value of s=-5.
Now we will integrate over a whole range of y-values. Click on Initialconds Range and fill in as follows:
- Range over: y
- Steps: 20
- Start: 3.5
- End: -1.5
and then click OK. The curves will be nicely filled in! (The range command lets you range over a parameter or initial condition; it is very useful for regular series of variations.)

After you are done admiring your work, click on Graphics Freeze Remove to get rid of all these curves. Try to get figure 2.2.5. This should be enough for you to also do problems 1,7 on page 105. You may have to alter the window that you want to graph over.

Sensitivity; page 107

Lets try to get Figure 2.31 and 2.32. Fire up gen1d.ode and do the following:

Edit the function f(t,y) using File Edit Function and type in -y*sin(t)+c*t*cos(t) . Note that this is Kosher since we have declared a few extra parameters.
Set the initial conditions for s=-4 and y=-6 since we want to start at t=-4. (Click on Ok in the initial data window.
Change the graph limits using Window/zoom Window and make xlo=-4, xhi=-1, ylo=-10, yhi=10
Turn on the freeze to keep the curves, Graphics Freeze On-freeze .
Click on Initialconds Range and fill it in as:
- Range over: c
- Steps: 6
- Start: .7
- End: 1.3
and then click OK.

There you go: all the curves are drawn. Now to get Figure 2.3.2, we increase the the size of the window. (XPP has actually computed beyond the window since the defaut integration is for 20 time units.) Use the Window/zoom Window command to increase xhi=4 and yhi=15 . Voila, figure 2.3.2!