Bonjour,
J'ai recu plusieurs codes écrits en Matlab, et je voudrais les convertir en Mathematica. Je ne sais pas me servir de Mathlab et mon établissement ne le possede pas donc je dois absolument les convertir en Mathematica. Est ce que quelqu'un peut m'aider ? Merci d'avance
Voici les codes:
(1er code)
(Defining the constant surface for the Poincare section
Zerocross2.m)
function [v,ist,df] = zerocross2(t,x,varargin)% varargin absorbs unwanted
parameters flag,g,l1, etc
% event function for test_poincare.m
v = sin(x(2)); % event is x(2)=0 (ie when v is zero)
%v = x(2);
ist = 1; % if true, terminate evolution when this event occurs
df = 1; % increasing sense only
(2ème code)(Producing the trajectory of the outer bob, Poincare section, angular velocity graph and energy of the system
PoincareTrajecAngleEnergy.m)
clear
%Code mainly adapted from Alexander Barnett
%declare time step
T = 10;
%y0 = [0;0;-2;0]; % close to regular
%y0 = [pi;pi;0;0];
% y0 = [pi;pi;.5;0]; % chaotic
% y0=[0.5233;0;0.5233;0]%periodic
y0 = [0.2,0.2828,0,0]%perfect periodic with energy = 0.7809
%y0 = [0.2,-0.2828,0,0]%perfect quasiperiodic with energy = 0.7809
%y0 = [27.8;0;1.22;2.62] %perfect KAM scenario
%y0 = [0,0,2,20]
%y0 = [0,0,1,20];
y0 =[pi;0;.5;.5];%high energy
s=ode45(@doublependulum2,[0,T],y0,[],0,1,1,1,1);
t = 0:0.01:T; x = deval(s,t)';
x1=sin(x(:,1));
y1=-cos(x(:,1));
x2=x1+sin(x(:,2));
y2=y1-cos(x(:,2));
figure; plot(x2,y2); axis equal; xlabel('position');ylabel('vel ocity');title
('trajectory')
tmax = 1e2; % max time to wait until next intersection
ns = 1000; % how many intersection (iterations of P map)
tp = nan*(1:ns); yp = nan*zeros(numel(y0),ns); yi = y0; % init arrays
figure;
for n=1:ns
s = ode45(@doublependulum2, [0, tmax], yi,
odeset('Events',@zerocross2,'a bstol',1e-9),0,9.8,1,1,1,1);
if isempty(s.xe), disp('no intersection found!');
else tp(n) = s.xe(end); yi = s.ye(:,end); yp(:,n) = yi;
%plot3(mod(yi(1),2*pi),yi(3),y i(4),'+'); hold on; if mod(n,10)==0, drawnow; end
% note it's in 3d now
plot(yi(1),yi(4));hold on; if mod(n,10)==0, drawnow; end
end
end
xlabel('x_1'); ylabel('x_3'); zlabel('x_4'); axis vis3d;
title('Poincare section (sin(x_2)=0, increasing)');
figure;
plot(t,x(:,3),'magenta',t,x(:, 4),'black');title('Angular velocities');
energy = hamiltonian(y0)
(3ème code)(Calculates the lyapunov exponent, the average lyapunov exponent and traces the separation of close trajectories
Crudelyap.m)
clear
nruns =100;
%for calculating average later
h = zeros(1,nruns);
%indexing of h
k = 1;
for j=1:nruns
% y0 = [pi+0.1*rand(1);pi;.5;0];
% y0 = [0+0.1*rand(1);0;-2;0];
%y0 = [pi;pi;10.515;-2.17];%periodic with T = 0.5
%y0 = [];%quasiperiodic;
y0 = [pi+0.1*rand(1);pi;0.5;0];%chaotic
%y0 = [pi;pi;0;0];
%y0 = [0;0;-2;0];%close to regular
%y0 = [27.8;0;1.22;2.62];
%y0 = [pi;0;0.5;0.5];
%y0 = [0;0;2;20];
%y0 = [pi;0;5;5];%high energy
%y0 = [0.2;0.2828;0;0];%perfect periodic with energy = 0.7809
%y0 = [0.2+0.1*rand(1);-0.2828;0;0];
%y0 = [0.2;-0.2828;0;0];%perfect quasiperiodic with energy = 0.7809
T = 20;
o = odeset('abstol',1e-14);
s=ode45(@doublependulum2,[0,T],y0,o,0,1,1,1,1); % note changes
eps = 1e-9;
s2 = ode45(@doublependulum2,[0,T],y0+[eps;0;0;0],o,0,1,1,1,1); % note changes
t = 0:0.01:T; x = deval(s,t)';xe = deval(s2,t)';
x1=sin(x(:,1));
y1=-cos(x(:,1));
x2=x1+sin(x(:,2));
y2=y1-cos(x(:,2));
xe1 = sin(xe(:,1));
ye1=-cos(xe(:,1));
xe2 = xe1+sin(xe(:,2));
ye2= ye1-cos(xe(:,2));
%store the size
N = numel(ye2);
%figure; plot(t, [y2 ye2], '-');title('Trajectories of two intial points with a
distance 10^-9')
logsep = zeros(1,N);
%loop fills matrix logsep with natural log of the changing separation
%between the initial trajectories
for i = 1:N
logsep(i) = (log(abs(y2(i)-ye2(i)))-log(eps));
end
%figure; plot(t,logsep);title('lyapunov exponent');
t1=5;
t2=0.7;
%stores the exponent in a matric
h(k) = (logsep(t1/0.01)-logsep(t2/0.01))/(t1-t2);
k=k+1;
end
figure;
%average lyapunov exponent
plot(1:j,h)
title('average lyapunov exponent');
(4ème code)
clear
T = 10;
%y0 = [0;0;0;10]; %really good
%y0 = [pi;pi;0.5;0];
%y0 = [0;0;-2;0];
%y0=[0;0;2;20];
%y0 = [27.8;0;1.22;2.62] %perfect KAM scenario
y0 = [0.2,-0.2828,0,0];
y0 = [pi;0;.5;.5]
%y0 = [pi;pi;0.5;0]
s=ode45(@doublependulum2,[0,T],y0,[],0,1,1,1,1);
t = 0:0.01:T; x = deval(s,t)';
%find energy for the inital conditions
energy = hamiltonian(y0)
%creates a plot with four figures of the variables plotted against each
%other
figure;
subplot(222)
plot(x(:,3),x(:,4));title('ang ular velocities')
subplot(221)
plot(x(:,1),x(:,2));title('ang les')
subplot(223)
plot(x(:,1),x(:,3));title('pos itions and momentum of 1')
subplot(224)
plot(x(:,2),x(:,4));title('pos itions and momentum of 1')
-----