%   Script STATION.M
%
%   Zerlegung der erzwungenen Bewegung
%   in Übergangsverhalten und stationäres Verhalten
%
% J. Lunze
% 18.2.1996
% für 2. Auflage: 5.3.1999
% für 3. Auflage: 3.4.2001
% für 4. Auflage: 18.9.2003
% für 5. Auflage: 18.4.2005
% für 6. Auflage: 4.12.2006
% für 7. Auflage: 19.6.2008
% für 8. Auflage: 20.1.2010
% für 9. Auflage: 14.2.2012
% für 10. Auflage: 27.1.2014
% für 12. Auflage: 23.07.2019
%
clear
close all
%%    System 1. Ordnung     xdot = ax + bu, x(0)=0
%                           y = x
%     Eingangsgröße   u(t) = cos wt
%
figure(1);
a=-0.1;
b= 1;
omega=0.4;
T = [0:0.1:50];
Tsize=size(T);
deltat = T(1,Tsize(2));
deltat = deltat/(Tsize(2)-1);
for k=1:Tsize(2)
    t = (k-1)*deltat;
    U(k,1) = 2*cos(omega*t);
    YU(k,1) = 2*a*b/(omega*omega + a*a)*exp(a*t);
    YS(k,1) = 2*b/(omega*omega + a*a)*(omega*sin(omega*t) - a*cos(omega*t));
end;
Y = YU + YS;
subplot(4, 1, 1);
plot(T, U, 'b');
latextitle('Station\"ares und \"Ubergangsverhalten eines PT$$_1$$-Gliedes');
axis([0 50 -2.2 2.2]);
xleer;
yticks([-2 0 2]);
latexylabel('$$u$$');
rotateY;
%
subplot(4, 1, 2);
plot(T, YU, 'b', [0 50],[0 0],'k:');
latexylabel('$$y_{\ddot{u}}$$');
rotateY;
axis([0 50 -1.2 0.2]);
xleer;
yticks([-1 -0.5 0]);
%
subplot(4, 1, 3);
plot(T, YS, 'b', [0 50],[0 0],'k:');
xleer;
latexylabel('$$y_{\rm s}$$');
rotateY;
axis([0 50 -5.5 5.5]);
xleer;
yticks([-5 0 5]);
%
subplot(4, 1, 4);
plot(T, Y, 'b', [0 50],[0 0],'k:');
latexxlabel('$$t$$');
latexylabel('$$y$$');
rotateY;
axis([0 50 -5.5 5.5]);
xticks([0 10 20 30 40 50]);
yticks([-5 0 5]);


%%    System 1. Ordnung     xdot = -x + u, x(0)=0
%                           y = x
%     Eingangsgröße   u(t) = sin wt
%
figure(2);
a=-1;
b= 1;
omega=10;
T = [-2:0.005:10];
Tsize=size(T);
deltat = T(1,Tsize(2));
deltat = deltat/(Tsize(2)-1);
U=[];
YU=[];
YS=[];
k0=400;
for k=1:k0
  U(k,1)=0;
  YU(k,1)=0;
  YS(k,1)=0;
end;
k1=b/(omega*omega+a*a);
k2=omega*k1;
for k=k0:Tsize(2)
    t = (k-k0-1)*deltat;
    U(k,1) = sin(omega*t);
    YU(k,1) = k2*exp(a*t);
    YS(k,1) = -k1*(a*sin(omega*t) + omega*cos(omega*t));
