// Studio dell'equazione del pendolo:
//
// (*)    x" = - (g/L) sen(x)
//
clear;
//
// Costanti:
//
g = 9.81; // m/s^2
L = 1; // m
p = sqrt(g/L); // 1/s
//
// 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) ;
         - p^2 * sin(x(1)) ];
endfunction
//
// Condizioni iniziali e istante finale;
//
t0 = 0;
u0(1) = %pi/8;
u0(2) = 0;
tf = 3;
//
// Metodo di RK12
//
exec('/media/TDK NERA/Lez Veicoli/RK12_pv.sci');
//
// Valore max stima errore locale
//
E = 1d-6;
//
[T,U,PASSO] = RK12_pv(u0,t0,tf,F,E);
//
// Soluzione caso linearizzato
//
function y = SE(t)
    y = [     u0(1) * cos(p * t);
         -p * u0(1) * sin(p * t)];
endfunction
//
// Grafici
//
clf();
subplot(221);
XE = SE(T);
plot2d(T',[U(1,:)',XE(1,:)'],style = [5,2]);
xgrid(); legend('U(1,:)','y(t)','in_lower_left');
subplot(223);
plot2d(T',[U(2,:)',XE(2,:)'],style = [5,2]);
xgrid(); legend('U(2,:)','y''(t)'); xlabel('t');
subplot(222);
plot2d(T(1:$-1)',PASSO',style = 5);
xlabel('t');ylabel('PASSO');
xgrid();
subplot(224);
plot2d(U(1,:),U(2,:),style = 5);
plot2d(XE(1,:),XE(2,:),style = 2);
xgrid();
xlabel('x');ylabel('x''');
legend('NL','LIN');
title('Traiettoria nel piano delle fasi')
