Home › Forums › Get Help › How do I ...? with Maple

bifurcation

Posted on 2008-08-17 04:59 By carmina (8)
Categories:

How could I plot bifurcation for a nonlinear system using maple? For example I have the following nonlin. syst:
x'=a*b*((x^2+c^2)/(x^2+1))-d*x
y'=a*e*((y^3+f^3)/(x^3+1)*(g*x+1)-y
for a, b c, d, e, f,g be given ; I want to get bifurcation by some parameter, e for instance.
An example: for a=7, b=1.4, c=0.2, d=0.1, e=12, f=0.5 and g=2.8 I obtain only one equilibria point, namely an atractor. Futhermore, for different values of e betwen 0.1 and 2, I always get one stable point, so a hypebolic one. This really means thai I have no bifurcation by e? For plotting bifurcation I should have at least nonhyperbolic critical point? How could I get it? Please help me! Thanks!

‹ Bug in coulditbe with abs? math symbols in plot ›

A try

Comment on 2008-08-17 14:04 from Kanellopoulos (6)

Hi! First of all the second equation is incomplete, a parenthesis missing...

I use DE(plot) for such systems but you need initial or boundary conditions..I tried something like this:

 

>with(DEtools):

>a:=7: b:=1.4: c:=0.2: d:=0.1: e:=2.1; f:=0.5:g:=2.8:

> DEplot([D(x)(t)=a*b*((x(t)^2+c^2)/(x(t)^2+1))-d*x(t)
 ,D(y)(t)=a*e*((y(t)^3+f^3)/(x(t)^3+1))*(g*x(t)+1)-y(t)],
 [x(t),y(t)],t=-5..5,[[x(0)=0,y(0)=0]],stepsize=.05);
 

Giorgos

Thanks Giorgios! First of

Comment on 2008-08-17 14:40 from carmina (8)

Thanks Giorgios! First of all, as you noted, a correction is needed
x'=a*b*((x^2+c^2)/(x^2+1))-d*x
y'=a*e*((y^3+f^3)/(x^3+1)*(g*x+1))-y

initial cond. x(0)=y(0)=0.01 With the above source Maple plots a trajectory. O.K. But how colud I identify the bifurcations of the system by some parameter, say e? If there exist..
I obtained below the phase plane, hope it is correct
> sys:=diff(x(t),t)=7*(0.2)*(((x(t))^2+(0.2)^2)/((x(t))^2+1))-(0.1)*(x(t)),diff(y(t),t)=7*15*(((y(t))^3+(0.5)^3)/(((y(t))^3+1)*((0.1)*x(t)+1)))-y(t):
> p2:=DEplot({sys},[x(t),y(t)],t=-1..20,[[0,1.7,52],[0,5,40],[0,9.5,62],[0,10,30],[0,18,57],[0,18,28],[0,28,45],[0,24,21],[0,36,38],[0,35,17],[0,40,18],[0,40,25],[0,60,10],[0,63,1],[0,64,63]],stepsize=0.1,x=0..50,y=0..80,linestyle=6,
linecolor=black,color=red,arrows=SLIM,title=phase_portrait):
> p2:=subs(THICKNESS(3)=THICKNESS(2),[p2]):
> display(p2,font=[TIMES,8],insequence=true);

What are the functions for plotting or looking for bifurcations for nonlinear systems? Thanks again!
carmina

Robert Israel's picture

Bifurcations

Comment on 2008-08-17 16:39 from Robert Israel (Maple Rating 4 1432)

The right side of the system:

> F := [a*b*((x^2+c^2)/(x^2+1))-d*x, a*e*((y^3+f^3)/(x^3+1)*(g*x+1))-y];

The Jacobian:

> J := VectorCalculus[Jacobian](F, [x,y]);

= RTABLE(45781648,MATRIX([[2*a*b*x/(x^2+1)-2*a*b*(x^2+c^2)/(x^2+1)^2*x-d, 0], [-3*a*e*(y^3+f^3)/(x^3+1)^2*(g*x+1)*x^2+a*e*(y^3+f^3)/(x^3+1)*g, 3*a*e*y^2/(x^3+1)*(g*x+1)-1]]),Matrix)

