ex13.m

function ex13()
% Loesung von Aufgabe 13

L = 2;              % m
alpha = 25*pi/180;  % rad
g = 9.81;           % m/s^2

% q0 und T bei kleiner Auslenkung
q0 = L*tan(alpha)
omega0 = sqrt(g/L)*sqrt(cos(alpha)^3 / (1 + sin(alpha)^2))
T = 2*pi/omega0

dx = q0 - [0.1 0.5 0.93]; % Anfangswerte

func = @(t,y) dgl13(t, y, L, alpha);        % Parameter in dgl13 einsetzen
[t1, y1] = ode45(func, [0 6], [dx(1) 0]);
[t2, y2] = ode45(func, [0 6], [dx(2) 0]);
[t3, y3] = ode45(func, [0 6], [dx(3) 0]);

plot(t1, y1(:,1), t2, y2(:,1), t3, y3(:,1));
hold on;
plot([0 6], [q0 q0], 'k-');   % Gleichgewichtslage 
plot([T T], [0 2.1], 'k:');   % Schwingungsdauer, linearisiert 
hold off;

% Beschriftungen
title('Schwingungen bei verschiedenen Auslenkungen', 'FontSize', 16);
ylabel('q [m]', 'FontSize', 14);
xlabel('t [s]', 'FontSize', 14);
axis([0 6 0 2.1]);
legend('0.1 m', '0.5 m', '0.93 m', 0)

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

%----------------------------------------------------------------------
function dydt = dgl13(t, y, L, alpha)
% rechte Seite der DGL bei Aufgabe 13

dy1 = y(2);
g = 9.81;               % einfach konstant
fak = 1./sqrt(L^2 + y(1).^2);
dy2 = (g*sin(alpha) - L^2*y(1).*fak.^4.*y(2).^2 - g*fak.*y(1)) ...
     ./ (1 + y(1).^2.*fak.^2);

dydt = [dy1 dy2]';