Hi there,

I'm struggling to understand the dynamics of a 3-D system of differential equations. I am interested in a numerical simulation. A phase diagram would be nice too. I've tried the standard bag of tricks, but remain quite lost.

At the very least, I'm interested in behavior in a neighborhood of a stationary state, if one exists. I have therefore tried to numerically solve for a stationary state of the system (x, y, z). All my variables x, y, z are strictly positive.  I'm in trouble finding the stationary state. In a nutshell, what happens is (see the equations below) that a non-trivial stationary state in the -x- and -y- equations results in the indeterminate form 0/0 for the -z- equation. I look at a graph of the 3 planes associated with the stationary state, and it looks like there may exist a stationary state, but the numerical procedures I have attempted have so far failed to locate it. I cannot be sure, but there is part of me that believes that a stationary state does exist. So my first problem is in numerically evaluating the stationary state of the 3-D system.

The difficulty comes from the fact that in the -z- equation the term dy/dx = (dy/dt) / (dx/dt)  appears. I tried to apply L'Hopital Rule to calculate the stationary state. I treated y=Y(x) and z=Z(x) as functions of x and searched for the limit of m(x) = dy/dx, as both dy/dt and dx/dt go to zero, that is I linearized the numerator and denominator around the stationary state value of x. For some reason, this procedure did not work. I was able to find an expression for m(x) but by inspection of the graph it did not seem to be close to the stationary state. It looked like the linearization was not precise enough to compute an accurate limit of dy/dx at the stationary state.

The dynamic system is a boundary-value problem (BVP). x(0)=x0 is given, while y(T) and z(T) are presumed to land very close to the stationary state after a long (infinite) period of time T, or y(T)=y* and z(T)=z*, where T=infinity.  To simulate the dynamics of the system around the presumed stationary state (which is currently eluding me), I thought of reducing the 3-D system { x(t), y(t), z(t) } into { Y(x), Z(x) } and to treat  { Y(x), Z(x) } as an initial value problem (IVP) with initial values given by the stationary state values, i.e. Y(x*)=y* and Z(x*)=z*. So to follow through with this method (and likewise with a shooting method), I need to be sure a stationary state exists and, if so, I need a precise estimate of the stationary state.

I have plotted the 3 null-clines associated with the stationary state. By inspection it looks like a stationary state exists. When I zoom towards it, it is also clear that the null-clines associated with the -x- and -z- equations are nearly parallel at the stationary state. And this is likely to be the source of my numerical troubles.

QUESTIONS: Does the stationary state exist? Can brute-force yield an accurate estimate of the stationary state (if it exists)?

The worksheet below sets up the model, tries to calculate a stationary state (the value obtained looks wrong), plots the null-clines associated with the stationary state, zooms near a point that looks like a stationary state, and lastly sets up the transformation of the system from BVP to IVP, to be simulated once a stationary state has been found.

Is my system somewhat ill-behaved ?  How can I analyze the dynamics? Any advice will be greatly appreciated!

Thanks for your attention,

Patrick.

The Maple Worksheet is below:

> restart: with(plots): with(DEtools):
> a:= 0.01: b:= 0.05: c:= 0.10:
> f:= x-> 20*x^(0.5):
> g:= x-> x^(-0.5):
> xdot:= diff(x(t),t) = f(x(t))-y(t);
> ydot:= diff(y(t),t) = (g(x(t))-z(t))*y(t);
> zdot:= diff(z(t),t) = (z(t)-a)*(c-z(t))+(z(t)-b)*(g(x(t))-z(t))*y(t)/(f(x(t))-y(t));
> xd:=eval(rhs(xdot),{x(t)=x, y(t)=y, z(t)=z}):
> yd:=eval(rhs(ydot),{x(t)=x, y(t)=y, z(t)=z}):
> zd:=eval(rhs(zdot),{x(t)=x, y(t)=y, z(t)=z}):
Solve for a stationary state  (The steady-state values obtained below don't seem right)
> ss:= fsolve({xd,yd/y,zd*(f(x)-y)}, {x,y,z});
> xss:= rhs(ss[1]); yss:= rhs(ss[2]); zss:= rhs(ss[3]);
> eval(xdot,{x(t)=xss,y(t)=yss,z(t)=zss}); eval(ydot,{x(t)=xss,y(t)=yss,z(t)=zss}); eval(zdot,{x(t)=xss,y(t)=yss,z(t)=zss});
Plot the 3-D nullclines
> plotopts:=style=patchcontour, shading=none,axes=boxed:
> xmin:=0: xmax:=2*xss: ymin:=0: ymax:=2*yss: zmin:=a: zmax:=c:
> p1:=implicitplot3d({xd},x=xmin..xmax,y=ymin..ymax,z=zmin..zmax, numpoints=5000, plotopts, colour=brown): display3d({p1}):
> p2:=implicitplot3d({yd/y},x=xmin..xmax,y=ymin..ymax,z=zmin..zmax, plotopts, colour=green): display3d({p2}):
> p3:=implicitplot3d({zd*(f(x)-y)},x=xmin..xmax,y=ymin..ymax,z=zmin..zmax, numpoints=5000, plotopts, colour=blue): display3d({p3}, orientation=[-90,20]):
Plot a small sphere to locate the computed steady state
> points:=[[xss,yss,zss]]:
> p4:=pointplot3d(points,colour=black,symbol=circle,filled=true,symbolsize=40):
> display3d({p1,p2,p3,p4},axes=boxed,view=[xmin..xmax,ymin..ymax,zmin..zmax],orientation=[60,-70],projection=0.5);
Zoom closer to the steady state
> p1:=implicitplot3d({xd},x=300..330,y=340..370,z=0.055..0.058, numpoints=5000, plotopts, colour=brown): display3d({p1}):
> p2:=implicitplot3d({yd/y},x=300..330,y=340..370,z=0.055..0.058, plotopts, colour=green): display3d({p2}):
> p3:=implicitplot3d({zd*(f(x)-y)},x=300..330,y=340..370,z=0.055..0.058, numpoints=5000, plotopts, colour=blue): display3d({p3}, orientation=[-90,20]):
> display3d({p1,p2,p3,p4},axes=boxed,view=[300..330,340..370,0.055..0.058],orientation=[50,-30]);
Transforming the BVP into an IVP
> ODE1:= diff(Y(x),x) = (g(x)-Z(x))*Y(x) / (f(x)-Y(x)):
> ODE2:= diff(Z(x),x) = (Z(x)-a)*(c-Z(x))+(Z(x)-b)*(g(x)-Z(x))*Y(x)/(f(x)-Y(x)) / (f(x)-Y(x)):
> SYS := {ODE1, ODE2}:
> VAR := {Y(x), Z(x)}:
> INI := {Y(xss)=yss, Z(xss)=zss}:
> SOL := dsolve (SYS union INI, VAR, type=numeric);
> SOL := dsolve (SYS union INI, VAR, type=numeric, method=classical);
> pY:= odeplot(SOL, [x,Y(x)], 0..xss):
> pZ:= odeplot(SOL, [x,Z(x)], 0..xss):
> display({pY});
> display({pZ});
>