A bifurcation can occur when the real part of an eigenvalue of the Jacobian at an equilibrium point crosses 0.  In this case the Jacobian is lower triangular, so the eigenvalues are the diagonal elements (and are real).

> S:= solve([op(F), J[2,2]], [x,y,e]);

 =[[x = RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z), y = RootOf(2*_Z^3-1,label = _L1)*f, e = 2/3*(d+a*b*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z)^2+a*b*c^2-d*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z))/d/a/f^2*RootOf(2*_Z^3-1,label = _L1)/(g*RootOf(-a*b*_Z^2-a*b*c^2+d*_Z^3+d*_Z)+1)]]

 Try it for some values of the other parameters:

 > pars:= { a=7, b = 0.2, c=0.2,d=0.1,f=1,g=0.1};

> evalf([allvalues(eval(S,pars))]);

[[[x = 13.93110365, y = .7937005260, e = 85.43224375]], [[x = .3444817452e-1+.1975123882*I, y = .7937005260, e = .7499096375e-1-.2003550007e-2*I]], [[x = .3444817452e-1-.1975123882*I, y = .7937005260, e = .7499096375e-1+.2003550007e-2*I]], [[x = 13.93110365, y = -.3968502630+.6873648188*I, e = -42.71612188+73.98649342*I]], [[x = .3444817452e-1+.1975123882*I, y = -.3968502630+.6873648188*I, e = -.3576035666e-1+.6594585468e-1*I]], [[x = .3444817452e-1-.1975123882*I, y = -.3968502630+.6873648188*I, e = -.3923060707e-1+.6394230467e-1*I]], [[x = 13.93110365, y = -.3968502630-.6873648188*I, e = -42.71612188-73.98649342*I]], [[x = .3444817452e-1+.1975123882*I, y = -.3968502630-.6873648188*I, e = -.3923060707e-1-.6394230467e-1*I]], [[x = .3444817452e-1-.1975123882*I, y = -.3968502630-.6873648188*I, e = -.3576035666e-1-.6594585468e-1*I]]]

We're only interested in real values of x,y and e, so the bifurcation should be at e= 85.43224375.

For e = 85.4:

> F854:= eval(F,{op(pars),e=85.4});
   with(plots): with(DEtools):
   display([implicitplot(F854,x=13.92 .. 13.95,y=0.7 ..0.9,gridrefine=2,colour=[green,blue]),
   DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(F854,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92 .. 13.95, y=0.7 .. 0.9)]);

Here we see two equilibrium points, one a saddle and the other a stable node.

For e = 85.5:

> F855:= eval(F,{op(pars),e=85.5});
  display([implicitplot(F855,x=13.92 .. 13.95,y=0.7 ..0.9,gridrefine=2,colour=[green,blue]),
  DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(F855,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92 .. 13.95, y=0.7 .. 0.9)]);

The two equilibrium points have collided and disappeared.

Honestly i dont know if

Comment on 2008-08-17 17:22 from Kanellopoulos (6)

Honestly i dont know if maple  has fixed commands  for  making  bifurcations...If i were you i would try to find the stable points for many different values of e (with all the others fixed) and then make a pointplot...As in Logistic map...Also your system seems a bit strange to me, the fisrt equation is only dependent of x(t) so you can solve the first by its own, and then you will have only one ode....I cant think anything else...

Good luck!

Giorgos

 

 

Robert, there were of very

Comment on 2008-08-17 19:12 from carmina (8)

