// Studio dell'equazione del pendolo:
//
// (*)    x" = - (g/L) sen(x)
//
clear;
clf();
//
// Costanti:
//
g = 9.82; // 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,u)
    // u(1) = x
    // u(2) = x'
    y = [             u(2) ;
         - p^2 * sin(u(1)) ];
endfunction
//
// Funzione G2 che restituisce la derivata seconda...
//
function y = G2(t,u)
    y = [       - p^2 *  sin(u(1)) ;
         - p^2 * u(2) .* cos(u(1)) ];
endfunction
//
// Condizioni iniziali e istante finale;
//
t0 = 0;
u0(1) = %pi/8;
u0(2) = 0;
tf = 3;
//
// Metodo di Eulero esplicito
//
exec('/home/ciampa/Lavoro/PaginaPersonale/public_html_pagine_dm/CalcNum/Lezioni/15-16/LMV/Lez/Lez Veicoli/TS_1_pv.sci');
//
// Valore max stima errore locale
//
E = 1d-5;
//
[T,U,PASSO] = TS_1_pv(u0,t0,tf,F,G2,E);
//
// Soluzione caso linearizzato
//
om = sqrt(g/L);
//
function y = SE(t)
    y = [      u0(1) * cos(om * t);
         -om * u0(1) * sin(om * t)];
endfunction
//
// Grafici
//
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')
