Matlab

mkdir matlab startup.m addpath '/usr/cld/name/matlab/SIT' -begin 
clear, clf

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

[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;

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

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)];

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)