//Hier ist das Scilab-Programm, das die Differentialgleichung numerisch mit dem Euler-Integrationsverfahren löst:

// Set the initial condition
x0 = 1;

// Set the number of integration steps and the step size
n = 10;
dt = 0.1;

// Create arrays to store the time and solution values
t = 0:dt:dt*n;
x = zeros(1, n+1);
x(1) = x0;

// Perform the Euler integration
for i = 1:n
    x(i+1) = x(i) + dt * x(i);
end

// Plot the results
plot(t, x)
xtitle('Numerical Solution of dx/dt = x', 'Time (s)', 'x')
//Das Programm gibt das Ergebnis als Grafik aus, die die numerische Lösung von x(t) gegen die Zeit t darstellt.

Code 0-1: Antwort von ChatGPT

A=[1,2;3,4];
b=[19;43];

xy = inv(A)*b

Code 0-2: Lösung mit Hilfe von Scilab.

clear;

C=1; //  N/M
m=1; //  kg
x0=1; // m
v0=0; // m/s

dt=0.01; // s
t = 0;
x = x0;
v = v0;

x_merk = x0;
v_merk = v0;
t_merk = t;
n=1;
for t=0:dt:5

    x_neu = x + v*dt;
    v_neu = v + (-(C/m)*x)*dt;

    x = x_neu;
    v = v_neu;

    n=n+1;
    t=t+dt;
    x_merk(1,n)=x;
    v_merk(1,n)=v;
    t_merk(1,n)=t;
end

plot(t_merk,x_merk,t_merk,v_merk);

Code 0-3: Feder-Masse-Schwinger mit Eulerintegrationsverfahren gelöst.

x0 = 1;

dt=0.0001;

x=x0;

x_merk = x;
t_merk = 0;
n=1;

for t=0:dt:5

   x_neu = x + x*dt;

   x = x_neu;

   n=n+1;
   x_merk(1,n)=x;
   t_merk(1,n)=t;
end

y = exp(t_merk);

plot(t_merk,x_merk,'red',t_merk,y,'blu--');

Code 0-4: Musterlösung zu Aufgabe 3.

Aufgabe1:
Simulieren Sie die resultierende dynamische Euler-Gleichung mit ODE in Scilab.
Verwenden Sie folgende Parameter:
l=1m (h=0,5m)
m=1kg
r=0,1m
phi_null = 10 Grad, 30 Grad, 90 Grad (als Tests), phi_punkt_null=0 rad/s 

Aufgabe2:
Ein Antriebsmoment MA zusätzlich in die dynamischen Gleichungen des Körperpendels einbauen und
die gesamte Herleitung damit nochmal anpassen, bis man wieder die Eulergleichung ohne Unbekannte hat.

Aufgabe3:
Bestimmen Sie die Eigenwerte der linearisierten dynamischen Gleichungen des Körperpendels.

Aufgabe4:
Legen Sie einen Zustandsregler für das linearisierte System des Körperpendels so aus, dass 
wenn k der Betrag bei dem ursprünglich konjugiert komplexen Eigenwertpaar auf der imaginären Achse ist,

Code 0-5: dann die Wunscheigenwerte bei -k+k*i und -k-k*i liegen.

clear;

m=1;
l=1;
h=l/2;
r=1.1;
g=9.81;

J=0.25*m*r*r + (1/12)*m*l*l;

Q=(m*g*h)/(J+m*h*h);

function f = rechteSeite(t,y)
    phi = y(1,1);
    w = y(2,1);
    f(1,1) = w;
    f(2,1) = -Q*sin(phi);
endfunction

t = 0:0.01:15;
y0 = [170*%pi/180,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:)',t,y(2,:)');

Code 0-6: Lösung zu Aufgabe 1.

clear;

m=1;
l=1;
h=l/2;
r=1.1;
g=9.81;

J=0.25*m*r*r + (1/12)*m*l*l;

Q=(m*g*h)/(J+m*h*h);


//Zustandsregler auslegen:

A = [0,1;-Q,0]
EW = spec(A)
K  = abs(EW(1))
EWwunsch = [-K+%i*K,-K-%i*K]

