clear;
// Question 2.1.1
a=1; // pour la bifurcation de Hopf simple
NomFic='vdp05.eps'
TypeDeLigne=-1;
function y=F(x); y=(x^3)/3-a*x; endfunction;
function f=f(x,y);f=(y-F(x));endfunction;
function g=g(x,y);g=-x;endfunction;
// Question 2.1.2
function ww=w(t,v);ww=[f(v(1),v(2));g(v(1),v(2))];endfunction;
function ww=VdP(t,v);ww=[f(v(1),v(2));g(v(1),v(2))];endfunction;
// Question 2.1.3
tMax=50;
tt=0:0.1:tMax;
// 
xMin=-2.5;xMax=2.5;yMin=-3;yMax=+3;
NbTraj=20;
xx=xMin:0.2:xMax;yy=yMin:0.2:yMax;
//fchamp(w,0,xx,yy)
// la cubique
plot2d(xx,F(xx),style=TypeDeLigne,rect=bornes)
//// une trajectoire de VpP(a)
// question 2.1.4
M0=[0;+0.5];
// Question 2.1.5
M=ode(M0,0,tt,w);
// Question 2.1.6
xt=M(1,:);yt=M(2,:);
// Question 2.1.7
xset("window",1);
// Question 2.1.8
bornes=[xMin,yMin,xMax,yMax];
plot2d(xt,yt,style=1,rect=bornes);
// Question 2.1.9
xPlus=max(xt);xMoins=min(xt);
disp(xMoins,"xMoins=",xPlus,"xPlus=");
disp(xPlus-xMoins,"amplitude=")
// Question 2.1.10
M=ode(M0,0,-tt,w); // la trajectoire négative
xt=M(1,:);yt=M(2,:);plot2d(xt,yt);
// Question 2.1.11
function J=jac(x,y);
    acc=0.0000001; //un "petit" accroissement
    J(1,1)=(f(x+acc,y)-f(x,y))/acc;
    J(2,1)=(g(x+acc,y)-g(x,y))/acc;
    J(1,2)=(f(x,y+acc)-f(x,y))/acc;
    J(2,2)=(g(x,y+acc)-g(x,y))/acc;
endfunction;
J=jac(0,0)
    
J
[vcp,vlp]=spec(J)
// Question 2.1.12
//// une autre trajectoire de VpP(a)
M1=[0;-0.5];
M=ode(M1,0,tt,w);
xt=M(1,:);yt=M(2,:);
plot2d(xt,yt,style=1)
M=ode(M1,0,-tt,w); // la trajectoire négative
xt=M(1,:);yt=M(2,:);plot2d(xt,yt)
// Question 2.1.13
// une trajectoire extérieure de VpP(a)
M2=[0;+2.5];tt=0:0.01:tMax;
M=ode(M2,0,tt,w);
xt=M(1,:);yt=M(2,:);plot2d(xt,yt,style=1)
M=ode(M2,0,-tt,w); // la trajectoire négative
xt=M(1,:);yt=M(2,:);plot2d(xt,yt,rect=bornes)
// Question 2.1.14
// une autre trajectoire extérieure de VpP(a)
M3=[0;-2.5];tt=0:0.01:tMax;
M=ode(M3,0,tt,w);
xt=M(1,:);yt=M(2,:);plot2d(xt,yt,style=1)
M=ode(M3,0,-tt,w); // la trajectoire négative
xt=M(1,:);yt=M(2,:);plot2d(xt,yt,rect=bornes)
xs2eps(1,NomFic);// sauvegarde du dessin en fichier eps
// Questions 2.2
// Bifurcation de Hopf
a=0.01;
tMax=5000;
tt=0:0.1:tMax;
xset("window",2);
M0=[0;+0.01];
M=ode(M0,0,tt,w);
xt=M(1,:);yt=M(2,:);
plot2d(xt,yt,style=1)
xPlus=max(xt);xMoins=min(xt);
disp(a,"a=")
disp(xPlus,"xplus=",xMoins,"xmoins=");

// dessin
Xplus=ones(11);Xmoins=ones(11);
for i=0:10 do
    if i==0 then a=0.01;
                tMax=5000;
                tt=0:0.1:tMax;
                M0=[0;+0.01];
                M=ode(M0,0,tt,w);
                xt=M(1,:);yt=M(2,:);
                Xplus(i+1)=max(xt);
                Xmoins(i+1)=min(xt);
                disp(Xplus(i+1),Xmoins(i+1),a,"a=",i);
           else a=0.1*i;
                tMax=50;
                tt=0:0.1:tMax;
                M0=[0;+0.5];
                M=ode(M0,0,tt,w);
                xt=M(1,:);yt=M(2,:);
                Xplus(i+1)=max(xt);
                Xmoins(i+1)=min(xt);
                disp(Xplus(i+1),Xmoins(i+1),a,"a=",i);               
    end
end
aa=zeros(11);
aa(1)=0.01;
for i=1:10 do aa(i+1)=i*0.1;
end
xset("window",3);
plot2d(aa,Xplus,style=-1);
plot2d(aa,Xmoins,style=-1);
xs2eps(3,'Hopf.eps');// sauvegarde du dessin en fichier eps



