## Approaching a Singular Critical Point in a 3-D Dyn...

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?

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 file details

## Looking for help with an equivalent statement- I t...

My professor wants us to verify the following

(p \/ q) V r = p V (q \/ r)

using the following Maple equivalent statement

with(Logic):

Equivalent(a &and (a &or b),a);

The issue I have is I really do not know that the statement is asking for and where it would go.

Thanks in advance for any help.

## problem with an exponential integral - brunch cut?...

```Int(exp(-1/100+1/1000*I*u)/(4*u^2+1)*exp(1/2*I*u),u = 0 .. 1)

The integrand has positive real and imaginary parts over the
range (just use plot it) and numerical evaluation gives it as
0.53434219089626 + 0.98249392969436e-1 * I (using Digits:=14).

A symbolic integration and using evalf gives the same.

Now writes this as

tstData:=[a=0, b=1, m=1/2, b0 = -1/100, b1 = 1/1000];

J:=Int(exp(b0+b1*u*I)/(4*u^2+1)*exp(u*m*I),u = a .. b);
eval(J,tstData...```

## Can I know the latest prime factorization algorith...

I'm student from Malaysia.

I prefer to study about prime factorization.

Can I know the two latest prime factorization algorithm which are widely use?

pleqse help te toprogram the pade approximant of a function with maple

## scatter plots...

Can anyone tell me how to work out correlation coefficients!!! I need it in very simple terms as I just dont get it ie if 2 variables has a correlation coefficient of .690 what does it mean?  Can someone explain it to me?  With thanks

## intersection of line segment and triangle...

I want to plot in 3D line segment (with 2 point) intersecting triangle(3 d with 3 points).

need help

## matrix and pictures...

Hello

I am a french student and I want to now if maple can transform a picture to a matrix.
thank you

## Gaussian Copula...

Does anyone know how to step-by-step set up a Gaussian Copula in Maple?

I know how to simulate cross correlated random variables by using the Cholesky Decomposition but the requirement that the

correlation matrix must be positive define (all eigenvalues +) is such a pain!

I know that copulas can easily be estimated in MATLAB but I have not seen one in Maple. Any ideas?

## how to change the view of the worksheet?...

when using maple in the worksheet mode, the cursor has the following form
As you can see, cursor shows a vertical line or bracket, to join the group calculations.
I've seen in examples, the cursor has the following configuration.

## Financial Numerics: Estimating Garch(1,1)

by: Maple
```It is some time ago, that I was fighting with that model and how
to estimate the 3 parameters (until recognizing that one may want
restrictions for them).
Essentially one uses the statistics of the underlying data to get
a reasonable starting guess.
For an estimation then one can use the Optimization package.
The gradients involved are best coded as a floating point library,
which here is done through a DLL (it should work on all Windows OS),
code in C is included for that.
Meanwhile all can be done using 'Compiler:-compile' or the Watcom
compiler delivered with Maple.
I left that older sheet as it is - in Maple 9 (may be one should
brush it up for concurrent Maple versions, they are grumbling a bit
about the code) and hope it is of interest despite of that.

Maple sheet: www.mapleprimes.com/files/102_Garch.mws
DLL + code:  www.mapleprimes.com/files/102_Garch.dll.zip
as a pdf:    www.mapleprimes.com/files/102_Garch.mws.pdf

```

## Chemical Engineering...

How do I plot a elliptical cross section with equation (((x-x[0])^2/1.75^2)+(y-y[0])^2/1.35^2))=1? The ellipse is centred at x[0], y[0].

## Help with Maximize function for GARCH(1,1)...

Hi,

I would like to find the optimal values of the parameters  [ w,  a ,  b ]  in a GARCH(1,1)

I however have some problems specifying the Maximize function in the Optimization package

due to all the individual loops.....Any suggestions?

I have attached the workshet and the data

## Maple having difficulty with straight-forward inte...

Hi There,

Here's the problem: I'm trying to tiple-integrate rst e with the limits...
t=r-s...r+s
s=r...∞
r=0...∞

This is a straighforward integral with the use of the gamma function

0 xe dr = Γ(n+1)/a

and my copy of Maple v11 crashes before it yields the answer.