end;
Y = YU + YS;
subplot(8, 1, [1,2]);
tend = Tsize(2)*deltat;
plot(T, U, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latextitle('Station\"ares und \"Ubergangsverhalten eines PT$$_1$$-Gliedes');
xleer;
latexylabel('$$u$$');
rotateY;
axis([-2 tend -1.2 1.2]);
yticks([-1 0 1]);
%
subplot(8, 1, [3,4]);
plot(T, YU, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latexylabel('$$y_{\ddot{u}}$$');
rotateY;
yumax=max(abs(YU));
axis([-2 tend -yumax yumax]);
xleer;
yticks([-0.07 0.07]);
%
subplot(8, 1, [5,6]);
plot(T, YS, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latexylabel('$$y_{\rm s}$$');
rotateY;
ymax = max(Y);
tend = Tsize(2)*deltat;
axis([-2 tend -ymax ymax]);
xleer;
yticks([-0.1 0.1]);

subplot(8, 1, [7,8]);
plot(T, Y, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latexxlabel('$$t$$');
latexylabel('$$y$$');
rotateY;
axis([-2 tend -ymax ymax]);
yticks([-0.1 0.1]);
xticks([-2 0 2 4 6 8 10]);



%%  Stationäres Verhalten bei Resonanz
%   
figure(3)
subplot(3,1,[1,2]);
omega=3;
A=[0  omega;
  -omega 0];
b=[1; 1];
c=[1  1];
d=0;
System=ss(A, b, c, d);
T=linspace(0, 8, 500);
U=sin(omega*T);
Y=lsim(System, U, T, [0; 0]);
plot(T, Y, 'b', [0 10], [0 10], 'k:', [0 10], [0 -10], 'k:', [0 10], [0 0], 'k:');
latexxlabel('$$t$$');
latexylabel('$$y$$');
rotateY;
axis([0 8 -8 8]);
xticks([0 2 4 6 8]);
yticks([-8 -4 0 4 8]);
latextitle('Ausgangsgr\"osse bei Resonanz');


%%  Übergangsverhalten und stationäres Verhalten
%   eines instabilen Systems erster Ordnung
%   
figure(4)
a=+1;
b= 1;
omega=10;
T = [-2:0.005:10];
Tsize=size(T);
deltat = T(1,Tsize(2));
deltat = deltat/(Tsize(2)-1);
U=[];
YU=[];
YS=[];
k0=400;
for k=1:k0
  U(k,1)=0;
  YU(k,1)=0;
  YS(k,1)=0;
end;
k1=b/(omega*omega+a*a);
k2=omega*k1;
for k=k0:Tsize(2)
    t = (k-k0-1)*deltat;
    U(k,1) = sin(omega*t);
    YU(k,1) = k2*exp(a*t);
    YS(k,1) = -k1*(a*sin(omega*t) + omega*cos(omega*t));
end;
Y = YU + YS;
subplot(8, 1, [1,2]);
tend = Tsize(2)*deltat;
plot(T, U, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latextitle('Station\"ares und \"Ubergangsverhalten eines instabilen Systems');
xleer;
latexylabel('$$u$$');
rotateY;
axis([-2 tend -1.2 1.2]);
yticks([-1 0 1]);
%
subplot(8, 1, [3,4]);
plot(T, YU, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latexylabel('$$y_{\ddot{u}}$$');
rotateY;
axis([-2 tend -2 2]);
xleer;
yticks([-2 0 2]);
%
subplot(8, 1, [5,6]);
plot(T, YS, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latexylabel('$$y_{\rm s}$$');
rotateY;
tend = Tsize(2)*deltat;
axis([-2 tend -ymax ymax]);
xleer;
yticks([-0.1 0.1]);
%
subplot(8, 1, [7,8]);
plot(T, Y, 'b', [-2 tend], [0 0],'k:', [-2 0],[0 0], 'b');
latexxlabel('$$t$$');
latexylabel('$$y$$');
rotateY;
axis([-2 tend -2 2]);
yticks([-2 0 2]);
xticks([-2 0 2 4 6 8 10]);

%%  Beispiel für die Erregung "in der Nullstelle"
%   G(s) = (s^2 + 4)/(s+1)(s+2)(s+3)
%   u(t) = sin 2t
%   
figure(5)
z=[1];                            % Nenner gegen Erregung gekürzt
n=poly([-1 -2 -3]);
System2=tf(z, n);
T = [-2:0.005:10];
Tsize=size(T);
deltat = T(1,Tsize(2));
deltat = deltat/(Tsize(2)-1);
U=[];
k0=400;
for k=1:k0
  U(k,1)=0;
end;
for k=k0:Tsize(2)
    t = (k-k0-1)*deltat;
    U(k,1) = sin(2*t);
end;
U1=zeros(1, length(T));        % für Simulation:
U1(1,400)=1/deltat;            % impulsförmige Erregung
Y=lsim(System2, U1, T);
subplot(6,1,[1,2]);
plot(T, U, 'b');
xleer;
latexylabel('$$u$$');
rotateY;
yticks([-1 1]);
axis([-2, 10, -1.2, 1.2]);
%
subplot(6,1,[3,4])
plot(T, Y, 'b');
latexxlabel('$$t$$');
latexylabel('$$y$$');
rotateY;
yticks([0 0.1]);
xticks([-2 0 2 4 6 8 10]);
axis([-2, 10, -0.02 0.1]);


%% Herstellung der Bilddateien *.eps
%j

Epsbildc15('station');
