ex8.m

% Lösung von Aufgabe 8

% Massen- und Steifigkeitsmatrix
M = diag([6 6 1])*10^3;               % in kg
C = [6 -3 0;-3 4 -1; 0 -1 1]*10^6;    % in N/m

% charakteristische Gleichung mit Symbolic Toolbox
syms om2;
A = -om2*M + C;
eq = det(A);
charEq = sym2poly(eq)';
charEq = charEq/charEq(1)    % zum Vergleichen auf 1 normiert

% Nullstellen
om2 = roots(charEq)

% Daraus die Eigenfrequenzen, aufsteigend sortiert
om = sort(sqrt(om2))

% 1. Eigenvektor bestimmen
om1 = om(1);
A = -om1^2*M + C
% Lösung der homogenen Gleichung (Länge = 1)
x1 = null(A)
% der Vergleichbarkeit wegen: 1. Element = 1
x1 = x1/x1(1)

% analog 2. Eigenvektor:
format long
A = -om(2)^2*M + C
format short
x2 = null(A)
% klappt nicht wegen numerischer Ungenauigkeit!
% hier hilft eine Option von null
x2 = null(A, 'r')
x2 = x2/x2(1)

% 3. Eigenvektor ginge hier mit null, alternativ direkt durch Lösen des
% Gleichungssystems
A = -om(3)^2*M + C;
xRest = A(1:end-1, 2:end) \ (-A(1:end-1,1))
x3 = [1; xRest]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Anfangsbedingungen anpassen
x0 = [0, 0, 0]';
v0 = [0, 0, 10]';

% Modalmatrix
Phi = [x1, x2, x3]
% Hilfsmatrix m
m = Phi'*M*Phi
% Inverse von Phi
% PhiInv = inv(m)*Phi'*M   % ok, es geht aber auch ganz ohne inv
PhiInv = diag(1./diag(m))*Phi'*M

% Anfangsbedingungen in Hauptkoordinaten
yv0 = PhiInv*v0
yhat = abs(yv0./om)
phi = sign(yv0)*pi/2

% Vorfaktoren der x-Terme
VF = Phi*diag(yhat.*sign(phi))