kill(all);

assume(omega>0);
assume(zeta>0);
assume(r>0);

m*acc + c*vel + k*disp = f

con

acc  = diff(u(t),t,2),
vel  = diff(u(t),t,1)
disp =      u(t)

normalizzo dividendo per m e chiamo

c/m = 2*zeta*omega
k/m = omega^2
q(t)= f(t)/m


scrivo il residuo dell'equaz. di equilibrio dinamico

res: 1 * diff(u(t),t,2) + 2*zeta*omega * diff(u(t),t,1) + omega^2 *u(t) -q(t);

ho forzante armonica modulata in pulsazione r*omega, con r arbitrario

q(t) = qr*cos(r*omega*t)-qi*sin(r*omega*t)

q(t):= qr*cos(r*omega*t)-qi*sin(r*omega*t);

modulo forzante

qm:sqrt(qr^2+qi^2);

suppongo la risposta armonica, e modulata sulla stessa pulsazione della forzante

u(t):= ur*cos(r*omega*t)-ui*sin(r*omega*t);

u(t):= ur*cos(r*omega*t)-ui*sin(r*omega*t);

valuto il residuo sotto queste nuove ipotesi, valutando le derivate

res:ev(res,diff)$
res:expand(res);

estraggo i coefficienti dei fattori coseno

res_cos:coeff(res,cos(omega*r*t));

estraggo i coefficienti dei fattori seno

res_sin:coeff(res,sin(omega*r*t));

controllo non vi siano residui

fullratsimp( res - res_cos*cos(r*omega*t) - res_sin*sin(r*omega*t));

annullo i coeff. dei fattori seno e coseno trovando la soluzione in termini di ur e ui 

linsolve([res_cos=0,res_sin=0],[ur,ui]),globalsolve=true;

modulo risposta

um:fullratsimp(sqrt(ur^2+ui^2));

amplificazione rispetto a risposta statica

um/(qm/omega^2);
factor(%);
define(
    amp(r,zeta) , 
    %
);

amp(1,zeta);

grafico amplificazione dinamica

wxplot2d(
[
amp(r,10^(-2  )),
amp(r,10^(-1.5)),
amp(r,10^(-1  )),
amp(r,10^(-0.5)),
amp(r,10^( 0  )),
amp(r,10^( 0.5))
],
[r,0,2],
[logy],
[legend,"0.01","0.0316","0.1","0.316","1","3.16"]

);

grafico fase

atan2(ui*denom(ur),ur*denom(ur))-atan2(qi,qr),qr=1,qi=0;
define(
    fase(r,zeta),
    %
);

wxplot2d(
[
fase(r,10^(-2  )),
fase(r,10^(-1.5)),
fase(r,10^(-1  )),
fase(r,10^(-0.5)),
fase(r,10^( 0  )),
fase(r,10^( 0.5))
],
[r,0,2],
[y,-%pi,0],
[legend,"0.01","0.0316","0.1","0.316","1","3.16"]
);