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];