Matlab

Starting command: matlab &
mkdir matlab startup.m addpath '/usr/cld/name/matlab/SIT' -begin 
clear, clf
% all text right of % is just a comment
Summe.m
% function declaration
function z=Summe(x,y);
z=x+y

1D mapping

  1. Plot of a phase space diagram for the harmonic oscillator.
    x=0:0.123:7*pi; 
    for i=1:10
     plot(i*sin(x),':bd')
     hold on
    end
    xlabel('q_\alpha'); title('phase space diagram')
    
  2. Plot of the orbit and the return map for the logistic map using a self-made function: x=Orbit(r,steps).
  3. Plot an orbit and a cobweb of the Bernoulli map.
  4. Plot the Feigenbaum diagram.
  5. Plot the map, an orbit, the auto correlation function and the power spectrum of the logistic (Bernoulli) map.
    
    Orbit1Dmap.m
    
    % 1D map orbit %%%%%%%
    r=2.79999; itend=100;
    r=1.99999; itend=10000;
    rand('state',sum(100*clock))
    x=rand(1);
    % transition phase
    for i=1:30
    x=x*r*(1-x);
    end;
    % rhs
    x0=0:0.01:1;
    rhs=x*r*(1-x)
    subplot(2,2,1)
    plot(x0,rhs)
    x(1)=x;
    for i=1:itend
    % Logistic map %%%%%%%
    % x(i+1)=x(i)*r*(1-x(i));
    end;
    
    subplot(2,2,2)
    iplot=100;
    if itend < iplot , iplot=itend, end;
    plot(x(1:iplot),'d-');
    xlabel('t')
    ylabel('x')
    title('Logistic map orbit')
    fid=fopen('Orbit1Dmap.dat','w');
     for i=1:length(x)
      fprintf(fid,'%f\n',x(i));
     end
    fclose(fid);
    [c,lags]=xcorr(x-mean(x),10);
    c=c/c(fix(length(c)/2)+1);
    % c=xcf(x,x,10)
    subplot(2,2,3)
    plot(lags,c,'d-');
    title('Autocorrelation function')
    xlabel('\tau')
    ylabel('C')
    subplot(2,2,4)
    Hs=spectrum.welch;
    Fs=1; % sample frequency;
    psd(Hs,x,'Fs',Fs)
    
    
  6. Plot the Lyapunov exponent vs. the control parameter for the logistic map.
MapDerivate.m

function [x]=MapDerivate(x)
% logarithm of the map derivative
% r is a global parameter
global r;
global w;
% logistic map
x=log(abs(r.*(1-(x+x))));

map.m

% 1D map analysis
%
clear   % removes all variables from the workspace
global r;
global w;
w=0.8;
clf;    % deletes from the current figure all graphics objects
%r_min=1.0; r_max=3.5;
% r_min=2.9; r_max=3.6;
 r_min=3.56; r_max=3.58;
 r_min=3.5; r_max=3.6;
 r_min=2.5; r_max=3.5;
r_anz=100;
r=[r_min:(r_max-r_min)/r_anz:(r_max)];
iterations=1000;
transition=10;
x0=2*rand-1;
x0=0.881;
x(1,:)=Map(x0);

% Transition to the attractor
for i=1:transition
x(1,:)=Map(x(1,:));
x(1,:)=Map(x(1,:));
end;

lyap=0;
for i=2:iterations
lyap=lyap+MapDerivate(x(i-1,:));
x(i,:)=Map(x(i-1,:));
end;
lyap=lyap/(iterations-1)/log(2);

subplot(2,1,1)
plot(r,x,'b.','MarkerSize',1)
hold on;
xlabel('r'); ylabel('x');
axis([r_min r_max -1 1]);

subplot(2,1,2)
plot(r,lyap,'b-')
hold on;
xlabel('r'); ylabel('\lambda');
plot([r_min r_max],[0 0],'r:')
y_min=max(-3,min(lyap));
y_max=max(-3,max(lyap));
axis([r_min r_max -3 3]);

% 1D map analysis
%
clear   % removes all variables from the workspace
global r;
global w;
w=0.8;
clf;    % deletes from the current figure all graphics objects
%r_min=1.0; r_max=3.5;
% r_min=2.9; r_max=3.6;
 r_min=3.56; r_max=3.58;
 r_min=3.5; r_max=3.6;
 r_min=2.5; r_max=3.5;
r_anz=100;
r=[r_min:(r_max-r_min)/r_anz:(r_max)];
iterations=1000;
transition=10;
x0=2*rand-1;
x0=0.881;
x(1,:)=Map(x0);


% Transition to the attractor
for i=1:transition
x(1,:)=Map(x(1,:));
end;

