Kulawik Bartosz
Gr. 23 IIID

Analiza i Identyfikacja Sygnałów
Sprawozdanie nr 6
Elementy Algebry Liniowej

1. Ćwiczenie 1
   Dla układu jak na rysunku dobierz tak parametry macierzy M, C i K aby pierwsza częstość drgań własnych znajdowała się nieco poniżej 5 Hz, a druga powyżej 50 Hz.
   n=2;
   m1=14;
   m2=1;
   al1 = 9;
   al2 = 3;
   al3 = 4;
   k1 = 10000;
   k2 = 84000;
   k3 = 10000;
   % Mas
   M = [m1, 0;
   0, m2];
   % Ws. tłumienia
   C = [al1+al2, -al2;
   -al2, al2+al3];
   % Ws. sztywności
   K = [k1+k2, -k2;
   -k2, k2+k3];
   ZER = zeros(size(M));
   A = [ZER,M;M,C];
   B = [-M,ZER;ZER,K];
   % Rozwiązanie uogólnionego zagadnienia wlasnego
   [PHI,LAMBDA]=eig(-B,A);
   WD=imag(diag(LAMBDA))/2/pi;
   % Czestotliwosci drgan wlasnych [Hz]
   WW=sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2)/2/pi
   % Tlumienie modalne
   KSI=-real(diag(LAMBDA))./sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2);
   WW =
   
   50.1869
   50.1869
   5.6911
   5.6911

2. Ćwiczenie 2
   Dla dobranych w poprzednim zadaniu parametrów dokonaj syntezy widmowych funkcji przejścia w oparciu o wzór (9) oraz o wzór:
   
   % Liczba stopni swobody
   n=2;
   % Masy w układzie
   m1=14;
   m2=1;
   % Współczynniki tłumienia
   al1 = 9;
   al2 = 3;
   al3 = 4;
   % Współczynniki sztywności
   k1 = 10000;
   k2 = 84000;
   k3 = 10000;
   % Macierze współczynników
   % Mas
   M = [m1, 0;
   0, m2];
   % Ws. tłumienia
   C = [al1+al2, -al2;
   -al2, al2+al3];
   % Ws. sztywności
   K = [k1+k2, -k2;
   -k2, k2+k3];
   % Budowanie macierzy do równań stanu w oparciu o wzor (8)
   ZER = zeros(size(M));
   A = [ZER,M;M,C];
   B = [-M,ZER;ZER,K];
   % Rozwiązanie uogólnionego zagadnienia wlasnego
   [PHI,LAMBDA]=eig(-B,A);
   % Czestotliwosci drgan wlasnych tlumionych [Hz]
   WD=imag(diag(LAMBDA))/2/pi;
   % Czestotliwosci drgan wlasnych [Hz]
   WW=sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2)/2/pi
   % Tlumienie modalne
   KSI=-real(diag(LAMBDA))./sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2);
   % Estymacja współczynników skalujących
   AN=zeros(4);
   AAA=PHI'*A*PHI;
   for a=1:n
   AN(:,2*a)=AAA(:,2*a-1);
   AN(:,2*a-1)=AAA(:,2*a);
   end;
   QQ=inv(AN);
   Q=diag(QQ);
   % Synteza WFP zgodnie ze wzorem (9)
   for c=1:2
   jj=[1:2];
   f=[0:0.25:60];
   for b=1:length(f)
   for a=1:n
   htemp=0;
   for r=1:2*n
   htemp=htemp+((Q(r)*PHI(n+a,r)* PHI(n+jj(c),r))/(i*f(b)*2*pi- LAMBDA(r,r)));
   end;
   H{a,c}(b)=htemp*(-1)*(f(b)*2*pi)^2;
   end;
   end;
   end;
   plot(f,abs(H{1,1}))
   title('Widmowa funkcja przejscia 1masy')
   xlabel('Czestotliwosc')
   ylabel('amplituda')
   grid on
   figure
   plot(f,abs(H{1,2}),'r')
   title('Widmowa funkcja przejscia 2masy')
   xlabel('Czestotliwosc')
   ylabel('amplituda')
   grid on
   
   A=[ 0 1 0 0;
   -(k1+k2)/m1 -(al1+al2)/m1 k2/m1 al2/m1;
   0 0 0 1;
   k2/m2 al2/m2 -(k3+k2)/m2 -(al3+al2)/m2 ];
   B=[0;
   1/m1;
   0;
   1/m2 ];
   C=[ 1 0 0 0;
   0 0 1 0 ];
   D=[ 0; 0 ];
   [l, m]=ss2tf(A,B,C,D);
   l1=l(1,:);
   l2=l(2,:);
   s=j*2*pi*f;
   X1=(l1(1,1)*s.^4+l1(1,2)*s.^3+l1(1,3)*s.^2+l1(1,4)*s+l1(1,5))./(m(1,1)*s.^4+m(1,2)*s.^3+m(1,3)*s.^2+m(1,4)*s+m(1,5));
   X2=(l2(1,1)*s.^4+l2(1,2)*s.^3+l2(1,3)*s.^2+l2(1,4)*s+l2(1,5))./(m(1,1)*s.^4+m(1,2)*s.^3+m(1,3)*s.^2+m(1,4)*s+m(1,5));
   figure
   plot(f,abs(X1*(1/max(X1))))
   title('Masa nr1')
   xlabel('Czestotliwosc')
   ylabel('amplituda')
   grid on
   figure
   plot(f,abs(X2*(1/max(X2))))
   title('Masa nr2')
   xlabel('Czestotliwosc')
   ylabel('amplituda')
   grid on
   
   
   
   
   
   
   
   
   
   
   
   Otrzymane wykresy:
   Metoda 1
   
   
   
   
   
   
   
   
   Metoda 2
   
   
   

