Question: Approaching a Singular Critical Point in a 3-D Dynamic System: One more attempt

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

Please Wait...