// Studio dell'equazione del pendolo:
//
// (*)    x" = - (g/L) sen(x)
//
// (CI)   x(t0) = x0   ,   x'(t0) = v0
//
clear;
f1 = scf(1)
clf(1);
//
// Costanti:
//
g = 9.82; // m/s^2
L = 1; // m
m = 1; // kg
//
// Funzione che definisce l'equazione differenziale (il sistema
// di due equazioni di ordine uno equivalente all'equazione (*))
//
function y = F(t,x)
    // x(1) = x
    // x(2) = x'
    y = [              x(2) ;
         -(g/L) * sin(x(1)) ];
endfunction
//
// Condizioni iniziali e istante finale;
//
t0 = 0;
x0 = %pi/4; // rad
v0 = 0; // rad/s
u0(1) = x0;
u0(2) = v0;
tf = 3; // s
//
// Metodo RK12 (Eulero esplicito + Heun)
//
exec('/home/maurizio/Documenti/scilab-2025.1.0/Lavoro/ODE/RK_12_pv.sci');
//
// Valore max stima errore locale
//
E = 1d-3;
//
//[T,U,PASSO] = RK12_pv(u0,t0,tf,F,E,TAU = 5d-1);
[T,U,PASSO] = RK12_pv(u0,t0,tf,F,E);
//
// Energia meccanica
//
function y = Energia(x)
    y = m*g*L*(1 - cos(x(1))) + (1/2)*m*L^2 * x(2)^2;
endfunction
//
EN = zeros(T);
for k = 1:length(T),
    EN(k) = Energia(U(:,k));
end
//
subplot(411)
plot(T,EN);
xgrid();
ylabel('energia');
title('RK12')
//
subplot(412)
plot(T,U(1,:));
xgrid();
ylabel('angolo');
//
subplot(413)
plot(T,U(2,:));
xgrid();
ylabel('velocità angolare');
//
subplot(414)
plot(T(1:$-2),PASSO(1:$-1));
xgrid(); xlabel('t'); ylabel('PASSO');
//
printf('\nMetodo RK12   E = %3.2e  (numero istanti integrazione = %d)\n',E,length(T));
printf('\nVariazione percentuale dell''energia: %3.2f\n',100*(max(EN) - min(EN))/EN(1));
printf('\nVariazione percentuale ampiezza: %3.2e\n',100*(max(U(1,:)) + min(U(1,:)))/ u0(1));