B=[0;1/(J+m*h*h)]

R=ppol(A,B,EWwunsch)

//Simulation des geregelten Systems:

function f = rechteSeite(t,y)
    phi = y(1,1);
    w = y(2,1);
    MA = R(1)*phi + R(2)*w;
    f(1,1) = w;
    f(2,1) = -Q*sin(phi) - (1/(J+m*h*h))*MA;
endfunction

t = 0:0.01:15;
y0 = [170*%pi/180,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:)',t,y(2,:)');

Code 0-7: Lösung zu Aufgabe 4 nur in Scilab.

clear;

m = 1;
l = 1;
h = 1/2;
r = 0,1;
g = 9.81;
J = 1/4*m*r^2+1/12*m*l^2;

Q = (m*g*h)/J;

A = [0,1; Q,0];
B = [0; h/J];

EW = spec(A);
k = abs(EW(1));

P = [-k+k*%i;-k-k*%i];

R = ppol(A,B,P);

function f = rechts(t,y)
    phi = y(1,1);
    w = y(2,1);
    F_a = -R*y;
    f(1,1) = w
    f(2,1) = Q*phi+(h/J)*F_a;
endfunction

t = 0:0.01:2;
y0 = [0.1;0];
t0 = 0;

y = ode(y0,t0,t,rechts);

plot(t,y(1,:)',t,y(2,:)');

Code 0-8: Studentische Lösung, Teil 1. bis 5.

R = ppol(A,B,EW)

Code 0-9: Verwendung von ppol zur Auslegung eines Zustandsreglers unter Vorgabe der Wunscheigenwerte für A* durch EW.

clear;

m = 1;
l = 1;
h = 1/2;
r = 0.1;
g = 9.81;
J = 1/4*m*r^2+1/12*m*l^2;

Q = (m*g*h)/J;

A = [0,1,0,0; 0,0,0,0; 0,0,0,1; 0,0,Q,0];
B = [0; 1/m; 0; h/J];

EW = spec(A);
k = abs(EW(1));

P = [-k+k*%i; -k-k*%i;-2*k+k*%i; -2*k-k*%i];

R = ppol(A,B,P);

function f = rechts(t,y)
    x = y(1,1);
    v = y(2,1);
    phi = y(3,1);
    w = y(4,1);
    F_a = -R*y;
    sp = sin(phi);
    cp = cos(phi);
    f(1,1) = v;
    f(2,1) = (1/m)*F_a;
    f(3,1) = w;
//f(4,1) = Q*phi+(h/J)*F_a; //Lineare Variante 
    f(4,1) = (-w^2*m*h^2*sp*cp+m*g*h*sp+h*cp*F_a)/(J+m*h^2*sp^2);
endfunction

//Verwendung der Matrixform (nur für linearisierte Variante möglich):
//function dx=rechts(t, x)
//    u = -R*x;
//    dx = A*x + B*u;
//endfunction

t = 0:0.01:1;
y0 = [0.1; 0; 0.1; 0];
t0 = 0;

y = ode(y0,t0,t,rechts);

plot(t,y(1,:)',"r", t,y(2,:)',"b", t,y(3,:)',"g", t,y(4,:)',"m");

Code 0-10: Studentische Lösung(en)

################################

SNIPPET Zufallswerte in Scilab:

dx = floor(grand(1, 1, "unf", 0, 2.999999999)) - 1;  //liefert -1 oder +1 zufällig

################################

MODELL.SCE:

R1=10.153333;
R2=0.3433333;

m = 1;  //kg
g = 9.81; //m/s^2
h = 0.5; //m
r = 0.1; //m
l = 2*h; //m
J = 0.25*m*r^2 + (1/12)*m*l; // kg*m^2

function f = rechteSeite(t,y)
    phi = y(1,1);
    om  = y(2,1);
    f(1,1) = om;
    FA = -R1*phi -R2*om;
    N = J + m*h*h*sin(phi)*sin(phi);
    f(2,1) = (-m*h*h*om*om*sin(phi)*cos(phi) + h*m*sin(phi)*g + h*cos(phi)*FA )/N;
