Дифференциальное неоднородное уравнение второго порядка с постоянными коэффициентами

clear 
clc 
close all
 
syms s F F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 F13 F14 EI h mpi
syms F15 F16 F17 F0011 F0012 F0013 F0014 F0015 F0016 F0017 
syms F21 F22 F23 F24 F25 F26 F27    A a q v tt t
 
%%  1   DATI DI INPUT
 
mpi=1;              %massa del piano
h=3;                %altezza interpiano
EI=10;              %rigidezza flessionale
Q0=15;              %spostameno modale iniziale
V0=5;               %velocita' modale iniziale
 
ag=0.8;             %accelerazione al suolo
omega=5;            %frequenza del terremoto
 
H=[h;2*h;3*h;4*h;5*h;6*h;7*h];
 
%%  2   MATRICE DI MASSA E DI RIGIDEZZA'
 
kpi=24*EI/h^3;      %rigidezza flessionale del singolo piano
m=mpi*eye(7);      %matrice delle masse
k=([[2*kpi,-kpi,0,0,0,0,0];[-kpi,2*kpi,-kpi,0,0,0,0];[0,-kpi,2*kpi,-kpi,0,0,0];[0,0,-kpi,2*kpi,-kpi,0,0];[0,0,0,-kpi,2*kpi,-kpi,0];[0,0,0,0,-kpi,2*kpi,-kpi];[0,0,0,0,0,-kpi,kpi]]);
%kp     -   matrice di rigidezza del sistema
 
%%  SISTEMA NON SMORZATO - OSCILLAZIONI LIBERE
%%  3   POLINOMIO CARATTERISTICO + AUTOVALORI + AUTOVETTORI:    sistema non smorzato
 
Pn=k-s^2*m;    %polinomio caratteristico senza smorzamento per q=b1*con OMEGAn*t + b2*sin OMEGAn*t
 
ssn=eig(Pn);      %autovalori del polinomio caratteristico
 
sn1=solve(ssn(1),s);
sn1=double(sn1);
 
sn2=solve(ssn(2),s);
sn2=double(sn2);
 
sn3=solve(ssn(3),s);
sn3=double(sn3);
 
sn4=solve(ssn(4),s);
sn4=double(sn4);
 
sn5=solve(ssn(5),s);
sn5=double(sn5);
 
sn6=solve(ssn(6),s);
sn6=double(sn6);
 
sn7=solve(ssn(7),s);
sn7=double(sn7);
 
Sn=[sn1(1);sn2(1);sn3(1);sn4(1);sn5(1);sn6(1);sn7(1)];   %tutti gli autovalori 
Sn=unique(Sn);
% KSI=cell(1,14);
n1=length(Sn);
% F=F*ones(1,n1);
% F=[F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14];   %per caso smorzato
F=[F1;F2;F3;F4;F5;F6;F7];           %per caso non smorzato
 
% P1=[P(1,:);P(1,:);P(2,:);P(2,:);P(3,:);P(3,:);P(4,:);P(4,:);P(5,:);P(5,:);P(6,:);P(6,:);P(7,:);P(7,:)];
% %caso smorzato
PPn=[Pn(1,:);Pn(2,:);Pn(3,:);Pn(4,:);Pn(5,:);Pn(6,:);Pn(7,:)];  %CASO NON SMORZATO
 
ff=PPn*F;    %determinante polinomio caratter*vettore modale
 
for i=1:n1
    FFn{i}=subs(ff,s,Sn(i));
end
 
%%  4   MATRICE MODALE
 
%%  4.1 FORMAZIONE DEL VETTORE MODALE DEL PRIMO MODO naturale
 
Fn011=1;
 
Fn0012=subs(FFn{1}(1),F1,Fn011);
Fn0012=solve(Fn0012,F2);
 
Fn0013=subs(FFn{1}(2),F1,Fn011);
Fn0013=subs(Fn0013,F2,Fn0012);
Fn0013=solve(Fn0013,F3);
 
Fn0014=subs(FFn{1}(3),F2,Fn0012);
Fn0014=subs(Fn0014,F3,Fn0013);
Fn0014=solve(Fn0014,F4);
 
Fn0015=subs(FFn{1}(4),F3,Fn0013);
Fn0015=subs(Fn0015,F4,Fn0014);
Fn0015=solve(Fn0015,F5);
 
Fn0016=subs(FFn{1}(5),F4,Fn0014);
Fn0016=subs(Fn0016,F5,Fn0015);
Fn0016=solve(Fn0016,F6);
 
Fn0017=subs(FFn{1}(6),F5,Fn0015);
Fn0017=subs(Fn0017,F6,Fn0016);
Fn0017=solve(Fn0017,F7);
 
% fprintf('Vettore modale del PRIMO modo naturale\n')
FFn1=[Fn011;Fn0012;Fn0013;Fn0014;Fn0015;Fn0016;Fn0017];     
FFn1=double(FFn1);     %VETTORE MODALE PRIMO MODO naturale
 
