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