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)
')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
-----