3. Ćwiczenie 3
   Dla układu z zadania 1 przeprowadź analizę modalną zastępując rozkład na wartości własne rozkładem na wartości szczególne.
   % Liczba stopni swobody
   n=2;
   % Masy w układzie
   m1=14;
   m2=1;
   % Współczynniki tłumienia
   al1 = 9;
   al2 = 3;
   al3 = 4;
   % Współczynniki sztywności
   k1 = 10000;
   k2 = 84000;
   k3 = 10000;
   % Macierze współczynników
   % Mas
   M = [m1, 0;
   0, m2];
   % Ws. tłumienia
   C = [al1+al2, -al2;
   -al2, al2+al3];
   % Ws. sztywności
   K = [k1+k2, -k2;
   -k2, k2+k3];
   % Budowanie macierzy do równań stanu w oparciu o wzor (8)
   ZER = zeros(size(M));
   A = [ZER,M;M,C];
   B = [-M,ZER;ZER,K];
   % Rozwiązanie uogólnionego zagadnienia wlasnego
   [PHI,LAMBDA]=eig(-B,A);
   % Czestotliwosci drgan wlasnych tlumionych [Hz]
   WD=imag(diag(LAMBDA))/2/pi;
   % Czestotliwosci drgan wlasnych [Hz]
   WW=sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2)/2/pi
   % Tlumienie modalne
   KSI=-real(diag(LAMBDA))./sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2);
   % Estymacja współczynników skalujących
   AN=zeros(4);
   AAA=PHI'*A*PHI;
   for a=1:n
   AN(:,2*a)=AAA(:,2*a-1);
   AN(:,2*a-1)=AAA(:,2*a);
   end;
   QQ=inv(AN);
   Q=diag(QQ);
   % Synteza WFP zgodnie ze wzorem (9)
   for c=1:2
   jj=[1:2];
   f=[0:0.25:60];
   for b=1:length(f)
   for a=1:n
   htemp=0;
   for r=1:2*n
   htemp=htemp+((Q(r)*PHI(n+a,r)* PHI(n+jj(c),r))/(i*f(b)*2*pi- LAMBDA(r,r)));
   end;
   H{a,c}(b)=htemp*(-1)*(f(b)*2*pi)^2;
   end;
   end;
   end;
   plot(f,abs(H{1,1}))
   hold on
   plot(f,abs(H{1,2}),'r')
   grid on
   
   % definicja sygnału wymuszenia
   t=[0:1/80:10];
   % definicja wektora siły o charakterze losowym
   f_zadane_1=10*rand(size(t));
   % rozdzielczosc czestotliwosciowa
   nn=round(120/0.25);
   % transformata Fouriera
   f_zadane_1_freq=(fft(f_zadane_1,nn)*2)/nn;
   dl=length(f_zadane_1_freq);
   f_zadane_freq{1}=f_zadane_1_freq(1:ceil(dl/2)+1);
   f_zadane_freq{1}(1:5)=0;
   % wyliczenie odpowiedzi na wymuszenie w masie nr 1
   frfs=H(:,1);
   for a=1:length(f)
   tfrfs=[];
   tforces=[];
   for b=1:length(f_zadane_freq)
   tforces(b)=f_zadane_freq{b}(a);
   end;
   for c=1:size(frfs,2)
   for b=1:size(frfs,1)
   tfrfs(b,c)=frfs{b,c}(a);
   end;
   end;
   tresps=tfrfs*tforces';
   for b=1:size(frfs,1)
   resp{b}(a)=tresps(b);
   end;
   end;
   % wyliczenie wektora siły wymuszającej f
   for a=1:length(resp{1})
   
   % tworzenie tymczasowego wektora widm odpowiedzi i WFP dla kolejnych
   % czestotliwosci
   tresps=[];
   for b=1:length(resp)
   tresps(b)=resp{b}(a);
   end;
   for c=1:size(frfs,2)
   for b=1:size(frfs,1)
   tfrfs(b,c)=frfs{b,c}(a);
   end;
   end;
   tresps=tresps.';
   % rozwiazanie rownania f=(H^-1)*p metoda rozkladu na wartosci szczegolne
   [u,sig,v]=svd(tfrfs);
   r=rank(tfrfs);
   beta=u'*tresps;
   for b=1:r
   en(b)=beta(b)/sig(b,b);
   end;
   for b=r+1:size(v,1)
   e2(b)=0;
   end;
   if r>0
   esizen=size(en,2);
   esize2=size(e2,2);
   en=[en,zeros(1,size(v,1)-esizen)];
   e2=[zeros(1,size(v,1)-esize2),e2];
   x1=v*en';
   x2=v*e2';
   ftemp=x1+x2;
   else
   esize2=size(e2,2);
   e2=[zeros(1,size(v,1)-esize2),e2];
   x2=v*e2';
   ftemp=x2;
   end;
   for b=1:length(ftemp)
   Identforce{b}(a)=ftemp(b);
   end;
   end;
   figure
   plot(f,abs(Identforce{1}))
   title('Widmo wymuszenia')
   xlabel('czestotliwosc')
   ylabel('amplituda')
   grid on
   figure
   plot(f,abs(resp{1}))
   title('Widmo odpowiedzi ukladu dla 1masy')
   xlabel('czestotliwosc')
   ylabel('amplituda')
   grid on
   figure
   plot(f,abs(resp{2}))
   title('Widmo odpowiedzi ukladu dla 2masy')
   xlabel('czestotliwosc')
   ylabel('amplituda')
   grid on
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   Otrzymane wykresy:
   
   
   
   W przypadku rozkładu macierzy względem wartości szczególnych obserwujemy wzrost amplitudy dla wartości częstotliwości około 5Hz oraz 50Hz. Natomiast w odróżnieniu od sposobu rozkładu macierzy względem wartości własnych przebieg widma jest „poszarpany”.

4. Ćwiczenie 4
   Rozwiąż układ równań (21), gdzie A będzie macierzą losowych współczynników o wymiarach 1000 x 1000, a wektor y to wektor losowych wyrazów wolnych o rozmiarze 1000 x 1. Zastosuj metodę klasyczną i metodę LU. Porównaj czasy obliczeń.
   
   W=rand(1000,1000);
   Y=rand(1000,1);
   tic
   X=W\Y;
   toc
   tic
   [L,U]=lu(W);
   z=L\Y;
   x=U\z;
   toc
   Czas obliczania metodą klasyczną: Elapsed time is 0.892621 seconds.
   Czas obliczania metodą LU: Elapsed time is 0.638388 seconds.