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')
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
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')
% 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)
% 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)
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,:));
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]);
[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)