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.


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.

To get figure 2.21, do the following:

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

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:

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!