kill(all);

Si definiscono positive le seguenti quantità: 
-omega = frequenza propria;
-zeta = smorzamento frazione critico;
- r = rapporto tra la pulsazione o la frequenza dell'eccitante e quella della frequenza propria.

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

Il sistema ad un grado di libertà può essere immaginato come una massa il cui unico gld è la traslazione verticale (u(t)).
Tale massa risulta ancorata a terra da una molla  di rigidezza K e da un elemento smorzante la cui caratteristica smorzante
 è C e su di essa è applicato un carico esterno F(t) definito come una smorzante armonica. 

Si effettua ora l'equilibrio del sistema in figura dopo aver prima definito l'accelerazione 
e la velocità della massa m come derivate del vettore u(t). 
L'equazione dinamica del sistema ad un grado di libertà è presentata nella seguente forma:

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

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

normalizzando dividendo per la massa verranno definiti i seguenti rapporti:
 
c/m    = 2*zeta*omega
k/m    = omega^2
f(t)/m = q(t)
u(t)   = xi(t)

Pertanto la nuova espressione sarà: 

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

risulta:
- omega  essere la pulsazione propria
- zeta=1 essere lo smorzamento di transizione tra 
         una amplificazione dinamica >=1 e una <=1 al variare di f

Si pone ora in forma di  residuo l'equaz. di equilibrio dinamico (il fattore 1 è  la massa resa unitaria dalla normalizzazione)

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

Si inserisce una forma armonica della forzante  modulata in pulsazione r*omega, 
con r arbitrario,in modo da non rinominare il parametro omega precedentemente definito.

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

in cui:
- qr = fr/m 
- qi = fi/m
in questo modo, modulando qr, qi posso dare arbitrario
modulo e arbitraria base.

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

modulo forzante, costante nel tempo

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

Si suppone la risposta armonica e modulata sulla stessa pulsazione della forzante. Ciò può essere
spiegato nel seguente modo:  nel lungo  termine ci si aspetta che il sistema che è soggetto all'azione pulsata r*omega 
si adegui a tale pulsazione e si metta a vibrare con tale frequenza. Pertanto, tra tutte le soluzioni possibili, si considera
una soluzione a lungo termine la cui forma specifica è un'oscillazione modulata come di seguito.

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

da ur e ui è possibile calcolare il modulo

xi(t):= xir*cos(r*omega*t)-xii*sin(r*omega*t);

Si valuta ora il residuo sotto queste nuove ipotesi, valutando le derivate 

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

Le derivate rimangono in forma simbolica perche non sono 
note la forma specifica di u(t);

Il residuo calcolato deve risultare nullo per ogni instante temporale, pertanto si
pone sotto la forma:  
A * cos(r*omega*t) + B *sin(r*omega*t)=0
ottenuto  raccogliendo tutti i termini coseno moltiplicati 
A e tutti i termini seno moltiplicati  per B, imponendo
A=0,B=0 la relazione sopra scritta sarà verificata per ogni t;    

Tramite il comando 'coeff' si estraggono i coefficienti dei fattori coseno del residuo

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

e i coefficienti dei fattori seno del residuo

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

Si controlla non vi siano residui di residuo tramite la funzione fullratsimp sottraendo all'espressione 
del residuo i coefficienti in coseno e in seno, il risultato dovrebbe essere zero

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

Si scrive ora il sistema di due equazioni annullando i coeff. dei fattori seno e coseno, 
si trova in questo modo la soluzione in termini di ur e ui e si associa rispettivamente ai
residui res_cos e res_sin le variabili xir e xii 

linsolve(
[ res_cos=0, res_sin=0 ],[xir,xii]),globalsolve=true;

Si trova il modulo risposta, funzione dell'eccitante, di omega, di r e di zeta;

omega = sqrt(k/m) è l'unica pulsazione propria del sistema,  

xim:fullratsimp(sqrt(xir^2+xii^2));

Prima di tutto si va a normalizzare la risposta appena trovata
rispetto ad una risposta statica; quest'ultima si ottiene quando si prende una massa 
unitaria supportata da una rigidezza (k/m), applico una forza esterna che in modulo è sqrt(qr^2+qi^2)
e lo spostamento statico sarà questa entità diviso omega quadro;

Si definisce quindi il fattore che lega il cedimento dinamico e il cedimento statico
come rapporto tra xim (cedimento din) e qm (cedimento stat)
La funzione factor permette di eliminare omega^2

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

Si definisce un fattore di amplificazione dinamico funzione del rapporto tra la pulsazione dell'eccitante 
e la pulsazione propria dell'oscillatore e dello smorzamento frazione critico zeta;
si valuta ora cosa succede alla risonanza, cioè quando il rapporto tra eccitante e pulsazione propria è nullo,
e si trova che l'amplificazione è 1/2*zeta, e se zeta tende a 0 la risposta tende a infinito

amp(1,zeta);

Si plotta ora il grafico dell'amplificazione dinamica in cui si hanno due parametri:
il rapporto r e smorzamento frazione critico zeta.

Al variare dello smorzamento frazione critico si ottengono cinque curve, in cui la curva nera è quella 
calcolata con valore critico 1
 
Sull'asse y i punti sono equidistanziati in scala logaritmica.

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"]
);

Dal grafico si evince che si può andare oltre il valore critico di smorzamento 
(curva celeste) e non si riscontra nessuna singolarità.

Per valori della frequenza dell'eccitante molto più piccoli di quelli della frequenza propria
la soluzione dinamica converge sulla statica; la soluzione dinamica inizia 
a divergere dalla statica man mano ci si avvicina all'ultima frequenza propria.
Qui il picco è limitato solo dalla presenza di smorzamento, per piccoli smorzamenti 
il picco è più o meno centrato sulla condizione risonanza. 

Si nota che quando lo smorzamento critico è:
 <1 la risposta dinamica, almeno per una frequenza, è leggermente amplificata  rispetto alla statica;
 =1 la risposta dinamica ()rappresentata dalla curva nera) non è mai maggiore di 1;
>1 la risposta crolla quasi subito sotto il livello 1 per cui la risposta dinamica è sotto la statica.

grafico fase; la moltiplicazione di ambo gli argomenti di atan2 
per lo stesso fattore non varia il risultato (ma semplifica la formula), alla fase della
risposta sottraggo la fase dell'eccitante;

Si plotta ora, in funzione di r e zeta, il grafico della fase ottenuta come differenza tra la fase della risposta, che è definita come l'arcotangente2 della parte seno  e della parte coseno, e la fase dell'eccitante definita in modo analogo, specificando un particolare eccitante con sola parte coseno (ponendo a zero la parte seno).

atan2(
xii*denom(xir),
xir*denom(xir)
)
-
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"]
);

Dal grafico si evince che la transizione fra fase e controfase diventa piu' smussata man mano che aumenta
lo smorzamento.

Si prosegue aprendo il foglio di calcolo: risposta_per_sovrapposizione_modale.wxmx