Pendolo con MATLAB

motoArmonico.m

function [t, theta] = motoArmonico(tMax, theta0)
    [t, vars] = ode45(@armonico, [0 tMax], [0 theta0]);
    theta = vars(:,2); % la prima colonna è la derivata prima
end

function dthdt = armonico(~,th) % questo è per il moto armonico semplice
    dthdt(1) = th(2);
    dthdt(2) = -th(1);
    dthdt = dthdt(:);
end

motoPendolo.m

function [t, theta] = motoPendolo(tMax, theta0)
    [t, vars] = ode45(@armonico, [0 tMax], [0 theta0]);
    theta = vars(:,2); % la prima colonna è la derivata prima
end

function dthdt = armonico(~,th)
    % per il moto di un pendolo fuori dal regime di piccole oscillazioni
    dthdt(1) = th(2);
    dthdt(2) = -sin(th(1));
    dthdt = dthdt(:);
end

pendolo.m

% plotta i valori di oscillazione di un pendolo e di un moto armonico
% per vedere di quanto varia l'errore dell'approssimazione
t0 = 0;
tf = 20;
thetas = [pi/8 pi/4; 3*pi/8 pi/2]; % valori a cui fare i diversi plot 
[m, n] = size(thetas);
clf;
layout = tiledlayout(m, n);
xlabel(layout, 't');
ylabel(layout, 'theta');
title(layout, 'Moto di un pendolo e moto armonico (l=g=1)');
layout.TileSpacing = 'Compact';
layout.Padding = 'Compact';
for i = 1:m
    for j = 1:n
        nexttile
        [tPendolo, thetaPendolo] = motoPendolo(tf, thetas(i,j));
        [tArmonico, thetaArmonico] = motoArmonico(tf, thetas(i,j));
        plot(tPendolo, thetaPendolo, '-', ...
            tArmonico, thetaArmonico, '-');
        legend('pendolo', 'armonico');
        title(['Moto per thetaMax = ' num2str(thetas(i,j))]);
    end
end

periodoPendolo.m

% mostra il periodo del pendolo in funzione dell'angolo massimo
global l;
l = 1;
global g;
g = 1; % siamo su un pianeta particolare dove g = 1
fplot(@periodo, [0 pi/2]);
title('Periodo in funzione di thetaMax (l=1, g=1)');
function t = periodo(thetaMax)
    global l;
    global g;
    t = 2*pi*sqrt(l/g)*(1+(thetaMax^2)/16);
end