A while ago, I tried to understand the dynamics of a tricky 3-D system. Back then I received great help from Robert Israel and Doug Meade. After leaving the problem aside for a while, I'm giving it another shot.

I have a 3-D system of ODEs in the variables (x,c,q) defined on an infinite time range. The initial value of x is given, x(0), but the other values c(0) and q(0) are free. The question is whether it is possible to find values of c(0) and q(0) such that the system will converge towards a critical point.

The first difficulty is in locating critical points, i.e. stationary states. Because the system below contains a singularity at a stationary state, the "critical" points must be understood as a limit that may or may not be approached as time goes to infinity. I have reduced the set of candidate critical points to one element (that for which q=b, where b is a parameter of the system, see below). Denote this critical/stationary point by (xs,cs,qs). The second difficulty is in establishing whether the critical point is reached by a stable (or center) manifold, i.e. whether it is possible to find values of c(0) and q(0) such that, given x(0), it is true that as T goes to infinity x(T)->xs, c(T)->cs, and q(T)=qs. My current assessment, based on the calculations below, is that the critical point may be approached if the initial value of x is below the stationary value but not if it is above it, i.e. x(0)<xs, but not if x(0)>xs. And, if I'm correct, the approach is along a one-dimensional manifold. I'm not 100 percent sure of myself, hence my message on mapleprimes.

The transformed system (u,v,r) appears to exhibit a stable manifold converging to the stationary state u=v=0 and r=b, for r(0) given. Does this indicate that, likewise, the system (x,c,q) has a stable manifold?

many thanks for your help.

My worksheet with comments follows. Unfortunately, today I do not seem to be able to insert figures into the code below.

**> restart: with(plots): plotsetup(default): with(plottools): with(DEtools): with(LinearAlgebra):**

**> Digits:=100: interface(displayprecision=5):**

Parameter Assignment:

**> A:=1: d:=1/10: s:=1:**

**> b1:=2/100: b2:=3/100: w2:=1/2: w1:=1-w2: b:=w1*b1+w2*b2: b; **

The dynamic system of interest:

**> xdot:=diff(x(t),t)=A*(x(t))^(1/2)-d*x(t)-c(t);**

**> cdot:=diff(c(t),t)=s*(A/2*(x(t))^(-1/2)-d-q(t))**

**> qdot:=diff(q(t),t)=-(q(t)-b1)*(q(t)-b2)+(q(t)-b)*rhs(cdot)/rhs(xdot);**

Computing the "stationary" state (xs,cs,qs):

**> xd:=eval(rhs(xdot),x(t)=x, c(t)=c, q(t)=q):**

**> cd:=eval(rhs(cdot),x(t)=x, c(t)=c, q(t)=q):**

**> qd:=eval(rhs(qdot),x(t)=x, c(t)=c, q(t)=q):**

**> qs:=b: xs:=fsolve(eval(cd,q=qs),x): cs:=fsolve(eval(xd,x=xs),c): ss:=evalf([xs,cs,qs]);**

**> R:=b;**

The Transformed Dynamic System

**> sys := convert([xdot,cdot,qdot],rational):**

**> PDEtools[dchange](x(t)=A^2/4/(r(t)+d)^(2),c(t)=A^2/2/(r(t)+d)-d*A^2/4/(r(t)+d)^2-u(t),q(t)=r(t)-v(t)/s,sys,[u(t),v(t),r(t)]):**

simplify(%) assuming u(t)>0 and v(t)>0 and r(t)>0:

newsys:= expand(solve(%,diff(u(t),t),diff(v(t),t),diff(r(t),t)));

**> eval(newsys, r(t) = R):**

sys2D:= select(has,%, diff);

Eigenvalues of the System

**> M:= [eval(rhs(newsys[2]),u(t)=u,v(t)=v,r(t)=r), **

eval(rhs(newsys[3]),u(t)=u,v(t)=v,r(t)=r),

eval(rhs(newsys[1]),u(t)=u,v(t)=v,r(t)=r)]:

**> us:=1e-20: vs:=1e-20: rs:=R:**

J:= eval(VectorCalculus[Jacobian](M,[u,v,r]),u=us, v=vs, r=rs);

**> LinearAlgebra[CharacteristicPolynomial](J,x):**

**> E,V:= LinearAlgebra[Eigenvectors](J):**

Re(E); Re(V):

Gives one negative eigenvalue and two positive ones.