Robert, there were of very help for me your explanations, but I have some questions. I'm a beginner in Maple, so I'm reading all commands step by step (working with Maple 9.5).
I got an error massage after the command
> S:= solve([op(F), J[2,2]], [x,y,e]);
I understand that you compute J[2,2] (a proper value) at the equilibrium point. OK, I deduce I should try the same algorithm for J[1,1], isn't it? Now remain only two questions: why doesn't work the above command?maybe because working with Maple 9.5 :(
and what means _Z? Thanks again you all!
carmina

Robert Israel's picture

Maple 9.5

Comment on 2008-08-18 00:09 from Robert Israel (Maple Rating 4 1432)

Maple 9.5 needs sets rather than lists in solve.  So it should be

 S:= solve({op(F), J[2,2]}, {x,y,e});

Also, Maple 9.5 does not allow the gridrefine option in implicitplot.

Using J[1,1] will not work, because F[1] and J[1,1] do not depend on y or e.

_Z is used as the variable in RootOf notation.  That is, RootOf(f(_Z)) stands for any value _Z such that f(_Z) = 0.

It worked! Some questions 1.

Comment on 2008-08-18 06:34 from carmina (8)

It worked! Some questions
1. Where do appear the initial conditions x(0)=y(0)=0.01?
2. the bifurcation is considered at e=85.4 or at e=85.5? Where should I stop, where are still two equillibria points or at least one or no more one?

Well, I tried for my real system (actually I thought is a slight modification but now I think it is significantly modified)
So
x'=a*b*((x^2+c^2)/(x^2+1))-d*x
y'=a*e*((y^3+f^3)/(y^3+1)*(g*x+1))-y (a modification appears in the down term, for the first bracket I have y instead of x ) I followed the above approach with the following changes
pars:= { a=7, b = 0.2, c=0.2, d=0.1, f=0.5, g=2.8};
and I got
>evalf([allvalues(eval(S,pars))]);

[{e = -0.1155507917e-1-.2044789251*I, x = 0.3444817452e-1-.1975123882*I, y = -.5608636000-.9714442513*I},
{e = -.1713064041-.1122464547*I, x = 0.3444817452e-1+.1975123882*I, y = -.5608636000-.9714442513*I},
{e = -0.3143281871e-2-0.5444323904e-2*I, x = 13.93110365, y = -.5608636000-.9714442513*I},
{e = -.1713064041+.1122464547*I, y = -.5608636000+.9714442513*I, x = 0.3444817452e-1-.1975123882*I},
{e = -0.1155507917e-1+.2044789251*I, x = 0.3444817452e-1+.1975123882*I, y = -.5608636000+.9714442513*I},
{x = 13.93110365, y = -.5608636000+.9714442513*I, e = -0.3143281871e-2+0.5444323904e-2*I}, {e= .1828614833+0.9223247044e-1*I, y = 1.121727200, x = 0.3444817452e-1-.1975123882*I},
{e = .1828614833-0.9223247044e-1*I, x = 0.3444817452e-1+.1975123882*I, y = 1.121727200},
{e = 0.006286563743, x = 13.93110365, y = 1.121727200},
{e = -0.1491201246e-1-.2638832880*I, y = -.2228705874-.3860231810*I, x =.3444817452e-1-.1975123882*I},
{x = 0.3444817452e-1+.1975123882*I, e = -.2210736249-.1448558256*I, y = -.2228705874-.3860231810*I},
{x = 13.93110365, y = -.2228705874-.3860231810*I, e = -0.4056454986e-2-0.7025986134e-2*I},
{e = -.2210736249+.1448558256*I, y = -.2228705874+.3860231810*I, x = .3444817452e-1-.1975123882*I},
{x = 0.3444817452e-1+.1975123882*I, e = -0.1491201246e-1+.2638832880*I, y = .2228705874+.3860231810*I},
{x = 13.93110365, y = -.2228705874+.3860231810*I, e = -0.4056454986e-2+0.7025986134e-2*I},
{e = .2359856373+.1190274624*I, y = .4457411749, x = 0.3444817452e-1-.1975123882*I},
{e= .2359856373-.1190274624*I, x = 0.3444817452e-1+.1975123882*I, y = .4457411749},
{e = 0.008112909971, x = 13.93110365, y = .4457411749}]

In this case I have two possibilities:
{e = 0.006286563743, x = 13.93110365, y = 1.121727200},
{e = 0.008112909971, x = 13.93110365, y = .4457411749}]

for e=0.006 I obtain no equilibria points. With
display([implicitplot(e006,x=13.92 .. 13.96,y=1.1..1.3,colour=[green,blue]),
DEplot(zip(`=`,[D(x)(t),D(y)(t)],eval(D006,{x=x(t),y=y(t)})),[x(t),y(t)],t=1..20,x=13.92..13.96, y=1.1..1.3,arrows=slim)]);
I can see the green line; If I take x=13.92..13.95 no line I can see, how do I interprete this fact?