figure (1)
plot(FFn1,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('PRIMO modo naturale di vibrare')
grid on;
hold on
 
%%  4.2 FORMAZIONE DEL VETTORE MODALE DEL SECONDO MODO naturale
 
Fn021=1;
 
Fn0022=subs(FFn{2}(1),F1,Fn021);
Fn0022=solve(Fn0022,F2);
 
Fn0023=subs(FFn{2}(2),F1,Fn021);
Fn0023=subs(Fn0023,F2,Fn0022);
Fn0023=solve(Fn0023,F3);
 
Fn0024=subs(FFn{2}(3),F2,Fn0022);
Fn0024=subs(Fn0024,F3,Fn0023);
Fn0024=solve(Fn0024,F4);
 
Fn0025=subs(FFn{2}(4),F3,Fn0023);
Fn0025=subs(Fn0025,F4,Fn0024);
Fn0025=solve(Fn0025,F5);
 
Fn0026=subs(FFn{2}(5),F4,Fn0024);
Fn0026=subs(Fn0026,F5,Fn0025);
Fn0026=solve(Fn0026,F6);
 
Fn0027=subs(FFn{2}(6),F5,Fn0025);
Fn0027=subs(Fn0027,F6,Fn0026);
Fn0027=solve(Fn0027,F7);
 
% fprintf('Vettore modale del SECONDO modo naturale\n')
FFn2=[Fn021;Fn0022;Fn0023;Fn0024;Fn0025;Fn0026;Fn0027];     
FFn2=double(FFn2);     %VETTORE MODALE SECONDO MODO naturale
 
figure (2)
plot(FFn2,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('SECONDO modo naturale di vibrare')
grid on;
hold on
 
%%  4.3 FORMAZIONE DEL VETTORE MODALE DEL TERZO MODO naturale
 
Fn031=1;
 
Fn0032=subs(FFn{3}(1),F1,Fn031);
Fn0032=solve(Fn0032,F2);
 
Fn0033=subs(FFn{3}(2),F1,Fn031);
Fn0033=subs(Fn0033,F2,Fn0032);
Fn0033=solve(Fn0033,F3);
 
Fn0034=subs(FFn{3}(3),F2,Fn0032);
Fn0034=subs(Fn0034,F3,Fn0033);
Fn0034=solve(Fn0034,F4);
 
Fn0035=subs(FFn{3}(4),F3,Fn0033);
Fn0035=subs(Fn0035,F4,Fn0034);
Fn0035=solve(Fn0035,F5);
 
Fn0036=subs(FFn{3}(5),F4,Fn0034);
Fn0036=subs(Fn0036,F5,Fn0035);
Fn0036=solve(Fn0036,F6);
 
Fn0037=subs(FFn{3}(6),F5,Fn0035);
Fn0037=subs(Fn0037,F6,Fn0036);
Fn0037=solve(Fn0037,F7);
 
% fprintf('Vettore modale del TERZO modo naturale\n')
FFn3=[Fn031;Fn0032;Fn0033;Fn0034;Fn0035;Fn0036;Fn0037];     
FFn3=double(FFn3);     %VETTORE MODALE TERZO MODO naturale
 
figure (3)
plot(FFn3,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('TERZO modo naturale di vibrare')
grid on;
hold on
 
%%  4.4 FORMAZIONE DEL VETTORE MODALE DEL QUARTO MODO naturale
 
Fn041=1;
 
Fn0042=subs(FFn{4}(1),F1,Fn041);
Fn0042=solve(Fn0042,F2);
 
Fn0043=subs(FFn{4}(2),F1,Fn041);
Fn0043=subs(Fn0043,F2,Fn0042);
Fn0043=solve(Fn0043,F3);
 
Fn0044=subs(FFn{4}(3),F2,Fn0042);
Fn0044=subs(Fn0044,F3,Fn0043);
Fn0044=solve(Fn0044,F4);
 
Fn0045=subs(FFn{4}(4),F3,Fn0043);
Fn0045=subs(Fn0045,F4,Fn0044);
Fn0045=solve(Fn0045,F5);
 
Fn0046=subs(FFn{4}(5),F4,Fn0044);
Fn0046=subs(Fn0046,F5,Fn0045);
Fn0046=solve(Fn0046,F6);
 
Fn0047=subs(FFn{4}(6),F5,Fn0045);
Fn0047=subs(Fn0047,F6,Fn0046);
Fn0047=solve(Fn0047,F7);
 
% fprintf('Vettore modale del QUARTO modo naturale\n')
FFn4=[Fn041;Fn0042;Fn0043;Fn0044;Fn0045;Fn0046;Fn0047];     
FFn4=double(FFn4);     %VETTORE MODALE QUARTO MODO naturale
 
figure (4)
plot(FFn4,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('QUARTO modo naturale di vibrare')
grid on;
hold on
 
%%  4.5 FORMAZIONE DEL VETTORE MODALE DEL QUINTO MODO naturale
 
Fn051=1;
 
Fn0052=subs(FFn{5}(1),F1,Fn051);
Fn0052=solve(Fn0052,F2);
 
Fn0053=subs(FFn{5}(2),F1,Fn051);
Fn0053=subs(Fn0053,F2,Fn0052);
Fn0053=solve(Fn0053,F3);
 
Fn0054=subs(FFn{5}(3),F2,Fn0052);
Fn0054=subs(Fn0054,F3,Fn0053);
Fn0054=solve(Fn0054,F4);
 
Fn0055=subs(FFn{5}(4),F3,Fn0053);
Fn0055=subs(Fn0055,F4,Fn0054);
Fn0055=solve(Fn0055,F5);
 
Fn0056=subs(FFn{5}(5),F4,Fn0054);
Fn0056=subs(Fn0056,F5,Fn0055);
Fn0056=solve(Fn0056,F6);
 
Fn0057=subs(FFn{5}(6),F5,Fn0055);
Fn0057=subs(Fn0057,F6,Fn0056);
Fn0057=solve(Fn0057,F7);
 
% fprintf('Vettore modale del QUINTO modo naturale\n')
FFn5=[Fn051;Fn0052;Fn0053;Fn0054;Fn0055;Fn0056;Fn0057];     
FFn5=double(FFn5);     %VETTORE MODALE QUINTO MODO naturale
 
figure (5)
plot(FFn5,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('QUINTO modo naturale di vibrare')
grid on;
hold on
 
%%  4.6 FORMAZIONE DEL VETTORE MODALE DEL SESTO MODO naturale
 
Fn061=1;
 
Fn0062=subs(FFn{6}(1),F1,Fn061);
Fn0062=solve(Fn0062,F2);
 
Fn0063=subs(FFn{6}(2),F1,Fn061);
Fn0063=subs(Fn0063,F2,Fn0062);
Fn0063=solve(Fn0063,F3);
 
Fn0064=subs(FFn{6}(3),F2,Fn0062);
Fn0064=subs(Fn0064,F3,Fn0063);
Fn0064=solve(Fn0064,F4);
 
Fn0065=subs(FFn{6}(4),F3,Fn0063);
Fn0065=subs(Fn0065,F4,Fn0064);
Fn0065=solve(Fn0065,F5);
 
Fn0066=subs(FFn{6}(5),F4,Fn0064);
Fn0066=subs(Fn0066,F5,Fn0065);
Fn0066=solve(Fn0066,F6);
 
Fn0067=subs(FFn{6}(6),F5,Fn0065);
Fn0067=subs(Fn0067,F6,Fn0066);
Fn0067=solve(Fn0067,F7);
 
% fprintf('Vettore modale del SESTO modo naturale\n')
FFn6=[Fn061;Fn0062;Fn0063;Fn0064;Fn0065;Fn0066;Fn0067];     
FFn6=double(FFn6);     %VETTORE MODALE SESTO MODO naturale
 
figure (6)
plot(FFn6,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('SESTO modo naturale di vibrare')
grid on;
hold on
 
%%  4.7 FORMAZIONE DEL VETTORE MODALE DEL SETTIMO MODO naturale
 
Fn071=1;
 
Fn0072=subs(FFn{7}(1),F1,Fn071);
Fn0072=solve(Fn0072,F2);
 
Fn0073=subs(FFn{7}(2),F1,Fn071);
Fn0073=subs(Fn0073,F2,Fn0072);
Fn0073=solve(Fn0073,F3);
 
Fn0074=subs(FFn{7}(3),F2,Fn0072);
Fn0074=subs(Fn0074,F3,Fn0073);
Fn0074=solve(Fn0074,F4);
 
Fn0075=subs(FFn{7}(4),F3,Fn0073);
Fn0075=subs(Fn0075,F4,Fn0074);
Fn0075=solve(Fn0075,F5);
 
Fn0076=subs(FFn{7}(5),F4,Fn0074);
Fn0076=subs(Fn0076,F5,Fn0075);
Fn0076=solve(Fn0076,F6);
 
Fn0077=subs(FFn{7}(6),F5,Fn0075);
Fn0077=subs(Fn0077,F6,Fn0076);
Fn0077=solve(Fn0077,F7);
 
% fprintf('Vettore modale del SETTIMO modo naturale\n')
FFn7=[Fn071;Fn0072;Fn0073;Fn0074;Fn0075;Fn0076;Fn0077];     
FFn7=double(FFn7);     %VETTORE MODALE SETTIMO MODO naturale
 
figure (7)
plot(FFn7,H,'bo-','MarkerSize',2)
xlabel('amplificazione')
ylabel('altezza')
title('SETTIMO modo naturale di vibrare')
grid on;
hold on
 
%%  4.8 FORMAZIONE DELLA MATRICE MODALE E MATRICE SPETTRALE dei modi naturali
 
fprintf('Matrice modale dei modi naturali\n')
FFn=[FFn1,FFn2,FFn3,FFn4,FFn5,FFn6,FFn7];     
FFn=double(FFn)     %MATRICE MODALE naturale
 
fprintf('Matrice spettrale dei modi naturali\n')
OMEGAn=abs(Sn);
OMEGAn_quadr=OMEGAn.^2.*eye(n1)       %mattrice SPETTRALE naturale
 
fprintf('Periodi naturali di vibrare\n')
Tn=2*pi./OMEGAn
 
%%  5   MATRICI DIAGONALI DI MASSA E DI RIGIDEZZA' 
 
K=FFn'*k*FFn;
M=FFn'*m*FFn;
 
tt=[0:Tn(1)/10:3*Tn(1)]';
n2=length(tt);
 
SSn1=Sn(1);
SSn2=Sn(2);
SSn3=Sn(3);
SSn4=Sn(4);
SSn5=Sn(5);
SSn6=Sn(6);
SSn7=Sn(7);
 
Fn1=FFn1';
 
Fn1_1=Fn1(1);
Fn1_2=Fn1(2);
Fn1_3=Fn1(3);
Fn1_4=Fn1(4);
Fn1_5=Fn1(5);
Fn1_6=Fn1(6);
Fn1_7=Fn1(7);
 
Fn2=FFn2';
 
Fn2_1=Fn2(1);
Fn2_2=Fn2(2);
Fn2_3=Fn2(3);
Fn2_4=Fn2(4);
Fn2_5=Fn2(5);
Fn2_6=Fn2(6);
Fn2_7=Fn2(7);
 
Fn3=FFn3';
 
Fn3_1=Fn3(1);
Fn3_2=Fn3(2);
Fn3_3=Fn3(3);
Fn3_4=Fn3(4);
Fn3_5=Fn3(5);
Fn3_6=Fn3(6);
Fn3_7=Fn3(7);
 
Fn4=FFn4';
 
Fn4_1=Fn4(1);
Fn4_2=Fn4(2);
Fn4_3=Fn4(3);
Fn4_4=Fn4(4);
Fn4_5=Fn4(5);
Fn4_6=Fn4(6);
Fn4_7=Fn4(7);
 
Fn5=FFn5';
 
Fn5_1=Fn5(1);
Fn5_2=Fn5(2);
Fn5_3=Fn5(3);
Fn5_4=Fn5(4);
Fn5_5=Fn5(5);
Fn5_6=Fn5(6);
Fn5_7=Fn5(7);
 
Fn6=FFn6';
 
Fn6_1=Fn6(1);
Fn6_2=Fn6(2);
Fn6_3=Fn6(3);
Fn6_4=Fn6(4);
Fn6_5=Fn6(5);
Fn6_6=Fn6(6);
Fn6_7=Fn6(7);
 
Fn7=FFn7';
 
Fn7_1=Fn7(1);
Fn7_2=Fn7(2);
Fn7_3=Fn7(3);
Fn7_4=Fn7(4);
Fn7_5=Fn7(5);
Fn7_6=Fn7(6);
Fn7_7=Fn7(7);
 
%%  8.1 MATRICE DI SMORZAMENTO
 
%%  8.1.1   RAYLEIGHT
 
ksi=0.05;       % valore fisso - lo fissiamo noi - uguale per due modi e tutti i gdl
ksii1=ksi;
ksii7=ksi;
 
x0=ksi*(2*SSn1*SSn7)/(SSn1+SSn7);       %   coefficiente a0
 
x1=ksi*2/(SSn1+SSn7);       %   coefficiente a1
 
cc=x0*m+x1*k;                %   matrice di smorzamento proporzionale allo Rayleigh
CC=FFn'*cc*FFn;               
 
ksii3=x0/(2*SSn3)+(x1*SSn3)/2;           %   coef di smorz per gli altri modi di vibr vanno determinate con x0 e x1
ksii4=x0/(2*SSn4)+(x1*SSn4)/2;          
ksii5=x0/(2*SSn5)+(x1*SSn5)/2;          
ksii6=x0/(2*SSn6)+(x1*SSn6)/2;          
ksii2=x0/(2*SSn2)+(x1*SSn2)/2;          
 
KSII=[ksii1;ksii2;ksii3;ksii4;ksii5;ksii6;ksii7];
KSII0=abs(KSII);
Sn0=abs(Sn);
 
SSd1=SSn1*sqrt(1-ksi^2);
SSd2=SSn2*sqrt(1-ksi^2);
SSd3=SSn3*sqrt(1-ksi^2);
SSd4=SSn4*sqrt(1-ksi^2);
SSd5=SSn5*sqrt(1-ksi^2);
SSd6=SSn6*sqrt(1-ksi^2);
SSd7=SSn7*sqrt(1-ksi^2);
 
Sd=[SSd1;SSd2;SSd3;SSd4;SSd5;SSd6;SSd7];
 
%%  8.1.2  CAUGHEY
syms y0 y1 y2 y3 y4 y5 y6 
 
F0 = (y0/2)*SSn1^(-1)+y1*SSn1/2+y2*SSn1^3/2+y3*SSn1^5/2+y4*SSn1^7/2+y5*SSn1^9/2+y6*SSn1^11/2-ksi;
F1 = (y0/2)*SSn2^(-1)+y1*SSn2/2+y2*SSn2^3/2+y3*SSn2^5/2+y4*SSn2^7/2+y5*SSn2^9/2+y6*SSn2^11/2-ksi;
F2 = (y0/2)*SSn3^(-1)+y1*SSn3/2+y2*SSn3^3/2+y3*SSn3^5/2+y4*SSn3^7/2+y5*SSn3^9/2+y6*SSn3^11/2-ksi;
F3 = (y0/2)*SSn4^(-1)+y1*SSn4/2+y2*SSn4^3/2+y3*SSn4^5/2+y4*SSn4^7/2+y5*SSn4^9/2+y6*SSn4^11/2-ksi;
F4 = (y0/2)*SSn5^(-1)+y1*SSn5/2+y2*SSn5^3/2+y3*SSn5^5/2+y4*SSn5^7/2+y5*SSn5^9/2+y6*SSn5^11/2-ksi;
F5 = (y0/2)*SSn6^(-1)+y1*SSn6/2+y2*SSn6^3/2+y3*SSn6^5/2+y4*SSn6^7/2+y5*SSn6^9/2+y6*SSn6^11/2-ksi;
F6 = (y0/2)*SSn7^(-1)+y1*SSn7/2+y2*SSn7^3/2+y3*SSn7^5/2+y4*SSn7^7/2+y5*SSn7^9/2+y6*SSn7^11/2-ksi;
 
% Soluzione sistema di equazioni delle condizioni al contorno
% nelle 8 costanti di integrazione incognite
Sol = solve([F0,F1,F2,F3,F4,F5,F6],[y0,y1,y2,y3,y4,y5,y6]);
y0 = Sol.y0;   
y1 = Sol.y1;                        
y2 = Sol.y2;                        
y3 = Sol.y3;                        
y4 = Sol.y4;                        
y5 = Sol.y5;                        
y6 = Sol.y6;          
 
y0=double(y0);
y1=double(y1);
y2=double(y2);
y3=double(y3);
y4=double(y4);
y5=double(y5);
y6=double(y6);
 
c=y0*m+y1*k+y2*k*m^(-1)*k+y3*k^3*m^(-2)+y4*k^2*m^(-3)*k^2+y5*k^3*m^(-4)*k^2+y6*k^3*m^(-5)*k^3;
C=FFn'*c*FFn;
 
%%  8.2 SOLUZIONE DI EQUAZIONE DEL MOTO - sistema smorzato classicamente - oscillazioni libere
fprintf('SOLUZIONE DI EQUAZIONE DEL MOTO - sistema smorzato classicamente - oscillazioni libere\n')
fprintf('PRIMO modo - sistema smorzato classicamente - oscillazioni libere\n')
 
M1=M(1,1);  M2=M(2,2);  M3=M(3,3);  M4=M(4,4);  M5=M(5,5);  M6=M(6,6);  M7=M(7,7);
K1=K(1,1);  K2=K(2,2);  K3=K(3,3);  K4=K(4,4);  K5=K(5,5);  K6=K(6,6);  K7=K(7,7);
C1=C(1,1);  C2=C(2,2);  C3=C(3,3);  C4=C(4,4);  C5=C(5,5);  C6=C(6,6);  C7=C(7,7);
 
%%  9   SISTEMA CLASSICAMENTE SMORZATO - OSCILLAZIONI FORZATE
%%  9.1     FORZANTE + Qst
 
%   PRIMO modo
 
p0=ag*mpi;
p=p0*exp(1i*omega*ttt);
 
P01=p0/FFn1;
P02=p0/FFn2;
P03=p0/FFn3;
P04=p0/FFn4;
P05=p0/FFn5;
P06=p0/FFn6;
P07=p0/FFn7;
 
P1_1=p/Fn1_1;
P1_2=p/Fn1_2;
P1_3=p/Fn1_3;
P1_4=p/Fn1_4;
P1_5=p/Fn1_5;
P1_6=p/Fn1_6;
P1_7=p/Fn1_7;
 
P2_1=p/Fn2_1;
P2_2=p/Fn2_2;
P2_3=p/Fn2_3;
P2_4=p/Fn2_4;
P2_5=p/Fn2_5;
P2_6=p/Fn2_6;
P2_7=p/Fn2_7;
 
P3_1=p/Fn3_1;
P3_2=p/Fn3_2;
P3_3=p/Fn3_3;
P3_4=p/Fn3_4;
P3_5=p/Fn3_5;
P3_6=p/Fn3_6;
P3_7=p/Fn3_7;
 
P4_1=p/Fn4_1;
P4_2=p/Fn4_2;
P4_3=p/Fn4_3;
P4_4=p/Fn4_4;
P4_5=p/Fn4_5;
P4_6=p/Fn4_6;
P4_7=p/Fn4_7;
 
P5_1=p/Fn5_1;
P5_2=p/Fn5_2;
P5_3=p/Fn5_3;
P5_4=p/Fn5_4;
P5_5=p/Fn5_5;
P5_6=p/Fn5_6;
P5_7=p/Fn5_7;
 
P6_1=p/Fn6_1;
P6_2=p/Fn6_2;
P6_3=p/Fn6_3;
P6_4=p/Fn6_4;
P6_5=p/Fn6_5;
P6_6=p/Fn6_6;
P6_7=p/Fn6_7;
 
P7_1=p/Fn7_1;
P7_2=p/Fn7_2;
P7_3=p/Fn7_3;
P7_4=p/Fn7_4;
P7_5=p/Fn7_5;
P7_6=p/Fn7_6;
P7_7=p/Fn7_7;
 
qst1=P01/K;   
qst2=P02/K;   
qst3=P03/K;   
qst4=P04/K;   
qst5=P05/K;   
qst6=P06/K;   
qst7=P07/K;   
 
P01_1=P01(1);
P01_2=P01(2);
P01_3=P01(3);
P01_4=P01(4);
P01_5=P01(5);
P01_6=P01(6);
P01_7=P01(7);
 
bt1=omega/SSn1;
bt2=omega/SSn2;
bt3=omega/SSn3;
bt4=omega/SSn4;
bt5=omega/SSn5;
bt6=omega/SSn6;
bt7=omega/SSn7;
 
bt=[bt1,bt2,bt3,bt4,bt5,bt6,bt7];
 
qst1_1=qst1(1);
qst1_2=qst1(2);
qst1_3=qst1(3);
qst1_4=qst1(4);
qst1_5=qst1(5);
qst1_6=qst1(6);
qst1_7=qst1(7);
 
%%  9.2   SOLUZIONE DI EQUAZIONE DEL MOTO - sistema smorzato classicamente - oscillazioni forzate
fprintf('SOLUZIONE DI EQUAZIONE DEL MOTO - sistema smorzato classicamente - oscillazioni forzate\n')
fprintf('PRIMO modo - sistema smorzato classicamente - oscillazioni forzate\n')
 
%   ricerca costante di integrazione della risposta permanente
 
[A,qperm,AAA1,AAA2,AAA3,AAA4,AAA5,AAA6,AAA7,qperm1_1,qperm1_2,qperm1_3,qperm1_4,qperm1_5,qperm1_6,qperm1_7] = tes_8;
 
K1=K(1,1);
K2=K(2,2);
K3=K(3,3);
K4=K(4,4);
K5=K(5,5);
K6=K(6,6);
K7=K(7,7);
 
C1=C(1,1);
C2=C(2,2);
C3=C(3,3);
C4=C(4,4);
C5=C(5,5);
C6=C(6,6);
C7=C(7,7);
 
M1=M(1,1);
M2=M(2,2);
M3=M(3,3);
M4=M(4,4);
M5=M(5,5);
M6=M(6,6);
M7=M(7,7);
 
%   risposta permanente
%   PRIMO modo
 
qperm1_1=subs(qperm1_1);
qperm1_1=double(qperm1_1);
 
qperm1_2=subs(qperm1_2);
qperm1_2=double(qperm1_2);
 
qperm1_3=subs(qperm1_3);
qperm1_3=double(qperm1_3);
 
qperm1_4=subs(qperm1_4);
qperm1_4=double(qperm1_4);
 
qperm1_5=subs(qperm1_5);
qperm1_5=double(qperm1_5);
 
qperm1_6=subs(qperm1_6);
qperm1_6=double(qperm1_6);
 
qperm1_7=subs(qperm1_7);
qperm1_7=double(qperm1_7);
 
qperm1=[qperm1_1;qperm1_2;qperm1_3;qperm1_4;qperm1_5;qperm1_6;qperm1_7];
 
%   ricerca costante di integrazione della risposta transitoria
 
%%  9.2.1   PRIMO MODO
 
[HA1_1,HA1_2,HA1_3,HA1_4,HA1_5,HA1_6,HA1_7,Q1_1...
    ,V1_1,A1_1,Q1_2,V1_2,A1_2,Q1_3,V1_3...
    ,A1_3,Q1_4,V1_4,A1_4,Q1_5,V1_5,A1_5,Q1_6...
    ,V1_6,A1_6,Q1_7,V1_7,A1_7] = tes_9;
 
%%  CINEMATICA
 
%   CINEMATICA PRIMO MODO - PRIMO GDL
 
Q1_1=subs(Q1_1);
Q1_1=double(Q1_1);
 
V1_1=subs(V1_1);
V1_1=double(V1_1);
 
A1_1=subs(A1_1);
A1_1=double(A1_1);
 
%   CINEMATICA PRIMO MODO - SECONDO GDL
 
Q1_2=subs(Q1_2);
Q1_2=double(Q1_2);
 
V1_2=subs(V1_2);
V1_2=double(V1_2);
 
A1_2=subs(A1_2);
A1_2=double(A1_2);
 
%   CINEMATICA PRIMO MODO - TERZO GDL
 
Q1_3=subs(Q1_3);
Q1_3=double(Q1_3);
 
V1_3=subs(V1_3);
V1_3=double(V1_3);
 
A1_3=subs(A1_3);
A1_3=double(A1_3);
 
%   CINEMATICA PRIMO MODO - QUARTO GDL
 
Q1_4=subs(Q1_4);
Q1_4=double(Q1_4);
 
V1_4=subs(V1_4);
V1_4=double(V1_4);
 
A1_4=subs(A1_4);
A1_4=double(A1_4);
 
%   CINEMATICA PRIMO MODO - QUINTO GDL
 
Q1_5=subs(Q1_5);
Q1_5=double(Q1_5);
 
V1_5=subs(V1_5);
V1_5=double(V1_5);
 
A1_5=subs(A1_5);
A1_5=double(A1_5);
 
%   CINEMATICA PRIMO MODO - SESTO GDL
 
Q1_6=subs(Q1_6);
Q1_6=double(Q1_6);
 
V1_6=subs(V1_6);
V1_6=double(V1_6);
 
A1_6=subs(A1_6);
A1_6=double(A1_6);
 
%   CINEMATICA PRIMO MODO - SETTIMO GDL
 
Q1_7=subs(Q1_7);
Q1_7=double(Q1_7);
 
V1_7=subs(V1_7);
V1_7=double(V1_7);
 
A1_7=subs(A1_7);
A1_7=double(A1_7);
 
%   CINEMATICA PRIMO MODO - COMPLETO
 
QQ1=[Q1_1,Q1_2,Q1_3,Q1_4,Q1_5,Q1_6,Q1_7];
VV1=[V1_1,V1_2,V1_3,V1_4,V1_5,V1_6,V1_7];
AA1=[A1_1,A1_2,A1_3,A1_4,A1_5,A1_6,A1_7];
 
%%   CONTROLLO DI EQUILIBRIO
 
% EEQQ1_1=Fn1_1*m11*A1_1+Fn1_1*k11*Q1_1+Fn1_1*c11*V1_1+ag*m11*exp(1i*omega*ttt);
% EEQQ1_1=A1_1+SSn1^2*Q1_1+2*ksi*SSn1*V1_1+Fn1_1*ag*m11*exp(1i*omega*ttt)/M1;
EEQQ1_1=M1*A1_1+K1*Q1_1+C1*V1_1+P1_1;
 
if EEQQ1_1<1e-5
        fprintf('Equilibrio del PRIMO modo al PRIMO GDL - soddisfatto\n')
else
        fprintf('Equilibrio del PRIMO modo al PRIMO GDL - non soddisfatto\n')
end

%   CALCOLO SIMBOLICO   
%   TESINA
%   OSCILLAZIONI FORZATE - SISTEMA SMORZATO CLASSICAMENTE
%   ricerca della costante di integrazione A
 
function [A,qperm,AAA1,AAA2,AAA3,AAA4,AAA5,AAA6,AAA7,qperm1_1,qperm1_2,qperm1_3,qperm1_4,qperm1_5,qperm1_6,qperm1_7] = tes_8(~,~,~)
 
format short
 
syms m c b1 b2  tt   qperm11 vperm11 aperm11 aperm qperm vperm qst ksi p0 k  omega bt A ccr bt qst p0
syms bt1 bt2 bt3 bt4 bt5 bt6 bt7 ksi SSn qst1 qst2 qst3 qst4 qst5 qst6 qst7 C K m k c
syms SSn1 SSn2 SSn3 SSn4 SSn5 SSn6 SSn7 ttt P0 M C K FFn FFn1 FFn2 FFn3 FFn4 FFn5 FFn6 FFn7
syms qst1_1 qst1_2 qst1_3 qst1_4 qst1_5 qst1_6 qst1_7 m21 m22 m23 m24 m25 m26 m27
syms qst2_1 qst2_2 qst2_3 qst2_4 qst2_5 qst2_6 qst2_7 m11 m12 m13 m14 m15 m16 m17
syms qst3_1 qst3_2 qst3_3 qst3_4 qst3_5 qst3_6 qst3_7 m31 m32 m33 m34 m35 m36 m37
syms qst4_1 qst4_2 qst4_3 qst4_4 qst4_5 qst4_6 qst4_7 m41 m42 m43 m44 m45 m46 m47
syms qst5_1 qst5_2 qst5_3 qst5_4 qst5_5 qst5_6 qst5_7 m51 m52 m53 m54 m55 m56 m57
syms qst6_1 qst6_2 qst6_3 qst6_4 qst6_5 qst6_6 qst6_7 m61 m62 m63 m64 m65 m66 m67
syms qst7_1 qst7_2 qst7_3 qst7_4 qst7_5 qst7_6 qst7_7 c11 c12 c13 c14 c15 c16 c17
syms Fn1_1 Fn1_2 Fn1_3 Fn1_4 Fn1_5 Fn1_6 Fn1_7 c21 c22 c23 c24 c25 c26 c27
syms Fn2_1 Fn2_2 Fn2_3 Fn2_4 Fn2_5 Fn2_6 Fn2_7 c31 c32 c33 c34 c35 c36 c37
syms Fn3_1 Fn3_2 Fn3_3 Fn3_4 Fn3_5 Fn3_6 Fn3_7 c41 c42 c43 c44 c45 c46 c47
syms Fn4_1 Fn4_2 Fn4_3 Fn4_4 Fn4_5 Fn4_6 Fn4_7 c51 c52 c53 c54 c55 c56 c57
syms Fn5_1 Fn5_2 Fn5_3 Fn5_4 Fn5_5 Fn5_6 Fn5_7 c61 c62 c63 c64 c65 c66 c67
syms Fn6_1 Fn6_2 Fn6_3 Fn6_4 Fn6_5 Fn6_6 Fn6_7 c71 c72 c73 c74 c75 c76 c77
syms Fn7_1 Fn7_2 Fn7_3 Fn7_4 Fn7_5 Fn7_6 Fn7_7 k11 k12 k13 k14 k15 k16 k17
syms k21 k22 k23 k24 k25 k26 k27 k31 k32 k33 k34 k35 k36 k37 k41 k42 k43 k44 
syms k45 k46 k47 k51 k52 k53 k54 k55 k56 k57 k61 k62 k63 k64 k65 k66 k67
syms k71 k72 k73 k74 k75 k76 k77 m71 m72 m73 m74 m75 m76 m77
syms K21 K22 K23 K24 K25 K26 K27 K31 K32 K33 K34 K35 K36 K37 K41 K42 K43 K44 
syms K45 K46 K47 K51 K52 K53 K54 K55 K56 K57 K61 K62 K63 K64 K65 K66 K67
syms K71 K72 K73 K74 K75 K76 K77 C11 C12 C13 C14 C15 C16 C17
syms C21 C22 C23 C24 C25 C26 C27 C31 C32 C33 C34 C35 C36 C37
syms C41 C42 C43 C44 C45 C46 C47 C51 C52 C53 C54 C55 C56 C57
syms C61 C62 C63 C64 C65 C66 C67 C71 C72 C73 C74 C75 C76 C77
syms K11 K12 K13 K14 K15 K16 K17 K1 K2 K3 K4 K5 K6 K7 C1 C2 C3 C4 C5 C6 C7
syms M1 M2 M3 M4 M5 M6 M7 mpi P01 P02 P03 P04 P05 P06 P07 P01_1 P01_2 P01_3 P01_4 P01_5 P01_6 P01_7
 
Qperm11=A*exp(1i*omega.*tt);
Vperm11=diff(Qperm11,tt,1);
Aperm11=diff(Qperm11,tt,2);
 
%   EQS scritta sulla base delle coord reali
EQS=P0*exp(1i*omega*tt)+M*Aperm11+C*Vperm11+K*Qperm11;
 
% Condizioni al contorno
f3 = subs(EQS);
 
% Soluzione sistema di equazioni delle condizioni al contorno
% nelle 8 costanti di integrazione incognite
sol = solve(f3,A);
A=simplify(sol);
A=subs(A,C,FFn'*c*FFn);
A=subs(A,c,ksi*2*m*SSn);
A=subs(A,M,FFn'*m*FFn);
A=subs(A,m,k/SSn^2);
A=subs(A,omega,bt*SSn);
A=subs(A,P0,qst*K);
A=subs(A,K,FFn'*k*FFn);
 
A=subs(A,k,(FFn')^(-1)*K*FFn^(-1));
A=subs(A,c,(FFn')^(-1)*C*FFn^(-1));
A=subs(A,m,(FFn')^(-1)*M*FFn^(-1));
 
A=simplify(A);
 
% Sostituzione delle costanti nelle funzioni di spostamento
AAA1=subs(A,bt,bt1);
AAA1=subs(AAA1,qst,qst1);
AAA1=subs(AAA1,FFn,FFn1);
AAA1=subs(AAA1,K,K1);
AAA1=subs(AAA1,C,C1);
AAA1=subs(AAA1,M,M1);
AAA1=subs(AAA1,P0,P01);
 
AAA2=subs(A,bt,bt2);
AAA2=subs(AAA2,qst,qst2);
AAA2=subs(AAA2,FFn,FFn2);
AAA2=subs(AAA2,K,K2);
AAA2=subs(AAA2,C,C2);
AAA2=subs(AAA2,M,M2);
AAA2=subs(AAA2,P0,P02);
 
AAA3=subs(A,bt,bt3);
AAA3=subs(AAA3,qst,qst3);
AAA3=subs(AAA3,FFn,FFn3);
AAA3=subs(AAA3,K,K3);
AAA3=subs(AAA3,C,C3);
AAA3=subs(AAA3,M,M3);
AAA3=subs(AAA3,P0,P03);
 
AAA4=subs(A,bt,bt4);
AAA4=subs(AAA4,qst,qst4);
AAA4=subs(AAA4,FFn,FFn4);
AAA4=subs(AAA4,K,K4);
AAA4=subs(AAA4,C,C4);
AAA4=subs(AAA4,M,M4);
AAA4=subs(AAA4,P0,P04);
 
AAA5=subs(A,bt,bt5);
AAA5=subs(AAA5,qst,qst5);
AAA5=subs(AAA5,FFn,FFn5);
AAA5=subs(AAA5,K,K5);
AAA5=subs(AAA5,C,C5);
AAA5=subs(AAA5,M,M5);
AAA5=subs(AAA5,P0,P05);
 
AAA6=subs(A,bt,bt6);
AAA6=subs(AAA6,qst,qst6);
AAA6=subs(AAA6,FFn,FFn6);
AAA6=subs(AAA6,K,K6);
AAA6=subs(AAA6,C,C6);
AAA6=subs(AAA6,M,M6);
AAA6=subs(AAA6,P0,P06);
 
AAA7=subs(A,bt,bt7);
AAA7=subs(AAA7,qst,qst7);
AAA7=subs(AAA7,FFn,FFn7);
AAA7=subs(AAA7,K,K7);
AAA7=subs(AAA7,C,C7);
AAA7=subs(AAA7,M,M7);
AAA7=subs(AAA7,P0,P07);
 
qperm = subs(Qperm11);
qperm=subs(qperm,k,(FFn')^(-1)*K*FFn^(-1));
qperm=subs(qperm,c,(FFn')^(-1)*C*FFn^(-1));
qperm=subs(qperm,m,(FFn')^(-1)*M*FFn^(-1));
 
    %   PRIMO modo
    
qperm1=subs(qperm,bt,bt1);
qperm1=subs(qperm1,SSn,SSn1);               
qperm1=subs(qperm1,FFn,FFn1);               
qperm1=subs(qperm1,qst,qst1);  
qperm1=subs(qperm1,K,K1);
qperm1=subs(qperm1,C,C1);
qperm1=subs(qperm1,M,M1);
qperm1=subs(qperm1,P0,P01);
 
qperm1_1=subs(qperm1,qst1,qst1_1);
qperm1_1=subs(qperm1_1,FFn1,Fn1_1);
qperm1_1=subs(qperm1_1,P01,P01_1);
 
qperm1_2=subs(qperm1,qst1,qst1_2);
qperm1_2=subs(qperm1_2,FFn1,Fn1_2);
qperm1_2=subs(qperm1_2,P01,P01_2);
 
qperm1_3=subs(qperm1,qst1,qst1_3);
qperm1_3=subs(qperm1_3,FFn1,Fn1_3);
qperm1_3=subs(qperm1_3,P01,P01_3);
 
qperm1_4=subs(qperm1,qst1,qst1_4);
qperm1_4=subs(qperm1_4,FFn1,Fn1_4);
qperm1_4=subs(qperm1_4,P01,P01_4);
 
qperm1_5=subs(qperm1,qst1,qst1_5);
qperm1_5=subs(qperm1_5,FFn1,Fn1_5);
qperm1_5=subs(qperm1_5,P01,P01_5);
 
qperm1_6=subs(qperm1,qst1,qst1_6);
qperm1_6=subs(qperm1_6,FFn1,Fn1_6);
qperm1_6=subs(qperm1_6,P01,P01_6);
 
qperm1_7=subs(qperm1,qst1,qst1_7);
qperm1_7=subs(qperm1_7,FFn1,Fn1_7);
qperm1_7=subs(qperm1_7,P01,P01_7);
 
end

%   CALCOLO SIMBOLICO   
%   TESINA
%   OSCILLAZIONI FORZATE - SISTEMA SMORZATO CLASSICAMENTE
%   ricerca di: b1 + b2 + HA + cinematica + qtrans
%   PRIMO modo
 
function [HA1_1,HA1_2,HA1_3,HA1_4,HA1_5,HA1_6,HA1_7,Q1_1...
    ,V1_1,A1_1,Q1_2,V1_2,A1_2,Q1_3,V1_3...
    ,A1_3,Q1_4,V1_4,A1_4,Q1_5,V1_5,A1_5,Q1_6...
    ,V1_6,A1_6,Q1_7,V1_7,A1_7] = tes_9(~,~,~)
 
format short
 
syms B1 B2 tt qst ksi p0   omega bt qperm qtrans A ttt  c 
syms FFn m k Q0 V0 Fn1 Fn2 Fn3 Fn4 Fn5 Fn6 Fn7 mpi  SSn1 SSn SSd SSd1 Fn1_1 Fn1_2 Fn1_3 Fn1_4 Fn1_5 Fn1_6 Fn1_7 
syms bt1 bt2 bt3 bt4 bt5 bt6 bt7 ksi qst1 qst2 qst3 qst4 qst5 qst6 qst7 C K M1 M2 M3 M4 M5 M6 M7 M
syms qst1_1 qst1_2 qst1_3 qst1_4 qst1_5 qst1_6 qst1_7 K1 C1 K2 C2 K3 C3 K4 C4 K5 C5 K6 C6 K7 C7
syms k11 k12 k13 k14 k15 k16 k17 p  m22 p0 omega tt mpi m qst ag
syms m11 m12 m13 m14 m15 m16 m17 c11 c12 c13 c14 c15 c16 c17
syms P0 P01 P02 P03 P04 P05 P06 P07 P01_1 P01_2 P01_3 P01_4 P01_5 P01_6 P01_7
 
qtrans=exp(-SSn*ksi*tt).*(B1*cos(SSd*tt)+B2*sin(SSd*tt));
 
[А,qperm,~,~,~,~,~,~,~,~,~,~,~,~,~,~] = tes_8;
 
q11(tt)=qtrans+2*qperm;
% q11(tt)=qtrans+P0/K;
 
v11(tt)=diff(q11,tt,1);
a11(tt)=diff(q11,tt,2);
HD=X/qst;
 
% Condizioni al contorno
f1 = subs(q11,tt,0)-FFn*m*Q0;                           % q(t=0) = q0
f2 = subs(v11,tt,0)-FFn*m*V0;                           % v(t=0) = v0
 
% Soluzione sistema di equazioni delle condizioni al contorno
% nelle 8 costanti di integrazione incognite
Sol = solve([f2,f1],[B1,B2]);
B1 = Sol.B1;                          %   SOLUZIONE GENERALE
B1=subs(B1,m,mpi);
 
B1_1=subs(B1,FFn,Fn1);          
B1_1=subs(B1_1,qst,qst1);          
B1_1=subs(B1_1,bt,bt1);
B1_1=subs(B1_1,SSn,SSn1);             %   SOLUZIONE PER PRIMO MODO DI VIBRARE
B1_1=subs(B1_1,K,K1);
B1_1=subs(B1_1,C,C1);
B1_1=subs(B1_1,M,M1);
B1_1=subs(B1_1,P0,P01);
 
B1_1_1=subs(B1_1,qst1,qst1_1);
B1_1_1=subs(B1_1_1,Fn1,Fn1_1);
B1_1_1=subs(B1_1_1,P01,P01_1);
 
B1_1_2=subs(B1_1,qst1,qst1_2);
B1_1_2=subs(B1_1_2,Fn1,Fn1_2);
B1_1_2=subs(B1_1_2,P01,P01_2);
 
B1_1_3=subs(B1_1,qst1,qst1_3);
B1_1_3=subs(B1_1_3,Fn1,Fn1_3);
B1_1_3=subs(B1_1_3,P01,P01_3);
 
B1_1_4=subs(B1_1,qst1,qst1_4);
B1_1_4=subs(B1_1_4,Fn1,Fn1_4);
B1_1_4=subs(B1_1_4,P01,P01_4);
 
B1_1_5=subs(B1_1,qst1,qst1_5);
B1_1_5=subs(B1_1_5,Fn1,Fn1_5);
B1_1_5=subs(B1_1_5,P01,P01_5);
 
B1_1_6=subs(B1_1,qst1,qst1_6);
B1_1_6=subs(B1_1_6,Fn1,Fn1_6);
B1_1_6=subs(B1_1_6,P01,P01_6);
 
B1_1_7=subs(B1_1,qst1,qst1_7);
B1_1_7=subs(B1_1_7,Fn1,Fn1_7);
B1_1_7=subs(B1_1_7,P01,P01_7);
 
B2 = Sol.B2;
 
B2=subs(B2,m,mpi);
 
B2_1=subs(B2,SSn,SSn1);
B2_1=subs(B2_1,SSd,SSd1);
B2_1=subs(B2_1,qst,qst1);          
B2_1=subs(B2_1,FFn,Fn1);          
B2_1=subs(B2_1,m,mpi);
B2_1=subs(B2_1,bt,bt1);         %   SOLUZIONE PER PRIMO MODO DI VIBRARE
B2_1=subs(B2_1,K,K1);
B2_1=subs(B2_1,C,C1);
B2_1=subs(B2_1,M,M1);
B2_1=subs(B2_1,P0,P01);
 
B2_1_1=subs(B2_1,qst1,qst1_1);
B2_1_1=subs(B2_1_1,Fn1,Fn1_1);
B2_1_1=subs(B2_1_1,P01,P01_1);
 
B2_1_2=subs(B2_1,qst1,qst1_2);
B2_1_2=subs(B2_1_2,Fn1,Fn1_2);
B2_1_2=subs(B2_1_2,P01,P01_2);
 
B2_1_3=subs(B2_1,qst1,qst1_3);
B2_1_3=subs(B2_1_3,Fn1,Fn1_3);
B2_1_3=subs(B2_1_3,P01,P01_3);
 
B2_1_4=subs(B2_1,qst1,qst1_4);
B2_1_4=subs(B2_1_4,Fn1,Fn1_4);
B2_1_4=subs(B2_1_4,P01,P01_4);
 
B2_1_5=subs(B2_1,qst1,qst1_5);
B2_1_5=subs(B2_1_5,Fn1,Fn1_5);
B2_1_5=subs(B2_1_5,P01,P01_5);
 
B2_1_6=subs(B2_1,qst1,qst1_6);
B2_1_6=subs(B2_1_6,Fn1,Fn1_6);
B2_1_6=subs(B2_1_6,P01,P01_6);
 
B2_1_7=subs(B2_1,qst1,qst1_7);
B2_1_7=subs(B2_1_7,Fn1,Fn1_7);
B2_1_7=subs(B2_1_7,P01,P01_7);
 
% Sostituzione delle costanti nelle funzioni di spostamento
%%  GENERALE + PRIMO GDL
 
qtrans1 = subs(qtrans);
qtrans1=subs(qtrans1,bt,bt1);
qtrans1=subs(qtrans1,FFn,Fn1);
qtrans1=subs(qtrans1,qst,qst1);
qtrans1=subs(qtrans1,SSn,SSn1);
qtrans1=subs(qtrans1,SSd,SSd1);
qtrans1=subs(qtrans1,m,mpi);
qtrans1=subs(qtrans1,K,K1);
qtrans1=subs(qtrans1,C,C1);
qtrans1=subs(qtrans1,M,M1);              
qtrans1=subs(qtrans1,P0,P01);
 
qtrans1_1=subs(qtrans1,qst1,qst1_1);
qtrans1_1=subs(qtrans1_1,Fn1,Fn1_1);
qtrans1_1=subs(qtrans1_1,P01,P01_1);
 
Q1 = subs(q11);
Q1=subs(Q1,bt,bt1);
Q1=subs(Q1,FFn,Fn1);
Q1=subs(Q1,qst,qst1);
Q1=subs(Q1,SSn,SSn1);
Q1=subs(Q1,SSd,SSd1);
Q1=subs(Q1,m,mpi);
Q1=subs(Q1,K,K1);
Q1=subs(Q1,C,C1);
Q1=subs(Q1,M,M1);
Q1=subs(Q1,P0,P01);
 
Q1_1=subs(Q1,qst1,qst1_1);
Q1_1=subs(Q1_1,Fn1,Fn1_1);
Q1_1=subs(Q1_1,P01,P01_1);
 
V1 = subs(v11);
V1=subs(V1,bt,bt1);
V1=subs(V1,FFn,Fn1);
V1=subs(V1,qst,qst1);
V1=subs(V1,SSn,SSn1);
V1=subs(V1,SSd,SSd1);
V1=subs(V1,m,mpi);
V1=subs(V1,K,K1);
V1=subs(V1,C,C1);
V1=subs(V1,M,M1);
V1=subs(V1,P0,P01);
 
V1_1=subs(V1,qst1,qst1_1);
V1_1=subs(V1_1,Fn1,Fn1_1);
V1_1=subs(V1_1,P01,P01_1);
 
fprintf('V1_1=%s\n',V1_1)           %   SOLUZIONE PER PRIMO MODO DI VIBRARE PER PRIMO GDL!!!!!
 
A1 = subs(a11);
A1=subs(A1,bt,bt1);
A1=subs(A1,FFn,Fn1);
A1=subs(A1,qst,qst1);
A1=subs(A1,SSn,SSn1);
A1=subs(A1,SSd,SSd1);
A1=subs(A1,m,mpi);
A1=subs(A1,K,K1);
A1=subs(A1,C,C1);
A1=subs(A1,M,M1);
A1=subs(A1,P0,P01);
 
A1_1=subs(A1,qst1,qst1_1);
A1_1=subs(A1_1,Fn1,Fn1_1);
A1_1=subs(A1_1,P01,P01_1);
 
HD= subs(HD);
HV=bt.*HD;
 
HA=bt.^2.*HD;
HA=subs(HA,m,mpi);
 
HA1=subs(HA,bt,bt1);
HA1=subs(HA1,FFn,Fn1);
HA1=subs(HA1,qst,qst1);
HA1=subs(HA1,K,K1);
HA1=subs(HA1,C,C1);
HA1=subs(HA1,M,M1);
HA1=subs(HA1,P0,P01);
 
HA1_1=subs(HA1,Fn1,Fn1_1);
HA1_1=subs(HA1_1,qst1,qst1_1);
HA1_1=subs(HA1_1,P01,P01_1);
 
%%  SECONDO GDL
 
Q1_2=subs(Q1,Fn1,Fn1_2);
Q1_2=subs(Q1_2,qst1,qst1_2);
Q1_2=subs(Q1_2,P01,P01_2);
 
V1_2=subs(V1,Fn1,Fn1_2);
V1_2=subs(V1_2,qst1,qst1_2);
V1_2=subs(V1_2,P01,P01_2);
 
A1_2=subs(A1,Fn1,Fn1_2);
A1_2=subs(A1_2,qst1,qst1_2);
A1_2=subs(A1_2,P01,P01_2);
 
qtrans1_2=subs(qtrans1,Fn1,Fn1_2);
qtrans1_2=subs(qtrans1_2,qst1,qst1_2);
qtrans1_2=subs(qtrans1_2,P01,P01_2);
 
HA1_2=subs(HA1,Fn1,Fn1_2);
HA1_2=subs(HA1_2,qst1,qst1_2);
HA1_2=subs(HA1_2,P01,P01_2);
 
%%  TERZO GDL
 
Q1_3=subs(Q1,Fn1,Fn1_3);
Q1_3=subs(Q1_3,qst1,qst1_3);
Q1_3=subs(Q1_3,P01,P01_3);
 
V1_3=subs(V1,Fn1,Fn1_3);
V1_3=subs(V1_3,qst1,qst1_3);
V1_3=subs(V1_3,P01,P01_3);
 
A1_3=subs(A1,Fn1,Fn1_3);
A1_3=subs(A1_3,qst1,qst1_3);
A1_3=subs(A1_3,P01,P01_3);
 
qtrans1_3=subs(qtrans1,Fn1,Fn1_3);
qtrans1_3=subs(qtrans1_3,qst1,qst1_3);
qtrans1_3=subs(qtrans1_3,P01,P01_3);
 
HA1_3=subs(HA1,Fn1,Fn1_3);
HA1_3=subs(HA1_3,qst1,qst1_3);
HA1_3=subs(HA1_3,P01,P01_3);
 
%%  QUARTO GDL
 
Q1_4=subs(Q1,Fn1,Fn1_4);
Q1_4=subs(Q1_4,qst1,qst1_4);
Q1_4=subs(Q1_4,P01,P01_4);
 
V1_4=subs(V1,Fn1,Fn1_4);
V1_4=subs(V1_4,qst1,qst1_4);
V1_4=subs(V1_4,P01,P01_4);
 
A1_4=subs(A1,Fn1,Fn1_4);
A1_4=subs(A1_4,qst1,qst1_4);
A1_4=subs(A1_4,P01,P01_4);
 
qtrans1_4=subs(qtrans1,Fn1,Fn1_4);
qtrans1_4=subs(qtrans1_4,qst1,qst1_4);
qtrans1_4=subs(qtrans1_4,P01,P01_4);
 
HA1_4=subs(HA1,Fn1,Fn1_4);
HA1_4=subs(HA1_4,qst1,qst1_4);
HA1_4=subs(HA1_4,P01,P01_4);
 
%%  QUINTO GDL
 
Q1_5=subs(Q1,Fn1,Fn1_5);
Q1_5=subs(Q1_5,qst1,qst1_5);
Q1_5=subs(Q1_5,P01,P01_5);
 
V1_5=subs(V1,Fn1,Fn1_5);
V1_5=subs(V1_5,qst1,qst1_5);
V1_5=subs(V1_5,P01,P01_5);
 
A1_5=subs(A1,Fn1,Fn1_5);
A1_5=subs(A1_5,qst1,qst1_5);
A1_5=subs(A1_5,P01,P01_5);
 
qtrans1_5=subs(qtrans1,Fn1,Fn1_5);
qtrans1_5=subs(qtrans1_5,qst1,qst1_5);
qtrans1_5=subs(qtrans1_5,P01,P01_5);
 
HA1_5=subs(HA1,Fn1,Fn1_5);
HA1_5=subs(HA1_5,qst1,qst1_5);
HA1_5=subs(HA1_5,P01,P01_5);
 
%%  SESTO GDL
 
Q1_6=subs(Q1,Fn1,Fn1_6);
Q1_6=subs(Q1_6,qst1,qst1_6);
Q1_6=subs(Q1_6,P01,P01_6);
 
V1_6=subs(V1,Fn1,Fn1_6);
V1_6=subs(V1_6,qst1,qst1_6);
V1_6=subs(V1_6,P01,P01_6);
 
A1_6=subs(A1,Fn1,Fn1_6);
A1_6=subs(A1_6,qst1,qst1_6);
A1_6=subs(A1_6,P01,P01_6);
 
qtrans1_6=subs(qtrans1,Fn1,Fn1_6);
qtrans1_6=subs(qtrans1_6,qst1,qst1_6);
qtrans1_6=subs(qtrans1_6,P01,P01_6);
 
HA1_6=subs(HA1,Fn1,Fn1_6);
HA1_6=subs(HA1_6,qst1,qst1_6);
HA1_6=subs(HA1_6,P01,P01_6);
 
%%  SETTIMO GDL
 
Q1_7=subs(Q1,Fn1,Fn1_7);
Q1_7=subs(Q1_7,qst1,qst1_7);
Q1_7=subs(Q1_7,P01,P01_7);
 
V1_7=subs(V1,Fn1,Fn1_7);
V1_7=subs(V1_7,qst1,qst1_7);
V1_7=subs(V1_7,P01,P01_7);
 
A1_7=subs(A1,Fn1,Fn1_7);
A1_7=subs(A1_7,qst1,qst1_7);
A1_7=subs(A1_7,P01,P01_7);
 
qtrans1_7=subs(qtrans1,Fn1,Fn1_7);
qtrans1_7=subs(qtrans1_7,qst1,qst1_7);
qtrans1_7=subs(qtrans1_7,P01,P01_7);
 
HA1_7=subs(HA1,Fn1,Fn1_7);
HA1_7=subs(HA1_7,qst1,qst1_7);
HA1_7=subs(HA1_7,P01,P01_7);
 
end

clear 
clc 
close all
 
%%  1   DATI DI INPUT
 
mpi=1;              %massa del piano
h=3;                %altezza interpiano
EI=10;              %rigidezza flessionale
Q0=15;              %spostameno modale iniziale
V0=5;               %velocita' modale iniziale
ag=0.8;             %accelerazione al suolo
omega=5;            %frequenza del terremoto
%%  2   MATRICE DI MASSA E DI RIGIDEZZA'
 
kpi=24*EI/h^3;      %rigidezza flessionale del singolo piano
m=mpi*eye(7);      %matrice delle masse
k=([[2*kpi,-kpi,0,0,0,0,0];[-kpi,2*kpi,-kpi,0,0,0,0];[0,-kpi,2*kpi,-kpi,0,0,0];[0,0,-kpi,2*kpi,-kpi,0,0];[0,0,0,-kpi,2*kpi,-kpi,0];[0,0,0,0,-kpi,2*kpi,-kpi];[0,0,0,0,0,-kpi,kpi]]);
 
SSn1=-0.6233 + 0.0000i;
Fn1_3=2.8271 - 0.0000i;
Tn = [10.0807;3.4099;2.1074;1.5748;1.3025;1.1534;1.0773];
K1 =   33.7016 + 0.0000i;
C1 =  -5.4071 + 0.0000i;
M1 =   86.7508;
ksi=0.05;
SSd1=SSn1*sqrt(1-ksi^2);
tt=10;
p0=ag*mpi;
p=p0*exp(1i*omega*tt);
P01_3=p0*Fn1_3;
P1_3=p*Fn1_3;
qst1_3=P01_3/K1;   
bt1=omega/SSn1;
 
[AAAA1_3,qperm1_3,EQS1_3] = tes_1000;
 
AAAA1_3=subs(AAAA1_3);
AAAA1_3=double(AAAA1_3);
   
qperm1_3=subs(qperm1_3);
qperm1_3=double(qperm1_3);
 
EQS1_3=subs(EQS1_3);
EQS1_3=double(EQS1_3);
 
[BB1_1_3,BB2_1_3,QQ1_3,VV1_3,AA1_3,ff1,ff2] = tes_1001;
 
QQ1_3=subs(QQ1_3);
QQ1_3=double(QQ1_3);
 
VV1_3=subs(VV1_3);
VV1_3=double(VV1_3);
 
AA1_3=subs(AA1_3);
AA1_3=double(AA1_3);
 
ff1=subs(ff1);
ff1=subs(ff1);
 
ff2=subs(ff2);
ff2=subs(ff2);
 
BIL1_3=M1*AA1_3+K1*QQ1_3+C1*VV1_3+P1_3;

function [AAAA1_3,qperm1_3,EQS1_3] = tes_1000(~,~,~)
 
format short
 
syms tt omega AAAA ccr bt1  ksi  m k c SSn1 qst1_3 Fn1_3 
syms K1 C1 M1 mpi P01_3 
 
Qperm=AAAA*exp(1i*omega*tt);
Vperm=diff(Qperm,tt,1);
Aperm=diff(Qperm,tt,2);
 
%   EQS scritta sulla base delle coord reali
F1000=P01_3*exp(1i*omega*tt)+M1*Aperm+C1*Vperm+K1*Qperm;
 
% Condizioni al contorno
ff3 = subs(F1000);
 
% Soluzione sistema di equazioni delle condizioni al contorno
% nelle 8 costanti di integrazione incognite
sol = solve(ff3,AAAA);
AAAA=simplify(sol);
 
AAAA1_3=subs(AAAA,C1,Fn1_3*c*Fn1_3);
AAAA1_3=subs(AAAA1_3,c,ksi*2*m*SSn1);
AAAA1_3=subs(AAAA1_3,M1,Fn1_3*m*Fn1_3);
AAAA1_3=subs(AAAA1_3,m,k/SSn1^2);
AAAA1_3=subs(AAAA1_3,omega,bt1*SSn1);
AAAA1_3=subs(AAAA1_3,K1,Fn1_3*k*Fn1_3);
AAAA1_3=subs(AAAA1_3,k,(Fn1_3)^(-1)*K1*Fn1_3^(-1));
AAAA1_3=subs(AAAA1_3,c,(Fn1_3)^(-1)*C1*Fn1_3^(-1));
AAAA1_3=subs(AAAA1_3,m,(Fn1_3)^(-1)*M1*Fn1_3^(-1));
AAAA1_3=simplify(AAAA1_3);
 
qperm1_3 = subs(Qperm);
qperm1_3=subs(qperm1_3,omega,bt1*SSn1);
qperm1_3=subs(qperm1_3,k,(Fn1_3)^(-1)*K1*Fn1_3^(-1));
qperm1_3=subs(qperm1_3,c,(Fn1_3)^(-1)*C1*Fn1_3^(-1));
qperm1_3=subs(qperm1_3,m,(Fn1_3)^(-1)*M1*Fn1_3^(-1));
 
vperm1_3 = subs(Vperm);
vperm1_3=subs(vperm1_3,omega,bt1*SSn1);
vperm1_3=subs(vperm1_3,k,(Fn1_3)^(-1)*K1*Fn1_3^(-1));
vperm1_3=subs(vperm1_3,c,(Fn1_3)^(-1)*C1*Fn1_3^(-1));
vperm1_3=subs(vperm1_3,m,(Fn1_3)^(-1)*M1*Fn1_3^(-1));
 
aperm1_3 = subs(Aperm);
aperm1_3=subs(aperm1_3,omega,bt1*SSn1);
aperm1_3=subs(aperm1_3,k,(Fn1_3)^(-1)*K1*Fn1_3^(-1));
aperm1_3=subs(aperm1_3,c,(Fn1_3)^(-1)*C1*Fn1_3^(-1));
aperm1_3=subs(aperm1_3,m,(Fn1_3)^(-1)*M1*Fn1_3^(-1));
 
EQS1_3=P01_3*exp(1i*omega*tt)+M1*aperm1_3+C1*vperm1_3+K1*qperm1_3;
 
end

function [BB1_1_3,BB2_1_3,QQ1_3,VV1_3,AA1_3,ff1,ff2] = tes_1001(~,~,~)
 
format short
 
syms BB1_1_3 BB2_1_3 tt  ksi  omega  qperm qtrans AAAA  c 
syms FFn m k Q0 V0  SSn1 SSd1 Fn1_3  bt1 ksi M1
syms qst1_3 K1 C1 omega tt mpi  ag P01_3 
 
qtrans1_3=exp(-SSn1*ksi*tt).*(BB1_1_3*cos(SSd1*tt)+BB2_1_3*sin(SSd1*tt));
 
[~,qperm1_3,~] = tes_1000;
 
q11(tt)=qtrans1_3+qperm1_3;
v11(tt)=diff(q11,tt,1);
a11(tt)=diff(q11,tt,2);
 
% Condizioni al contorno
ff1 = subs(q11,tt,0)-Fn1_3*mpi*Q0;                           % q(t=0) = q0
ff2 = subs(v11,tt,0)-Fn1_3*mpi*V0;                           % v(t=0) = v0
 
% Soluzione sistema di equazioni delle condizioni al contorno
% nelle 8 costanti di integrazione incognite
Sol = solve([ff2,ff1],[BB1_1_3,BB2_1_3]);
BB1_1_3 = Sol.BB1_1_3;                          %   SOLUZIONE GENERALE
BB2_1_3 = Sol.BB2_1_3;
 
QQ1_3 = subs(q11);
QQ1_3=subs(QQ1_3,omega,bt1*SSn1);
 
VV1_3 = subs(v11);
VV1_3=subs(VV1_3,omega,bt1*SSn1);
 
AA1_3 = subs(a11);
AA1_3=subs(AA1_3,omega,bt1*SSn1);
 
end

clear, clc
close all
dt = 0.01;
tk = [0:dt:1000]; % интервал времени
q0 = [0, 0]; % начальные условия
[t Q] = ode45('dsys', tk, q0); % решение
q = Q(:,1); % берем значения q(t)
plot(t, q)
grid on
 
Y = fft(q)/length(t);
Fs = 1/dt;
f = Fs/2*linspace(-1,1,length(t));
% Plot single-sided amplitude spectrum.
figure, plot(2*pi*f, 2*abs(fftshift(Y))) 
grid on
xlim([0 7])