Tema - Pendulul Anarmonic

clear; close all; clc;

g = 9.80665; % acceleratia gravitationala (m/s^2)
lungime = 0.5; % lungimea firului (m)
omega0 = sqrt(g/lungime); % pulsatia oscilatiilor armonice (amplitudini mici)
T0 = 2*pi/omega0; % perioada proprie oscilatiilor armonice
ti = 0; tf = 15*T0; N = 700;
t = linspace(ti, tf, N); % sirul momentelor de timp
dt = t(2) - t(1); % pas de timp
theta_i = 45*pi/180; % unghiul initial (intre Oy si fir)

theta_an = theta_i * sin(omega0 * t + pi/2); % solutia armonica

figure(1); % reprezentarea legii unghiulare de miscare
plot(t, theta_an*180/pi, '-r'); hold on;
xlabel('t/s'); ylabel('theta/grade'); grid;
title('Legea unghiulara de miscare');

theta_num = zeros(1, N); % prealocare solutie numerica
theta_num(1) = theta_i; theta_num(2) = theta_i; % conditii initiale ale recurentei

aux1 = g/lungime*dt^2;
aux2 = dt/g;
for i=2 : 1 : N-1
    % Relatia de recurenta - oscilatie amortizata
    theta_num(i+1) = 2*theta_num(i) - theta_num(i-1) - aux1*sin(theta_num(i)) - aux2*(theta_num(i) - theta_num(i-1));
end

plot(t, theta_num*180/pi, '-', 'Color', [0.1 0.5 0.1]);
legend('Analitic', 'Numeric');

figure(2); % simularea dinamica a oscilatiei
xp = lungime*sin(theta_num);
yp = -lungime*cos(theta_num);

for i=1 : 1 : N
    plot([0 xp(i)], [0 yp(i)], '-k', xp(i), yp(i), 'ob');
    xlabel('x/m'); ylabel('y/m'); grid;
    title('Simularea Oscilatiei');
    axis equal;
    axis([-lungime lungime -lungime 0]); % limitele graficului
    Film(i) = getframe();
end

figure(3); % simularea in timp real
xlabel('x/m'); ylabel('y/m'); grid;
title('Simularea oscilatiei in timp real');
axis equal;
axis([-lungime lungime -lungime 0]); % stabileste limitele graficului [xmin xmax ymin ymax]
movie(Film, 1, N/tf);