% vdpco3.m
% integration de     dx/dt =  y + x - x**3 , dy/dt =  -x
% par les schemas d Euler expli, Euler impli,
%       Crank-Nicolson,  Heun et Runge Kutta
% la solution exacte est l attracteur de Van der Pol
% fd, iacs, 18 septembre 2004, 2, 16 janvier 2005, 22 janvier 2006.
% fd, tp pour la page web, 25 fevrier, 02 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlabel ('vdpco3. Cycle limite non lineaire avec cinq  schemas.') ;
title ('Oscillateur de Van der Pol') ;
n=4;  dt=1/n
xe=zeros(2,1);
xe(1,1)=0.1  ;  xe(2,1)=0 ; xi=xe ; xcn=xe ;  xhe=xe ; xrk=xe ;

niter = 6*n 
traje=zeros(niter+1,2) ; traje(1,:) = xe ; 
traji=zeros(niter+1,2) ; traji(1,:) = xi ; 
trajn=zeros(niter+1,2) ; trajn(1,:) = xcn ; 
trajr=zeros(niter+1,2) ; trajr(1,:) = xrk ;  
trajh=zeros(niter+1,2) ; trajh(1,:) = xhe ;  


hold off 
plot ( xe(1,1) ,  xe(2,1) , '*k;depart;')  ; hold on

	   
for i=1:niter ;
xne  = xe + dt * f(xe) ;   xni  = er(xi,dt) ; xncn  = cn(xcn,dt) ;
xnrk  = rk(xrk,dt) ;  xnhe  = he(xhe,dt) ;

traje(i+1,:) = xne ; traji(i+1,:) = xni  ; trajn(i+1,:) = xncn  ;
trajr(i+1,:) = xnrk  ; trajh(i+1,:) = xnhe  ;

xe=xne ; xi=xni  ;  xcn=xncn   ; xrk=xnrk ; xhe = xnhe ; 
end ;


plot ( traje(:,1) , traje(:,2), '*c;Euler explicite;') ; hold on 
plot ( traji(:,1) , traji(:,2), '*m;Euler implicite;') ; hold on 
plot ( trajn(:,1) , trajn(:,2), 'ob;Crank Nicolson;') ; hold on 
plot ( trajh(:,1) , trajh(:,2), 'dr;Heun;') ; hold on 
plot ( trajr(:,1) , trajr(:,2), '+k;Runge Kutta;') ; hold on 

axis ([ -.8   .7   -.5   1.4 ]) 
  title ('Oscillateur de Van der Pol : ecart entre Euler explicite, Euler implicite et Crank-Nicolson') ;

hold off 

%% print -dpdf 08-vdpco3.pdf
%% print -djpg 08-vdpco3.jpg
print('08-vdpco3.jpg', '-djpg', '-S891,630');




%%% veche(1,1)=xhe(1,1) ; veche(2,1)=xnhe(1,1) ; 
%%% veche(1,2)=xhe(2,1) ; veche(2,2)=xnhe(2,1) ; 
%%% xe=xne ; xi=xni ;  xcn=xncn   ;  xhe=xnhe   ;    xrk=xnrk   ;  end ; 
