Skip to content

Instantly share code, notes, and snippets.

@svycka

svycka/matlab.m

Last active August 29, 2015 14:07
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save svycka/5a76d6b66ea0c41b08fe to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save svycka/5a76d6b66ea0c41b08fe to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. svycka revised this gist Sep 30, 2014. 1 changed file with 5 additions and 4 deletions.
    9 changes: 5 additions & 4 deletions matlab.m
    View file
    Original file line number Diff line number Diff line change
    @@ -30,7 +30,8 @@
    % STANDUMO KOEFICIENTAI
    k=[k1, k1, k1, k1, k1];
    % KLAMPIOS TRINTIES KOEFICIENTAI
    c=[10, 10, 10, 10, 10];
    % c=[10, 10, 10, 10, 10];
    c=[ 1, m1, m1, m1, m2, m1];
    % MAZGAI (MASES), KURIUOS JUNGIA ELEMENTAI (SPYRUOKLES)
    % 3 PAREMTRAS SKIRTAS SPYRUOKLIU ATVAIZDAVIMO AUKSCIUI
    ind=[1, 3, 3;
    @@ -92,7 +93,7 @@
    %SKAITINIS INTEGRAVIMAS
    Urez=zeros(nmz,1); % rezultatu masyvas, bus pleciamas dinamiskai
    TT = 20; %Integravimo laikas
    dt=0.005; % skaitinio integravimo zingsnis
    dt=0.0005; % skaitinio integravimo zingsnis
    for t=0:dt:TT % skaitinio integravimo ciklas
    % PAGREICIU APSKAICIAVIMAS
    DDU=pagreitis(m,k,c,ind,F,tf, U,DU,t,@F_laiko_funkcija);
    @@ -169,11 +170,11 @@
    DDU= Ft(F, tf, t); % SUTEIKIAMI JEGU POVEIKIAI
    for iii=1:nel % PRIDEDAMOS ELEMENTU VEIKIAMOS JEGOS
    i=ind(iii,1);j=ind(iii,2);
    T=(U(j)-U(i))*k(iii); %+(DU(j)-DU(i))*c(iii);
    T=(U(j)-U(i))*k(iii)+(DU(j)-DU(i))*c(iii);
    DDU(i)=DDU(i)+T;
    DDU(j)=DDU(j)-T;
    end
    a = 0.01;
    a = 3;
    for iii=1:nmz
    DDU(i) = DDU(i)+ a * DU(iii) * m(i);
    end
  2. svycka created this gist Sep 30, 2014.
    208 changes: 208 additions & 0 deletions matlab.m
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,208 @@
    function pvz_2
    clc, close all, clear all

    %KONSTRUKCIJOS DUOMENYS
    % MAZGU (MASIU) KONSTANTOS
    m1 = 1; m2 = 3;
    % ELEMENTU (SPYRUOKLIU) KONSTANTOS
    k1 = 2000;
    % VEIKIANCIU JEGOS (N)
    f1=800; f2=400;
    % JEGU VEIKIMO PRADZIA (s)
    tf1 = 0.2; tf2 = 0.1;
    % MAZGU ITVIRTINIMO MASYVAS
    IS = [ 1, 0, 0, 0, 0, 0];
    % 2 MAZG0 KINEMATINIIS POVEIKIIS (m)
    % u2_deltaU = 0.5;
    % % 2 MAZGO KINEMATINIO POVEIKIO PRADZIA (s)
    % u2_t_start = 0.2;
    % % 2 MAZGO KINEMATINIO POVEIKIO TRUKME (s)
    % u2_deltaT = 0.1;
    % MAZGAI (MASES) - MATRICU PAVIDALU
    % MAZGU MASES (kg)
    m= [ 1, m1, m1, m1, m2, m1];
    % MAZGUS VEIKIAN?IOS J?GOS (N)
    F = [ 0, -f2, 0, 0, 0, -f1];
    % JEGU POVEIKIO PRADZIA (s)
    tf = [ 0, tf2, 0, 0, 0, tf1];

    % ELEMENTAI (SPYRUOKLES)
    % STANDUMO KOEFICIENTAI
    k=[k1, k1, k1, k1, k1];
    % KLAMPIOS TRINTIES KOEFICIENTAI
    c=[10, 10, 10, 10, 10];
    % MAZGAI (MASES), KURIUOS JUNGIA ELEMENTAI (SPYRUOKLES)
    % 3 PAREMTRAS SKIRTAS SPYRUOKLIU ATVAIZDAVIMO AUKSCIUI
    ind=[1, 3, 3;
    2, 4, 1;
    3, 5, 3;
    4, 5, 1;
    5, 6, 2];
    % ATVAIZDAVIMO DUOMENYS
    % MAZGU X KORDINATES
    x=[2, 2, 4, 4, 6, 8];
    % Mazgu Y KOORDINATES
    y=[3, 1, 3, 1, 2, 2];
    % MASIU DIAMETRAI
    rad1=0.2; rad2=1.5;
    rads = [rad1, rad1;
    rad1, rad1;
    rad1, rad1;
    rad1, rad1;
    rad1, rad2;
    rad1, rad1];

    nmz=length(m); % MAZGU SKAICIUS (PAVYZDYJE = 6)
    nel=length(k); % ELEMENTU SKAICIUS (PAVYZDYJE = 5)
    %PRADINIAI GREICIAI IR POSLINKIAI LYGUS 0
    U=zeros(nmz,1); % PRADINIAI POSLINKIAI = 0
    DU=zeros(nmz,1); % PRADINIAI GREICIAI = 0
    % U(2) = u2_deltaU;
    % atvaizduojami pradiniai duomenys
    % vaizdavimas(x, y, ind, U, rads, 0);

    K=zeros(nmz,nmz);
    for iii=1:nel
    i = ind(iii,1); j = ind(iii,2);
    Kel = [k1, -k1;
    -k1, k1];
    K([i,j],[i,j]) = K([i,j],[i,j])+Kel;
    end
    F = [ 0, 0, 0, 0, 0, -f1];
    K00 = K(find(~IS), find(~IS));
    U0 = U(find(~IS));
    F0 = F(find(~IS));
    K01 = K(find(~IS), find(IS));
    U1 = U(find(IS));

    % nezinomi poslinkiai
    U0 = K00 \ (F0' - K01*U1);
    U(find(~IS)) = U0;
    U(find(IS)) = U1;
    'poslinkiai atlikus skaitin analize'
    disp(U');
    vaizdavimas(x, y, ind, U, rads, 1);
    F = [ 0, -f2, 0, 0, 0, -f1];
    U = zeros(nmz,1);





    %SKAITINIS INTEGRAVIMAS
    Urez=zeros(nmz,1); % rezultatu masyvas, bus pleciamas dinamiskai
    TT = 20; %Integravimo laikas
    dt=0.005; % skaitinio integravimo zingsnis
    for t=0:dt:TT % skaitinio integravimo ciklas
    % PAGREICIU APSKAICIAVIMAS
    DDU=pagreitis(m,k,c,ind,F,tf, U,DU,t,@F_laiko_funkcija);
    % GREICIU ATNAUJINIMAS
    DU=DU+dt*DDU;
    % PRISKIRIAMOS GREICIU KRASTINES SALYGOS
    DU(find(IS))=0;
    % DU(2) = du_laiko_funkcija(u2_t_start, u2_deltaT, u2_deltaU, t);
    % POSLINKIU ATNAUJINIMAS
    U=U+dt*DU;
    % ISSAUGOMI TARPINIAI REZULTATAI
    jj=round(t/dt)+1;
    % POSLINKIU REIKSMIU ISSAUGOJIMAS
    Urez(:,jj)=U;
    % ATVAIZDUOJAMA KONSTRUKCIJOS BUSENA
    % figure(1); vaizdavimas(x, y, ind, U, rads,t);
    % pause(0.01);
    end

    'poslinkiai atlikus i?reik?tin? analize'
    disp(U');

    % ATVAIZDUOJAMAS POSLINKIU NUO LAIKO GRAFIKAS
    figure(2);hold on; grid on;
    axis([0 ,TT,-2 ,0.5]);
    spalva={'b-';'r-';'g-';'m-';'c-';'k-';};
    for i=1:nmz
    plot([0:dt:TT],Urez(i,:),spalva{i});
    end
    % plot([0.05, 0.05], [-0.2, 0.2], 'g--');
    % plot([0.1, 0.1], [-0.2, 0.2], 'k--');
    % plot([0.2, 0.2], [-0.2, 0.2], 'r--');
    legend('1 mazgas','2 mazgas', '3 mazgas', '4 mazgas', '5 mazgas', '6 mazgas');
    xlabel('Laikas (s)');
    ylabel('Poslinkiai (m)');
    end

    % JEGU LAIKO FUNKCIJA
    function Ft=F_laiko_funkcija(F, tf, t)
    nmz=length(F);
    Ft = zeros(nmz, 1);
    for i=1:nmz,
    if(tf(i) < t)
    Ft(i) = F(i);
    end
    end
    if t>0.1 && t<0.2
    Ft(2) = F(2);
    else
    Ft(2)=0;
    end
    end

    % POSLINKU LAIKO FUNKCIJA
    function du=du_laiko_funkcija(t_start, deltaT, deltaU, t)
    if t >= t_start % JEI POSLINKIS JAU PRASIDEJO
    if t <= t_start + deltaT, %JEI POSLINKIS ATLIEKAMAS
    % POSLINKIS VYKSTA, APSKAICIUOJAMAS MAZGO GREITIS
    omega=(pi/2)/deltaT;
    du = deltaU * omega * cos(omega *(t - t_start));
    else,
    du = 0; % POSLINKIS JAU IVYKO, MAZGO GREITIS = 0
    end
    else
    du = 0; % POSLINKIS DAR NEPRASIDEJO, MAZGO GREITIS = 0
    end
    return
    end

    % *********************************************
    function DDU=pagreitis(m,k,c,ind,F,tf,U,DU,t,Ft)
    nel=length(k);nmz=length(m);
    DDU=zeros(nmz,1);
    DDU= Ft(F, tf, t); % SUTEIKIAMI JEGU POVEIKIAI
    for iii=1:nel % PRIDEDAMOS ELEMENTU VEIKIAMOS JEGOS
    i=ind(iii,1);j=ind(iii,2);
    T=(U(j)-U(i))*k(iii); %+(DU(j)-DU(i))*c(iii);
    DDU(i)=DDU(i)+T;
    DDU(j)=DDU(j)-T;
    end
    a = 0.01;
    for iii=1:nmz
    DDU(i) = DDU(i)+ a * DU(iii) * m(i);
    end
    DDU=DDU./m';
    return
    end

    % *********************************************
    function vaizdavimas(x,y,ind,U,rads,t)
    %PARUOSIAMA FIGURA ATVAIZDAVIMUI
    clf; hold on; grid on;axis equal;
    axis([min(x)-1 ,max(x)+1,min(y)-1 ,max(y)+1]);
    title(['t =', num2str(t), '(s)']);
    xlabel('X');
    ylabel('Y');
    %MAZGU IR ELEMENTU SKAICIUS
    nmz=length(x);nel=size(ind,1);
    % ATVAIZDUOJAMI ELEMENTAI
    for i=1:nel
    line([x(ind(i,1))+U(ind(i,1)), ...
    x(ind(i,2))+U(ind(i,2)) ], ...
    [ind(i,3), ind(i,3)],...
    'Color','blue','Linewidth',2);
    end
    %ATVAIZDUOJAMI MAZGAI
    for i=1:nmz
    rectangle('Position',[x(i)+U(i)-rads(i,1), ...
    y(i)-rads(i,2),rads(i,1)*2,rads(i,2)*2]....
    ,'Curvature',[1,1],'FaceColor',[0.4 0.6 1]);
    end
    return
    end