Mutual information

Calculate the Shannon entropy, Renyi entropy (q=0.2 and q=2), mutual information (matlab: mi(x,10,40) ODER gnuplot and TISEAN: plot"<./mutual -b50 -D30 Random.dat"w lp) of a random telegraph sample [floor(rand(N,1)+0.5)], a periodic signal [x=floor((sin(2*pi/2.1*t)+1));] and trajectories of the logistic map for different control parameters depending on the sample size N. Useful matlab script: symbolics.m.

Modell-Raten

1. Gegeben sei folgende Komponente eines Orbits eines dynamischen Systems: H1. Die Zeit zwischen zwei Messungen betrage dt=0.03 Sekunden. Charakerisieren Sie das dynamische System.

2. Gegeben sei folgende Komponente eines Orbits eines dynamischen Systems: V2. Die Zeit zwischen zwei Messungen betrage dt=1 Sekunde. Charakerisieren Sie das dynamische System.

Nutzen Sie matlab oder TISEAN.
Einlesen einer Datei mit dem Namen V2.dat und Ablage der enthaltenen Daten als Feld x in matlab: x=importdata('V2.dat')

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

Dissipative discrete deterministic dynamical system, state space, forward orbits, cobweb, steady state, fixed point, cycle of period 1, cycles of period n, aperiodic behavior, locally unstable/neutral/indifferent/stable/super stable fixed points, control parameter, stability chart, bifurcation diagram. Statistical & ergodic properties of orbits. Density function, acf, power spectrum.
  1. Plot of a phase space diagram for the harmonic oscillator.
    clf,clear
    x=0:0.123:7*pi; 
    for i=1:10
     plot(i*cos(x),i*sin(x),':bd')
     hold on
    end
    xlabel('q_\alpha'); title('phase space diagram')
    
  2. Plot of orbits/trajectories and return maps/scatter plots of Gaussian noise process, a harmonic process and the logistic map.
  3. Plot the density estimations of orbits/trajectories.
    % density.m
    clear,clf
    N=500;
    x=randn(N,1);
    subplot(3,1,1)
    plot(x,'-d')
    x_max=max(x),x_min=min(x)
    classes=10
    class_range=x_max-x_min
    i=0,p=(1:classes+1)'*0;
    while (i < N)
        i=i+1;
        class=1+fix((x(i)-x_min)*classes/class_range);
        p(class)=p(class)+1;
    end
    subplot(3,1,2)
    bar(p)
    %axis('tight')
    subplot(3,1,3)
    hist(x,11)
    
  4. Plot the map, an orbit, the auto correlation function and the power spectrum of the logistic (Bernoulli) map.
    
    % Orbit1Dmap.m
    clear,clf
    % 1D map orbit %%%%%%%
    r=1.99999; itend=10000;
    r=3.78;
    disp(['Control parameter=',num2str(r)])
    disp(['Orbit length     =',num2str(itend)])
    %rand('state',sum(100*clock))
    xi=rand;
    disp(['Initial point    =',num2str(xi)])
    % transition phase
    for i=1:30
        % xi=mod(r*xi,1.);
        xi=xi*r*(1-xi);
    end;
    % rhs of the map
    x0=0:0.01:1;
    rhs=x0.*r.*(1-x0);
    % rhs=mod(r*x0,1.);
    subplot(2,2,1) % plot the rhs
    plot(x0,rhs),xlabel('x'),ylabel('f(x)')
    axis square
    title('rhs of the map')
    % preallocation saves time
    x=zeros(itend,1);
    x(1)=xi;
    disp(['Initial point after transient phase=',num2str(xi)])
    for i=1:itend-1
        % Bernoulli
        % x(i+1)=mod(r*x(i),1.);
        % Logistic map %%%%%%%
        x(i+1)=x(i)*r*(1-x(i));
    end;
    
    subplot(2,2,2) % Orbit of the map
    iplot=100;
    if itend < iplot , iplot=itend, end;
    plot(x(1:iplot),'d-');
    xlabel('t')
    ylabel('x')
    title('Logistic map orbit')
    % title('Bernoulli map orbit')
    
    % Write the orbit as ASCII column
    fid=fopen('Orbit1Dmap.dat','w'); % File allocation
     for i=1:length(x)
      fprintf(fid,'%f\n',x(i));
     end
    fclose(fid);
    % Remove the mean
    x=x-mean(x);
    % Fix the lag maximum of the acf
    lagmax=10;
    [c,lags]=xcorr(x,lagmax);  % Calculate the acf
    c=c/c(fix(length(c)/2)+1); % Normalize the acf
    % c=xcf(x,x,10)
    sig=lags*0+1.96/sqrt(itend); % Significance level of the acf
    zero=lags*0;                 % Zero line of the acf
    
    subplot(2,2,3) % acf
    plot(lags,c,'d-',lags,zero,'-');
    axis([ 0 lagmax -1.2 1.2])
    % xis('tight')
    hold on
    errorbar(lags,zero,sig)
    
    title('Autocorrelation function')
    xlabel('\tau')
    ylabel('C')
    subplot(2,2,4) % Power spectrum
    Hs=spectrum.welch;
    Fs=1; % sample frequency;
    psd(Hs,x,'Fs',Fs)
    
    
    
  5. Plot the Feigenbaum diagram.
    MapDerivate.m
    
    function [x]=MapDerivate(x)
    % logarithm of the map derivative
    % r is a global parameter
    global r;
    % logistic map
    x=log(abs(r.*(1-(x+x))));
    
    map.m
    
    % 1D map analysis
    %
    clear   				% removes all variables from the workspace
    global r; 				% control parameter
    clf;    				% deletes from the current figure all graphics objects
    r_min=3; r_max=4.;			% range of the control parameter
    r_anz=100;				% number of control parameters-1
    r=[r_min:(r_max-r_min)/r_anz:(r_max)];	% control parameter set
    iterations=1000;			% orbit length
    transition=10;				% length of the transient orbit
    x0=rand;				% randomly chosen initial point
    x=zeros(iterations,r_anz+1);		% time saving preallocation
    x(1,:)=Map(x0);
    
    % Transition to the attractor
    for i=1:transition
    x(1,:)=Map(x(1,:));
    end;
    
    lyap=MapDerivate(x(1,:));
    for i=2:iterations
    x(i,:)=Map(x(i-1,:));
    
  6. Plot the Lyapunov exponent vs. the control parameter for the logistic map.
    lyap=lyap+MapDerivate(x(i-1,:));
    end;
    lyap=lyap/(iterations)/log(2);
    
    subplot(2,1,1)
    plot(r,x,'b.','MarkerSize',1)
    hold on;
    xlabel('r'); ylabel('x');
    
    subplot(2,1,2)
    plot(r,lyap,'b-')
    hold on;
    xlabel('r'); ylabel('\lambda');
    plot([r_min r_max],[0 0],'r:')
    axis([r_min r_max -3 1]);
    
    

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)
    
    

    C++

    // g++ logmap.cpp
    // gnuplot
    // set xlabel "r"
    // plot "<./a.out" w d
    // set term postscript landscape 'Helvetica' 14
    // set output 'name.ps'
    // plot "<./a.out" w d
    // set term x11
    
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using namespace::std;
    int main ()
    {
    double x,r;
    int i;
    x=0.6;
    r=2.5;
    //ofstream Datei ("Feigenbaum.dat");
    for (r=2.5;r<=4;r=r+0.0025)
            {for (i=0;i<200;i++)
                    {x=r*x*(1-x);
                    //Datei<< r <<"\t"<< x<