lyap=0;
for i=2:iterations
lyap=lyap+MapDerivate(x(i-1,:));
x(i,:)=Map(x(i-1,:));
end;
lyap=lyap/(iterations-1)/log(2);

subplot(2,1,1)
plot(r,x,'b.','MarkerSize',1)
hold on;
xlabel('r'); ylabel('x');
axis([r_min r_max -1 1]);

subplot(2,1,2)
plot(r,lyap,'b-')
hold on;
xlabel('r'); ylabel('\lambda');
plot([r_min r_max],[0 0],'r:')
y_min=max(-3,min(lyap));
y_max=max(-3,max(lyap));
axis([r_min r_max -3 3]);


Direction field, Tangent field

  • Plot the direction field and some trajectories for the FitzHugh-Nagumo equations.
    [x,y]=meshgrid(-2.2:.2:2.2,-2.2:0.2:2.2);
    % harmonic oscillator
    quiver(y,-x); axis equal; axis off;
    % damped harmonic oscillator
    beta=0.1;\\
    quiver(y,-x-beta*y); axis equal; axis off;
    % van der Pol
    quiver(y,(1.-x*x)*y-x); axis equal; axis off;
    

    Eigenvalues, Eigenvectors

  • Characterize the fixed points of the damped harmonic oscillator equations and the van der Pol equations.
    a=[1,2,3;4,5,6;7,8,0]
    b=magic(4)  r=rand(3)
    Determinant=det(a)
    Inverse=inv(a)
     [s,d]=eig(a)
    Eigenvalues=Diagonal of d
    Eigenvectors=Columns of s
    ScalarProduct=s(:,1)'*s(:,2)
    Coefficients in the characteristic polynomial poly(a)
    \ / matrix left or right division
    A system of simultaneous equations can be solved using matrix division:
    A=[3,2,-1; -1 3 2; 1,-1,-1]  b=[10,5,-1]'
    A x = b    x=A\b
    b=[10,5,-1]
    x A = b    x=b/A
    

    ODE solver

    ode1solv.m
    
    % odesolv1 solves the 1D ode
    % ode1 contains the rhs of an 1D ode
    
    t0=0.0; tend=20;
    x01=2000000.55555;
    x01=1.0;
    
    options = odeset('RelTol',1e-7,'AbsTol',1e-8);
    [t,x]=ode45('ode1',[t0,tend],x01,options);
    plot(t,x)
    xlabel('t')
    ylabel('x')
    fid=fopen('X1.dat','w');
     for i=1:length(x)
      fprintf(fid,'%f\n',x(i));
     end
    fclose(fid);
    
    ode1.m 
    
    function yprime=ode1(x,y)
    % ode1(x,y) contains the rhs of an 1D ode
    % rhs ODE solver
    % yprime=-2*x*y;
    yprime=sin(x);
    % a=.01 ; b=100;
    % yprime=-b*x*log(a*x);
    
    ode2solv.m
    
    % odesolv2 solves the 2D ode system
    % ode2 contains the rhs of an 2D ode
    
    x0=0; xend=100;
    y01=0.1; y02=1;
    y01=300; y02=7000;  % Lotka-Voterra
    [x,y]=ode45('ode2',[x0,xend],[y01,y02]);
    subplot(1,2,1)
    plot(x,y)
    xlabel('t')
    ylabel('x, x''')
    subplot(1,2,2)
    plot(y(:,1),y(:,2))
    xlabel('x')
    ylabel('x''')
    fid=fopen('X2.dat','w');
     for i=1:length(x)
      fprintf(fid,'%g %e\n',y(i,1),y(i,2));
     end
    fclose(fid);
    
    ode2.m
    
    function yprime=ode2(x,y)
    
    % ode2(x,y) contains the rhs of an 2D ode system
    % rhs ODE solver
    
    yprime=zeros(2,1);
    
    % harmonic oscillator
    %yprime=[y(2); -y(1)];
    
    % damped harmonic oscillator
    beta=0.1;
    %yprime=[y(2); -y(1)-beta*y(2)];
    
    % van der Pol
    yprime=[y(2);(1. -y(1)^2)*y(2)-y(1)];
    
    % Lotka-Volterra
    a1=0.08; a2=1;b1=0.00001;b2=0.002;
    yprime=[y(1)*(-a1+b1*y(2));(a2 -b2*y(1))*y(2)];
    
    
    
    

    xppaut - http://www.math.pitt.edu/~bard/xpp/xpp.html

    Tutorial

    Tips

    All shortcuts (W F) are lowercase!
    No blanks in *.ode: x'=a*x ! No case sensitivity!
    
    Starting command: xppaut *.ode
    Quit (F Q Y)
    Erase the plotting window (E)
    
    Hardcopy: (V 2) (G F O) (G P)