for e=0.008 again I do not obtain any equilibria point, but for e=0.0081, I do; what is the conclusion? there is no bifurcation? :(( Where I'm wrong?

Thanks for your patience!

Robert Israel's picture

Some questions

Comment on 2008-08-18 15:18 from Robert Israel (Maple Rating 4 1432)

1) I did not use the initial conditions, because in my understanding of bifurcation theory initial conditions are not important.  Bifurcations occur when equilibria (or sometimes other invariant sets) change their character. 

2) The bifurcation occurs at approximately e = 85.43224375.  In order to see what happens there, I looked at one value of e on each side of this point. 

 

I'm back. Trying to look for

Comment on 2008-08-18 17:35 from carmina (8)

I'm back.
1. Trying to look for bifurcations by a, for the set of values
pars:= {F2 =1, E2=0.2, D2=0.5, D1=14, E1=0.15, F1=2.2} Using the command
> evalf([allvalues(eval(S,pars))]) I get the the following Error message, (in evala/preproc3) floats not handled yet. What's wrong?
2. Is it posiibly to find two bifurcations points for the same parameter?

Robert Israel's picture

Bug

Comment on 2008-08-19 02:54 from Robert Israel (Maple Rating 4 1432)

1)  It looks like a bug.  Can you give us the exact commands you are using?
I don't know what are F2, E2 etc.

2) Yes, it should be possible.

bug?

Comment on 2008-08-19 03:21 from carmina (8)

Sorry,
pars:= {b =1, c=0.2, d=0.5, e=14, f=0.15, g=2.2} (my original system has coefficients like E2, E2..etc) I studied the system using the commands you proposed above; it worked for e, g but as I see not for a.

Robert Israel's picture

bug?

Comment on 2008-08-19 15:32 from Robert Israel (Maple Rating 4 1432)

I didn't get that bug in Maple 9.5.  Perhaps you could upload your worksheet.  I did get something else interesting: two spurious solutions with y=0, which clearly don't satisfy J[2,2] = 0.  Disregarding those, I get three bifurcations:

[{a = -.9628740975, y = .1190550790, x = -.8162589472}, {a = 1.249060056, x = .1592138727, y = .1190550789}, {y = .1190550789, a = .9588448130, x = .8747161127}]

Here's my worksheet.

Download 4541_bifurcation.mws
View file details

here is

Comment on 2008-08-19 18:33 from carmina (8)

Ok. I apologize for the mistakes along my posts. In the attach you can see my real system. Another question: after this computation, how could I do to insert the code into another file (in latex for example, but not exporting the document, I would like with copy/paste) just as it is.  More precisely i want to appear 0.0012 not something times with "e", the Jacobian as a matrix and so on. Thanks a lot!

p.s I have Maple 9.5 and start a sheet running Maple 9.5 not Classic Worksheet 9.5. I note that it runs rather slowly.

View 9134_bifurcAcaseII.mw on MapleNet or Download 9134_bifurcAcaseII.mw
View file details

Robert Israel's picture

bug

Comment on 2008-08-19 23:39 from Robert Israel (Maple Rating 4 1432)

The problem (which still exists in Maple 12) seems to have to do with the nested RootOfs.  For a simpler example:

> Q:= RootOf((p-RootOf(_Z^3-1))*_Z^3-1);

> eval(Q, p = 2.0);

Error, (in evala/preproc3) floats not handled yet

Technically, this may be a weakness rather than a bug.  I think there are some difficult mathematical issues connected with the interaction between algebraic numbers and floats.

A work-around is to either use exact rationals instead of floats (e.g 22/10 instead of 2.2), or to substitute the parameter values before trying solve.  The second is probably better.  Thus:

> S:=solve(eval({op(F),J[2,2]},pars),{x,y,A});
    remove(has,{S},I);

{{A = -.9791368510,x=-0.7190485833,y=0.1193582308},{A=0.1380471481,x=0.01107628534,y=1.256721040}}

 

 

system

Comment on 2008-08-20 11:19 from carmina (8)

Indeed, I have obtained the same values following the hint above. Further it would be possibly to look for bifurcations by a parameter for nonlinear and non autonomous systems? for example if I would substitute A by a function   A(t) or at least a function constant on intervals and look again for bifurcations by F1 (here g). Or more, bifurcations by two parameters at once ? Should here be a meaning?

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}