u' = -u+f(a_{ee}u-a_{ie}v-z_{e} + I_{e}(t) )
v' = (-v+f(a_{ei}u-a_{ii}v-z_{i} + I_{i}(t) ))/tau
and here is the corresponding ODE file:
# wilson-cowan f(u)=1/(1+exp(-u)) par aee=10,aie=8,aei=12,aii=3,ze=.2,zi=4 par tau=1 i_e(t)=ie0+ie1*sin(w*t) i_i(t)=ii0+ii1*sin(w*t) par ie0=0,ie1=0,w=.25 par ii0=0,ii1=0 init u=.1,v=.05 u'=-u+f(aee*u-aie*v-ze+i_e(t)) v'=(-v+f(aei*u-aii*v-zi+i_i(t)))/tau @ total=40 @ xp=u,yp=v,xlo=-.1,xhi=1,ylo=-.1,yhi=1 doneI have added definitions for the inputs and the function f as well as default initial data and parameters. I also set it up so that the initial window is the phaseplane.
x''' + a x'' + b x' + c x = x^{2}
should be rewritten as
x' = y
y' = z
z' = x^{2} - a z -b y -c x
with the corresponding ODE file
# third.ode x'=y y'=z z'=x^2-a*z-b*y-c*x par a=1,b=2,c=3.7 init x=1,y=0,z=0 aux en=x^2+y^2+z^2 @ total=200 @ xp=x,yp=y,xlo=-2,xhi=5,ylo=-4,yhi=5 doneI have set the plotting parameters and integration time. I have also added another quantity en which is stored and can be plotted as well.
x_{n+1} = a x_{n-1}(1-x_{n})
This is a second order map and XPP can only solve first order equations. We write it as a pair of first order equations:
x_{n+1} = y_{n}
y_{n+1} = a y_{n}(1-x_{n})
and get the corresponding ODE file:
# delayed logistic map x'=y y'=a*y*(1-x) par a=2.2 init x=.2,y=.2 # set method to discrete @ meth=discrete @ total=100 @ ylo=-.1,yhi=1,xhi=100 doneHere the most important thing is that I tell XPP that this is a discrete dynamical system.
x'(t) = -g x(t) + b f(x(t-d))
f(x) = x/(1+x^{n})
Here is the ODE file:
# mackey-glass f(x)=x/(1+x^n) x'=-g*x+b*f(delay(x,d)) init x=.2 par g=.1,b=.2,n=10,d=6 @ delay=20 @ total=200 @ xhi=200,yhi=2 doneI have told XPP that the maximal delay to expect is 20.
Delay equations with delay d require initial data on an interval -d < t < 0. You can set these in the ODE file by adding the line:
x(0)=exp(t)*0.2which sets x(t)=0.2e^{t} for -d < t < 0.
NOTE this does not set the initial condition for x at t=0. You can do this by adding one more line after the delayed condition:
x(0)=exp(t)*0.2 init x=.2
x' = F(x,y)
0 = G(x,y)
can be solved in XPP. The DAE
x' = z
z' = y
0 = x+z+y+y^3
has an ODE file:
# dae.ode x'=z z'=y 0=x+z+y+y^3 solve y=-1 init x=2 aux ys=y @ yhi=2.5 doneYou have to tell it which variables to solve for. In this case, we solve for y as a function of x,z . You can give XPP a guess about what the value of y is at the initial conditions for x,z. I have added an additional quantity ys which is the solution to the algebraic condition.
u(t)=sin(t) + 0.5 cos(t)-0.5 e^{-t}(t + 1) + int[u(t') (t-t') e^{t'-t},t'=0..t]
and here is the ODE file
# volterra example 1 u(t)=sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{(t-t')*exp(t'-t)*u} aux utrue=sin(t) doneI have included the true solution as well. This is actually a convolution equation so that you and improve the speed by exploiting that:
# volt2.ode - uses convolution to speed up soln u(t)=sin(t)+.5*cos(t)-.5*t*exp(-t)-.5*exp(-t)+int{t*exp(-t)#u} aux utrue=sin(t) done
Integro-differential equations are also possible:
u'(t)=-2u+e^{-2t} + int[u^2(t')e^{t'-t},t'=0..t]
with u(0)=1 has a solution u(t)=e^{-t} as can be confirmed with the ODE file
# nonlinear integrodifferential volterra equation u'=-2*u+exp(-2*t)+int{exp(-t)#u^2} init u=1 aux utrue=exp(-t) doneFinally, integral equations with singularities are easily solved. For example
u(t)=exp(-t)-int[e^{-(t-s)}(t-s)^{-¼} u(s),s=0..t]
This has a removeable singularity in that the power of t is less than 1. The ODE file has the form:
# singular integral equation u(t)=exp(-t)-int[.25]{exp(-t)#sin(u)} done
WARNINGS
dx=z f(x) dt + s dW
The variable z randomly switches back and forth between 0 and 1 at rates r_{0},r_{1} respectively. dW is a Wiener process. The function f(x) is an asymmetric periodic potential with zero mean. Here is the ODE file
# a simple thermal ratchet wiener w par a=.8,s=.05,r0=.1,r1=.1 f(x)=if(x<a)then(-1)else(a/(1-a)) x'=z*f(mod(x,1))+s*w # z is two states markov z 2 {0} {r0} {r1} {0} @ meth=euler,dt=.1,total=2000,nout=10 @ xhi=2000,yhi=8,ylo=-8 doneThe rates in a Markov process need not be constant but can depend on any of the variables or parameters. Always use Euler's method for stochastic equations.
u' = k_{4} (v-u)(a + u^{2}) - k_{6}
v' = k_{1} -k_{6} u
m' = b m
with the condition that when u falls below 0.2 the mass is halved. We tell XPP to look for crossings of u and to cut the mass in half. Here is the ODE file
# tyson.ode init u=.0075,v=.48,m=1 p k1=.015,k4=200,k6=2,a=.0001,b=.005 u'= k4*(v-u)*(a+u^2) - k6*u v'= k1*m - k6*u m'= b*m global 1 {0.2-u} {m=.5*m} @ total=1000,nout=10 @ yp=m,xhi=1000,ylo=0,yhi=2 done
u'' = sin(t u)
u(0)= 1
u(2) = 0
We rewrite this as a first order system:
u' = v
v' = sin(t u)
u(0)= 1
u(2) = 0
Here is the ODE file for this:
# nonlinear boundary value problem u'=v v'=sin(t*u) bdry u-1 bdry u' init u=1,v=-1 @ total=2,bell=0,xhi=2 doneBoundary conditions take the form:
bdry F(u,u')where u' represents the value of the variable at the end of the interval. It is not the derivative in this context!! XPP tries to make the functions defined by each boundary condition vanish. The length of the interval is the total amount of time integrated. Note that in this example, I have turned off the stinking bell.
You solve boundary value problems in XPP by using the Bdryval menu item instead of the Initialconds item.
XPP can also solve systems with homoclinic orbits. See the manual for an example.
u_{t} = -u+vu^{2}+d_{u} u_{xx}
v_{t} = a - vu^{2}+d_{v} v_{xx}
which we can discretize yielding:
u'_{j} = -u_{j}+v_{j}u^{2}_{j}+(d_{u}/h^{2}) (u_{j+1}-2 u_{j}+u_{j-1})
v'_{j} = a - v_{j}u^{2}_{j}+(d_{v}/h^{2}) (v_{j+1}-2 v_{j}+v_{j-1})
with appropriate end conditions. We will make this into an ODE file:
# schnakenberg PDE 100 points f(u,v)=-u+v*u^2 g(u,v)=a-v*u^2 u0'=f(u0,v0)+duh*(u1-u0) u[1..99]'=f(u[j],v[j])+duh*(u[j+1]-2*u[j]+u[j-1]) u100'=f(u100,v100)+duh*(u99-u100) v0'=g(u0,v0)+dvh*(v1-v0) v[1..99]'=g(u[j],v[j])+dvh*(v[j+1]-2*v[j]+v[j-1]) v100'=g(u100,v100)+dvh*(v99-v100) par a=1.05,du=.2,dv=3,h=.2 !duh=du/(h*h) !dvh=dv/(h*h) init u0=1.5,u1=1.3,u2=1.1,v0=1.05,v1=1.05,v2=1.05 init u[3..100]=1.05,v[j]=1.05 @ dt=.1,nout=50,total=50,meth=cvode,tol=1e-5,atol=1e-5 doneNote that I have defined the kinetics as a pair of functions. The statements starting with the ! are derived parameters -- parameters which depend on other parameters. I use a stiff integrator CVODE . Plot this using array plots or look at the spatial profile with the transpose operation.
Take advantage of the banded quality of the system and speed up the stiff integrators my an order of magnitude. You have to rearrange the equations a bit. Here is the ODE file:
# schnakenberg PDE 100 points # interlaced f(u,v)=-u+v*u^2 g(u,v)=a-v*u^2 u0'=f(u0,v0)+duh*(u1-u0) v0'=g(u0,v0)+dvh*(v1-v0) %[1..99] u[j]'=f(u[j],v[j])+duh*(u[j+1]-2*u[j]+u[j-1]) v[j]'=g(u[j],v[j])+dvh*(v[j+1]-2*v[j]+v[j-1]) % u100'=f(u100,v100)+duh*(u99-u100) v100'=g(u100,v100)+dvh*(v99-v100) par a=1.05,du=5,dv=75,h=.2 !duh=du/(h*h) !dvh=dv/(h*h) init u0=1.5,u1=1.3,u2=1.1,v0=1.05,v1=1.05,v2=1.05 init u[3..100]=1.05,v[j]=1.05 @ dt=.1,nout=50,total=50,meth=cvode,tol=1e-5,atol=1e-5 @ bandlo=2,bandup=2 doneI tell XPP the system is banded and interlace the two variables. This runs significantly faster than the previous ODE!
u'_{j} = -u_{j} + F(a_{ee}K_{e,j}
-a_{ie} K_{i,j})
v'_{j} = (-v_{j} + F(a_{ei}K_{e,j}
-a_{ii} K_{i,j}))/tau
where
K_{e,j} = SUM(W_{e}(k)u_{j-k},k=-m..m)
K_{i,j} = SUM(W_{i}(k)u_{j-k},k=-m..m)
with some sort of "boundary" conditions at the edges. Here is an XPP file for a network of 201 cells and an exponential weight function:
# Wilson-Cowan network table we % 51 -25 25 .5*exp(-abs(t)/se)/se table wi % 51 -25 25 .5*exp(-abs(t)/si)/si special ke=conv(even,201,25,we,u0) special ki=conv(even,201,25,wi,v0) par se=4,si=2.5 f(u)=1/(1+exp(-u)) par aee=16,aie=10,aei=25,aii=3,ze=4,zi=10 par tau=4 init u[0..4]=1 u[0..200]'=-u[j]+f(aee*ke([j])-aie*ki([j])-ze) v[0..200]'=(-v[j]+f(aei*ke([j])-aii*ki([j])-zi))/tau @ total=60,meth=qualrk,dt=.25 @ yp=u100,xhi=60,ylo=0 @ yp2=u150,nplot=2 doneThe weights are defined as tables which give the number of points in the table, the x-range, and the functional form of the table. You can also read a file to fill a table. Consult the documentation for more. The special declaration convolves the values in the tables with the variables. The even means that the ends are reflected. Note that you must write ke([j]), that is the [ ] must be included. I have also told XPP to plot two things, u100 and u150.