solveVibrationODE.m
function [t, d, v] = solveVibrationODE(M, C, d0, v0, tEnd, dt)
% löst die (Nx2)- dimensionale Schwingungs-DGL für das Zeitintervall
% [0 tEnd]
% Parameter
% M Massenmatrix (2N x 2N)
% C Steifigkeitsmatrix (2N x 2N)
% d0 Anfangsauslenkungen (2 x N)
% v0 Anfangsgeschwindigkeiten (2 x N)
% tEnd Endzeit
% dt optional: Ausgabe-Schrittweite
% Ergebniswerte:
% t Zeitwerte (NT-Vektor, vom Solver bestimmt)
% d Auslenkungen (NT x 2 x N)
% v Geschwindigkeiten (NT x 2 x N)
% mache aus d0 und v0 einen langen Vektor
N = size(d0,2);
y0 = reshape([d0, v0], 4*N, 1);
% löse DGL
odeFunc = @(t,y) vibrationODE(t, y, M, C);
if (nargin < 6)
[t, y] = ode45(odeFunc, [0 tEnd], y0);
else
[t, y] = ode45(odeFunc, 0:dt:tEnd, y0);
end
% pflücke den Ergebnisvektor auseinander
NT = length(t);
d = y(:, 1:(2*N));
v = y(:, (2*N+1):end);
% und ordne ihn um
d = reshape(d, NT, 2, N);
v = reshape(v, NT, 2, N);
%----------------------------------------------
function dy = vibrationODE(t, y, M, C)
% Zustandsform der Schwingungsgleichung
% M x'' + C x = 0
% Dabei ist
% y = [x; v]
N = length(y)/2;
x = y(1:N);
v = y(N+1:end);
dy = [v; -inv(M)*C*x];