endfunction

################################

FEHLERFUNKTION.SCE

function z = berechneFehler(RR1,RR2)
    R1 = RR1;
    R2 = RR2;

    t = linspace(0,20,3000);
    y0 = [0.2,0.1]';
    t0 = 0;
    yy  = ode(y0,t0,t,rechteSeite);

    fehler = sum(abs(yy(1,:)));    
    z = fehler;    
endfunction

################################

OPTIMIERER.SCE

clear();

exec modell.sce
exec fehlerfunktion.sce

################################

SNIPPET ZUR ANIMATION

x = y(3,:);
phi = y(1,:);
y = h.*cos(phi);

zs = size(t);
spalten = zs(2);

pA = [x + h.*sin(phi);zeros(1,spalten)];
pB = [x - h.*sin(phi);y + h.*cos(phi)];

//Eine Animation soll eine Verbindungslinie von pA mit pB im zeitlichen Verlauf darstellen:
isoview;
clf;
// ABBRUCH MIT CTR C !!
for i=1:spalten
   clf;
   plot([-1.2,-1.1999,1.1999,1.2],[-1.2,1.2,1.2,-1.2]);  //Rahmen, um Skalierung konstant zu halten
   plot([pA(1,i),pB(1,i)],[pA(2,i),pB(2,i)]);
   sleep(10); // "Echtzeit"
end

Code 0-11: ...insbesondere die Programmstruktur in Scilab.

//s als Argument für ein Polynom definieren:
s = poly(0,"s"); 
//Übertragungsfunktion definieren:
G = syslin('c',[1],[1+2*s+s^2]);
t=[0:0.01:10]; 
u=ones(1,1001);
y=csim(u,t,G);
plot2d(t,y);
//plzr(G);

Code 0-12: Verwendung von Scilab zur Handhabung von Übertragungsfunktionen (Beispiel).

m = 1.0;
C = 1.0;
D = 1.0;
F = 10.0;
function f = rechteSeite(t,y)
    x = y(1,1);
    v = y(2,1);
    f(1,1) = v;
    f(2,1) = -(C/m)*x -(D/m)*v + F/m;
endfunction