If I set the stationary value of r (and consequently q) to something other than b, I get one very large eigenvalue and one very small one. Analytically: one is zero, the other is infinite. This led me to rule out all candidate stationary states except r=b. When r=b, it is possible that the terms (r-b)/u and v/u in the jacobian, both of the form 0/0, might tend to something finite. It is this possibility that I investigate here.

A first look at the transformed system

**> DEplot(sys2D,[u(t),v(t)], t=-1..1, u=-0.1..0.1, v=-0.1..0.1,arrows=medium); **

Computing the nullclines (I learned these tricks from Robert Israel)

**> up := subs(sys2D, u(t)=u, v(t)=v, diff(u(t),t));**

vp := subs(sys2D, u(t)=u, v(t)=v, diff(v(t),t));

**> Uu := solve(up, v);**

**> Uv := solve(vp, v);**

**> u1 := fsolve(Uu=abs(Uv[1]),u,avoid=u=0); v1:=eval(Uu,u=u1);**

Plotting the Phase Diagram with the NullClines

**> p1 := DEplot(sys2D,[u(t),v(t)],t=-1..1,u=-0.1 .. 0.1,v=-0.002..0.002, arrows=medium):**

**> p2 := plot([Uu,Uv], u=-0.1 .. 0.1, colour=[brown,blue], thickness=3):**

**> display(p1,p2);**

Zooming on the critical point

**> p1 := DEplot(sys2D,[u(t),v(t)],t=-1..1, u=-0.015 .. 0.01, v = -0.001 .. 0.001, arrows=medium):**

**> p2 := plot([Uu,Uv], u=-0.015 .. 0.01, colour=[brown,blue], thickness=3):**

**> display(p1,p2);**

The critical point can be approached from the north-east quadrant only

Simulating the 3-D Transformed system

**> T := 100:**

**> INI := u(T)=1e-15, v(T)=1e-15, r(T)=b:**

**> SYS := udot, vdot, rdot:**

**> VAR := u(t), v(t), r(t):**

**> Digits:=100: interface(displayprecision=5):**

**> SOL := dsolve(SYS, INI, VAR, type=numeric, range=0..T, abserr=1e-25, output=listprocedure);**

Reverting back from the transformed system (u,v,r) to the original system in (x,c,q)

**> r0:=rhs(SOL[2](0)); u0:=rhs(SOL[3](0)); v0:=rhs(SOL[4](0)); **

**> X := (u,v,r)-> A^2/4/(r+d)^(2): **

C := (u,v,r)-> A^2/2/(r+d)-d*A^2/4/(r+d)^2-u:

Q := (u,v,r)-> r-v/s:

Plotting the original system

**> x := t-> X(rhs(SOL[3](t)),rhs(SOL[4](t)),rhs(SOL[2](t))):**

**> c := t-> C(rhs(SOL[3](t)),rhs(SOL[4](t)),rhs(SOL[2](t))):**

**> q := t-> Q(rhs(SOL[3](t)),rhs(SOL[4](t)),rhs(SOL[2](t))):**

**> px := plot(x(t), t=0..T): display(px);**

**> pc := plot(c(t), t=0..T): display(pc);**

**> pq := plot(q(t), t=0..T): display(pq);**

Unfortunately, couldn't plot the implicit relation among the variables. To be fixed.

**> pxq := implicitplot(q(t),x(t), q=q(0)..q(T), x=x(0)..x(T)): display(pxq): #THIS IMPLICITPLOT FAILS**

pxq := implicitplot3d(t,q(t),x(t), t=0..T, q=q(0)..q(T), x=x(0)..x(T)): display(pxq, axes=boxed): #ALSO FAILS

Error, (in plots/implicitplot) invalid input: the following extra unknowns were found in the input expression: t

Error, (in plots:-display) expecting plot structures but received: pxq

Plotting the transformed system

**> pu := odeplot(SOL, [t,u(t)], 0..T): display(pu);**

**> pv := odeplot(SOL, [t,v(t)], 0..T): display(pv);**

**> pr := odeplot(SOL, [t,r(t)], 0..T): display(pr);**

**> puv := odeplot(SOL, [u(t),v(t)], 0..T): display(puv);**

**> pur := odeplot(SOL, [u(t),r(t)], 0..T): display(pur);**

**> pvr := odeplot(SOL, [v(t),r(t)], 0..T): display(pvr);**

A 3-D plot of the (u,v,r) stable path

**> p3D := odeplot(SOL, [u(t),v(t),r(t)], 0..T, colour=black, thickness=3): display(p3D, axes=boxed);**

This post was generated using the MaplePrimes File Manager

View 9249_mapleprimes08.mw on MapleNet or Download 9249_mapleprimes08.mw

View file details