ex22.m

function ex22()
% Lösung von Aufgabe 22
scale = [1 10 100 1e3 1e4 1e5 1e6];
solver = {@ode45, @ode113, @ode15s};
solvername = {'ode45', 'ode113', 'ode15s'};

% a. Lösen und Zeiten nehmen
fprintf('   Skala   Solver    Zeit    nSteps\n');
for J = 1:length(solver)
  for I = 1:length(scale);
    s = scale(I);
    [time nSteps] = getTimeSteps(s, solver{J});
    fprintf('%6.2e  %7s  %6.2f  %8d\n', ...
        s, solvername{J}, time, nSteps);
  end
end

% b. hin und zurück lösen und Unterschied ausgeben
fprintf('\n\n   Skala   Solver    Fehler\n');
for J = 1:length(solver)
  for I = 1:length(scale);
    s = scale(I);
    d = getDistance(s, solver{J});
    fprintf('%6.2e  %7s  %8.2e\n', s, solvername{J}, d);
  end
end

% zur Kontrolle Ergebnisse plotten
for I = 1:6
  s = scale(I+1);
  [d, t1, y1, t2, y2] = getDistance(s, @ode15s);
  subplot(2,3,I)
  plot(t1, y1(:,1), t2, y2(:,1))
  titleString = sprintf('s=%.0e', s);
  title(titleString)
end

F = getframe(gcf);
imwrite(F.cdata, 'bild73.png');

% Versuch mit höherer Genauigkeit
[d, t1, y1, t2, y2] = getDistance(1e2, @ode15s, 1e-5);
figure()
plot(t1, y1(:,1), t2, y2(:,1))


%-----------------------------------------------------------------
function [time nSteps] = getTimeSteps(s, solver)
% löst die Schwingungs-DGL und gibt alle Messwerte zurück
% verändert b_2 und c_2 um den Faktor s
% verwendet den als String angegebenen Solver

tic;
[t, y] = solver(@(t,y) f2d(t,y,s), [0 2], [0 0 0 0]);
time = toc;
nSteps = length(t);

%--------------------------------------------------------------------
function [d, t1, y1, t2, y2] = getDistance(s, solver, tol)
% löst die Schwingungs-DGL hin und zurück und bestimmt den Unterschied
% gibt auch beide Lösungen zurück zum Plotten

% optionales 3. Argument für den relativen Fehler
if nargin == 2
  options = odeset();              % nimmt Standardwerte
else
  options = odeset('RelTol', tol);
end
    
[t1, y1] = solver(@(t,y) f2d(t,y,s), [0 2], [0 0 0 0], options);
[t2, y2] = solver(@(t,y) f2d(t,y,s), [2 0], y1(end,:), options);

% Unterschied von x1 = y1: am Anfang ist x1 = 0, also
d = abs(y2(end,1));


%---------------------------------------------------------------------
function dydt = f2d(t, y, s)
% Schwingungs-DGL 2d, Skalenfaktor s für c2/b2 als Parameter

m1 = 5; m2 = 0.5;
c1 = 10; c2 = s*1;
b1 = 0.1; b2 = s*0.1;
F1 = 1;
Om = 1.414;
M = [m1, 0; 0, m2];
C = [c1+c2,-c2;-c2,c2];
B = [b1+b2,-b2;-b2,b2];
Fhat = [F1; 0];
dydt = [y(3:4); inv(M)*(Fhat*cos(Om*t) - B*y(3:4) - C*y(1:2))];