Вывод точечного графика тел в пространстве

clear all
m1= m0*(1+rand(1))
m2=m0*(1+0.8*rand(1))
m3=m0*(1+0.6*rand(1))
 
k1= k0*(1+rand(1))
k2= k0*(1+0.8*rand(1))
k3=m0*(1+0.6*rand(1))
k4=m0*(1+0.4*rand(1))
 
>> N=3; % число тел колебательной системы
>> m=[m1 m2 m3]; % массы тел колебательной системы
>> k=[k1 k2 k3 k4]; % жесткости пружин колебательной системы
>> R0=[-0.2 0 -0.3]; % смещения тел в момент времени t = 0
>> v0=[1 -3 0]; % скорости тел в момент времени t = 0
% вычисление элементов матрицы 
>> for alpha=1:N+1
        for beta=1:N
          omega(alpha,beta)=k(alpha)/m(beta);
        end;
     end;
>> i=1:N;
>> j=1:N;
% вычисление элементов матрицы OMEGA в соответствии с (8)
>> OMEGA(i,j)=0;
>> for i=1:N
       if i==1
        OMEGA(i,i)=omega(1,1)+omega(2,1);
        OMEGA(1,2)=-omega(2,1);
      end;
      if i>1
        if i<N
            OMEGA(i,i-1)=-omega(i,i);
            OMEGA(i,i)=omega(i,i)+omega(i+1,i);
            OMEGA(i,i+1)=-omega(i+1,i);
        else
            OMEGA(i,i-1)=-omega(i,i);
            OMEGA(i,i)=omega(i,i)+omega(i+1,i);
        end;
      end;
    end;
>> [Sigma,Teta]=eig(OMEGA); % вычисление собственных значений и 
% собственных векторов матрицы OMEGA
>> Teta=Teta^0.5; % вычисление собственных частот
>> for i=1:N
        for j=1:N
          SigmaV(j,i)=-Teta(i,i)*Sigma(j,i);
        end;
    end;
>> C1=Sigma^-1*R0'; % решение системы уравнений (18)
>> C2=SigmaV^-1*v0'; % решение системы уравнений (19)
>> C=(C1.^2+C2.^2).^0.5; % вычисление координат вектора C
>> clear alpha
% вычисление фазы нормальных колебаний в соответствие с (21), (22)
>> for i=1:N
        if C(i)==0 
          alpha(i)=0;
        else
          alpha(i)=atan(C2(i)./C1(i));
          if C1(i)<0
            alpha(i)=pi+alpha(i);
          end;
          if C1(i)>0
            if C2(i)<0
              alpha(i)=2*pi+alpha(i);
            end;
          end;
        end;
    end;
>> N=length(OMEGA);
>> N1=2^13; % число узлов временной сетки
>> j=1:N1;
>> Tmax=80; % правая граница временного интервала
>> t(j)=(j-1)/(N1-1)*Tmax; % координаты узлов временной сетки
% вычисление значений координат тел в узлах временной сетки
>> for j=1:N1
        s=zeros(N,1);
        for i=1:N
          s=s+C(i)*Sigma(:,i).*cos(Teta(i,i)*t(j)+alpha(i));
        end;
        X(:,j)=s;
    end;
% вычисление значений скоростей тел в узлах временной сетки
>> for j=1:N1
        s=zeros(N,1);
        for i=1:N
          s=s+C(i)*Sigma(:,i).*Teta(i,i)*sin(Teta(i,i)*t(j)+alpha(i));
        end;
        Xv(:,j)=-s;
    end;
% визуализация зависимостей мгновенных значений смещений и 
% скорости от времени
>> figure(1);figure(1);plot(t,X(1,:),'-k',t,X(2,:),'--k',t,X(3,:),':k')
>> figure(2);plot(t,Xv(1,:),'-k',t,Xv(2,:),'--k',t,Xv(3,:),':k')
% построение траектории движения тел на фазовой плоскости
>> figure(3);plot(X(1,:),Xv(1,:))
>> figure(4);plot(X(2,:),Xv(2,:))
>> figure(5);plot(X(3,:),Xv(3,:))
% вычисление спектров зависимостей смещений тел системы от времени
>> c1=fft(X(1,:));
>> c2=fft(X(2,:));
>> c3=fft(X(3,:));
>> j=2:N1/2;
% вычисление спектральной плотности смещений тел
>> Cm1(j-1)=abs(c1(j-1))/(N1/2);
>> Cm2(j-1)=abs(c2(j-1))/(N1/2);
>> Cm3(j-1)=abs(c3(j-1))/(N1/2);
>> Freq(j-1)=(j-1)/Tmax; % вычисление частот спектральных гармоник
% визуализация спектральных плотностей смещений тел
>> figure(6);semilogy(Freq,Cm1,'-k',Freq,10*Cm2,'--k',Freq,500*Cm3,':k')
>> axis([0 2.5 10^-3 2000])