t = linspace(0,30,3000);
y0 = [0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

plot(t,y(1,:)','gre');

// -----------------------------------------

// Laplace-Darstellung:


//s als Argument für ein Polynom definieren:
s = poly(0,"s"); 
//Übertragungsfunktion definieren:
G = syslin('c',[1],[m*s*s + D*s + C]);

u=10*ones(1,3000);  //entspricht F
y=csim(u,t,G);
plot(t,y,'blu--');
//plzr(G);

Code 0-13: Feder-Masse-Schwinger mit äußerer Kraft F im Zeit- und Laplacebereich simulieren.

Nachbesprechung Übung 2b letzte Woche
Verkettung zwei PT1 Laplace / Zeitbereich
Regeln Laplace
P-Regler Laplace / Zeitbereich
ÜBUNG: PD-Regler

Totzeitglied
Runge-Kutta-Integrator
Methode 1 nach Ziegler und Nichols 

Code 0-14: Themen

clear();

xtot=0;

function f=modell(y,t,dt)
x=y(1,1);
v=y(2,1);

Kkrit = 19.75;
xsoll = 1.0;
//e = xsoll - x;
e = xsoll - xtot;
u=Kkrit*e;

f(1,1)=v;
f(2,1)=-x-2*v+u;
endfunction

function yneu=ruku(y,t,dt)
k1=modell(y,t,dt);
k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
k4=modell(y+dt.*k3,t+dt);
yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

tmax = 20.0;
dt = 0.01;
schritte = ceil(tmax/dt);

yalt = [0,0]';
ysim = yalt;
t=0.0;
tt=t;

Ttot = 0.1;
anztot = round(Ttot/dt)
xtotarr = zeros(anztot,1);

for i=1:1:schritte
yneu=ruku(yalt,t,dt);
yalt=yneu;
ysim=[ysim,yalt];
tt =[tt,t];
t=t+dt;

xtot = xtotarr(modulo((i-1),anztot)+1);
xtotarr(modulo((i-1),anztot)+1)=yneu(1);
end

plot(tt,ysim(1,:));

Code 0-15: Beispielcode Methode 1 Ziegler / Nichols.

clear();

xtot=0;

function f=modell(y,t,dt)
x=y(1,1);
v=y(2,1);
eI=y(3,1); //Integral der Regeldifferenz

Kkrit = 19.75;
Tkrit = 1.45;

//PI-Regler:
K=0.45*Kkrit;
TN=0.85*Tkrit;

xsoll = 1.0;
//e = xsoll - x;
e = xsoll - xtot;
//u=Kkrit*e;
u=K*e + (K/TN)*eI;

f(1,1)=v;
f(2,1)=-x-2*v+u;
f(3,1)=e;
endfunction

function yneu=ruku(y,t,dt)
k1=modell(y,t,dt);
k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
k4=modell(y+dt.*k3,t+dt);
yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

tmax = 20.0;
dt = 0.01;
schritte = ceil(tmax/dt);

yalt = [0,0,0]'; //jetzt dritter Startwert eI0=0
ysim = yalt;
t=0.0;
tt=t;

Ttot = 0.1;
anztot = round(Ttot/dt)
xtotarr = zeros(anztot,1);

for i=1:1:schritte
yneu=ruku(yalt,t,dt);
yalt=yneu;
ysim=[ysim,yalt];
tt =[tt,t];
t=t+dt;

xtot = xtotarr(modulo((i-1),anztot)+1);
xtotarr(modulo((i-1),anztot)+1)=yneu(1);
end

plot(tt,ysim(1,:));

Code 0-16: Ergänzen des I-Anteils und Simulation mit eingebautem PI-Regler, der Mit Ziegler-Nichols ausgelegt wurde.

clear();

xtot=0;

function f=modell(y,t,dt)
x=y(1,1);
v=y(2,1);

Kkrit = 0.49;
Tkrit = 8.125;

C=0.1;
D=0.1;

//PI-Regler:
K=0.45*Kkrit;
I=0.85*Tkrit;

xsoll = 2.0;
//e = xsoll - x;
e = xsoll - xtot;
u=Kkrit*e;

f(1,1)=v;
f(2,1)=-C*x-D*v+u;
endfunction

function yneu=ruku(y,t,dt)
k1=modell(y,t,dt);
k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
k4=modell(y+dt.*k3,t+dt);
yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

tmax = 80.0;
dt = 0.01;
schritte = ceil(tmax/dt);

yalt = [0,0]';
ysim = yalt;
t=0.0;
tt=t;

Ttot = 0.2;
anztot = round(Ttot/dt)
xtotarr = zeros(anztot,1);

for i=1:1:schritte
yneu=ruku(yalt,t,dt);
yalt=yneu;
ysim=[ysim,yalt];
tt =[tt,t];
t=t+dt;

xtot = xtotarr(modulo((i-1),anztot)+1);
xtotarr(modulo((i-1),anztot)+1)=yneu(1);
end

plot(tt,ysim(1,:));

Code 0-17: Ziegler-Nichols, Kkrit finden.

clear();

xtot=0;

function f=modell(y,t,dt)
x=y(1,1);
v=y(2,1);
eI=y(3,1); //Integral der Regeldifferenz

Kkrit = 0.49;
Tkrit = 8.125;

C=0.1;
D=0.1;

//PI-Regler:
K=0.45*Kkrit;
TN=0.85*Tkrit;

xsoll = 2.0;
//e = xsoll - x;
e = xsoll - xtot;
//u=0.5*Kkrit*e;
u=K*e + 0.1*(K/TN)*eI;

f(1,1)=v;
f(2,1)=-C*x-D*v+u;
f(3,1)=e;
endfunction

function yneu=ruku(y,t,dt)
k1=modell(y,t,dt);
k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
k4=modell(y+dt.*k3,t+dt);
yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

tmax = 200.0;
dt = 0.01;
schritte = ceil(tmax/dt);

yalt = [0,0,0]'; //jetzt dritter Startwert eI0=0
ysim = yalt;
t=0.0;
tt=t;

Ttot = 0.2;
anztot = round(Ttot/dt)
xtotarr = zeros(anztot,1);

for i=1:1:schritte
yneu=ruku(yalt,t,dt);
yalt=yneu;
ysim=[ysim,yalt];
tt =[tt,t];
t=t+dt;

xtot = xtotarr(modulo((i-1),anztot)+1);
xtotarr(modulo((i-1),anztot)+1)=yneu(1);
end

plot(tt,ysim(1,:));

Code 0-18: Ziegler-Nichols, PI-Regler testen.

clear();

xtot=0;

function f=modell(y,t,dt)
x=y(1,1);
v=y(2,1);
eI=y(3,1); //Integral der Regeldifferenz

Kkrit = 0.49;
Tkrit = 8.125;

C=0.1;
D=0.1;

//PID-Regler:
K=0.6*Kkrit;
TN=0.5*Tkrit;
TV=0.12*Tkrit;

xsoll = 2.0;
//e = xsoll - x;
e = xsoll - xtot;
//u=0.5*Kkrit*e;
u=K*e + (K/TN)*eI - TV*v;

f(1,1)=v;
f(2,1)=-C*x-D*v+u;
f(3,1)=e;
endfunction

function yneu=ruku(y,t,dt)
k1=modell(y,t,dt);
k2=modell(y+0.5.*dt.*k1,t+0.5*dt);
k3=modell(y+0.5.*dt.*k2,t+0.5*dt);
k4=modell(y+dt.*k3,t+dt);
yneu=y+dt.*(k1+2.*k2+2.*k3+k4)./6;
endfunction

tmax = 200.0;
dt = 0.01;
schritte = ceil(tmax/dt);

yalt = [0,0,0]'; //jetzt dritter Startwert eI0=0
ysim = yalt;
t=0.0;
tt=t;

Ttot = 0.2;
anztot = round(Ttot/dt)
xtotarr = zeros(anztot,1);

for i=1:1:schritte
yneu=ruku(yalt,t,dt);
yalt=yneu;
ysim=[ysim,yalt];
tt =[tt,t];
t=t+dt;

xtot = xtotarr(modulo((i-1),anztot)+1);
xtotarr(modulo((i-1),anztot)+1)=yneu(1);
end

plot(tt,ysim(1,:));

Code 0-19: Ziegler-Nichols, PID-Regler testen.

//s als Argument für ein Polynom definieren:
s = poly(0,"s"); 
//Übertragungsfunktion definieren:
C=0.1;
D=0.1;

Kkrit = 0.49;
Tkrit = 8.125;

C=0.1;
D=0.1;

//PID-Regler:
K=0.6*Kkrit;
TN=0.5*Tkrit;
TV=0.12*Tkrit;

//K=0.45*Kkrit;
//TN=0.85*Tkrit;


G = syslin('c',[1],[s*s+C*s+D]);
R = syslin('c',[K*TN*TV*s*s+K*TN*s+K],[TN*s]);
H=G*R;
Q=H/(H+1);
t = linspace(0,30,3000);
u=2*ones(1,3000);  //entspricht F
y=csim(u,t,Q);
//plot(t,y,'blu--');
plzr(Q);

Code 0-20: Prüfen der Eigenwertlagen.

s  =   poly(0,"s");
Q  =   ... (Ihre Lösung)              

Code 0-21: Scilab-Code.

r1  =  ...
r2  =  ...

Code 0-22: Angabe der Regelmatrix, d.h. der Komponenten r1 und r2.

Das System ist ...

Code 0-23: Antwort.

s  =   poly(0,"s");
G  =  ... (Ihre Lösung)
t = linspace(0.0,1.0,10);
anz = size(t);
u = ones(anz(1),anz(2));
ystrecke   = csim(u,t,G);
disp(G);

Code 0-24: Scilab-Code.

s  =   poly(0,"s");
R  =  .... (Ihre Lösung)
disp(R);

Code 0-25: Scilab-Code