// Studio dell'equazione del pendolo:
//
// (*)    x" = - (g/L) sen(x)
//
// (CI)   x(t0) = x0   ,   x'(t0) = v0
//
clear;
f0 = scf(0)
clf(0);
//
// 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
//
// Funzione G2 che restituisce la derivata seconda...
//
function y = G2(t,x)
    y = [        -(g/L) *  sin(x(1)) ;
         - (g/L) * x(2) .* cos(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 di Eulero esplicito (TS(1))
//
exec('/home/maurizio/Documenti/scilab-2025.1.0/Lavoro/ODE/TS_1_pv.sci');
//
// Valore max stima errore locale
//
E = 1d-3;
//
[T,U,PASSO] = TS_1_pv(u0,t0,tf,F,G2,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(311)
plot(T,EN);
xgrid();
ylabel('energia');
title('TS(1)')
//
subplot(312)
plot(T,U(1,:));
xgrid();
ylabel('angolo');
//
subplot(313)
plot(T(1:$-2),PASSO(1:$-1));
xgrid(); xlabel('t'); ylabel('PASSO');
//
printf('\nE = %3.2e\n',E);
printf('\nVariazione percentuale dell''energia: %3.2f\n',100*(EN($) - EN(1))/EN(1));
printf('\nVariazione percentuale ampiezza massima: %3.2f\n',100*(max(U(1,:)) - u0(1))/u0(1));
printf('\nVariazione ampiezza massima: %3.2e\n',max(U(1,:)) - u0(1));
//halt('...premere un tasto per continuare...')
pause;
//
//
//
// Funzione G3 che restituisce la derivata terza...
//
function y = G3(t,x)
    y = [        -(g/L) * x(2) * cos(x(1)) ;
          (g/L)^2 * sin(x(1)) + (g/L) * x(2).^2 .* sin(x(1)) ];
endfunction
//
//
// ***************** Metodo TS(2)
//
exec('/home/maurizio/Documenti/scilab-2025.1.0/Lavoro/ODE/TS_2_pv.sci');
//
//
[T2,U2,PASSO2] = TS_2_pv(u0,t0,tf,F,G2,G3,E);
//
//
EN = zeros(T2);
for k = 1:length(T2),
    EN2(k) = Energia(U2(:,k));
end
//
f1 = scf(1);
clf(1);
subplot(311)
plot(T2,EN2);
xgrid();
ylabel('energia');
title('TS(2)')
//
subplot(312)
plot(T2,U2(1,:));
xgrid();
ylabel('angolo');
//
subplot(313)
plot(T2(1:$-2),PASSO2(1:$-1));
xgrid(); xlabel('t'); ylabel('PASSO');
//
printf('\nCon TS(2)  e  E = %3.2e\n',E);
printf('\nVariazione percentuale dell''energia: %3.2e\n',100*(max(EN2) - min(EN2))/EN2(1));
printf('\nVariazione percentuale ampiezza massima: %3.2e\n',100*abs(max(U2(1,ceil(size(U2,'c')/2):$)) - u0(1))/u0(1));
printf('\nVariazione ampiezza massima: %3.2e\n',abs(max(U2(1,ceil(size(U2,'c')/2):$)) - u0(1)));
