% pclirk.m
% Precision pour le Cycle LImite avec le schema de Runge Kutta
% integration de dx/dt + y = 0,  dy/dt - x = 0 
%                                    par LE Runge Kutta d ordre 4
% fd, iacs, mai, 08 juillet 2004, 25 fevrier 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xlabel ('pclirk. Ordre de convergence du schema de Runge et Kutta  d ordre 4') 
% title ('Cnam Paris. Travaux pratiques Iacs.') ;
m=9
pas=zeros(m,1);
err=zeros(m,1);
%
  n=32;
  for j=1:m;
dt=2*pi/n;
x=zeros(n+1,1);   y=zeros(n+1,1);
px1=zeros(n+1,1); py1=zeros(n+1,1);
px2=zeros(n+1,1); py2=zeros(n+1,1);
px3=zeros(n+1,1); py3=zeros(n+1,1);
px4=zeros(n+1,1); py4=zeros(n+1,1);
x(1:1)=1;
y(1:1)=0;
  for i=1:n;
px1(i) = - y(i)                         ; py1(i) = x(i) ; 
px2(i) = -( y(i) + 0.5 * dt * py1(i) )  ; py2(i) = ( x(i) + 0.5 * dt * px1(i) ) ; 
px3(i) = -( y(i) + 0.5 * dt * py2(i) )  ; py3(i) = ( x(i) + 0.5 * dt * px2(i) ) ; 
px4(i) = -( y(i) + dt * py3(i) )  ;       py4(i) = ( x(i) + dt * px3(i) ) ; 
x(i+1)  = x(i) + dt * ( px1(i) + 2*px2(i) + 2*px3(i)+ px4(i)  ) / 6 ;
y(i+1)  = y(i) + dt * ( py1(i) + 2*py2(i) + 2*py3(i)+ py4(i)  ) / 6 ;
end ;
% 
pas(j)=log(n)/log(2);
err(j)=log ( abs(x(n+1,1)-1) +  abs(y(n+1,1))) / log(2) ;
n=2*n; end;


%  calcul du coefficient avec une droite de regression
alpha=codre(m,pas,err)

plot(pas,err(:,1),'*;logarithme des ereurs;') ; 
axis ([ 4.8   13.2      -48.1     -10.2 ]) 
title('cycle limite    schema de Runge et Kutta : ordre 4')
print -dpdf 21-pclirk.pdf
print -djpg 21-pclirk.jpg
