单自由度系统的响应图，相位图，根轨迹图绘制

博大广阔

function free_vibration_with_vious_damping(m,k,c,x0,dx0,t)
%equation :                              M*ddx+C*dx+k*x=0
%again asuming :                         xt=c*exp(s*t),                    
%and then inserting this function into equation
%obtain the characteristic equation :    m*s^2+c*s+k=0  (the roots is s1 and s2
%these roots give two solution:          xt1=C1*exp(s1*t) and xt2=C2*exp(s2*t)
%thus the general solution is :          x(t)=C1*exp(s1*t)+C2*exp(s2*t)  
% where C1and C2 is determined by the initial conditions of the system
%howerer S12=(-c+/-sqrt(C^2-4mk))/2M;   (2)
%the radical in the equation(2)becomes zero as Critial Damping Cc=2*sqrt(k*m)
%damping radio e=C/Cc;
%then wei can rewrite the equation of (2) as :
%                         S12=(-e +/- sqrt(e^2-1))*wn=-e*wn +/- wd (wd=sqrt(1-e^2)*wn)
%now wei must discuss three situation

%Frist one named underddamped system (0<e<1)
%      S1=(-e+i*sqrt(1-e^2))*wn     and S2=(-e - i*sqrt(1-e^2))*wn
%at last X(t)=expt(-e*wn*t)*{ C1*cos(wd*t)+C2*sin(wd*t)}  
% where   C1=x0;C2=(dx0+e*wn*x0)/wd

%the second one e=1
% S1=S2=-wn
%because of the repeated rotos,the solution of equation of(1) is
%X(t)=(C1+C2*t)*exp(-wn*t)
% where C1=x0 ;C2=dx0+wn*x0

%the thrid situation overdamped system (e>1)
% S1=(-e+sqrt(e^2-1))*wn  S2=(-e-sqrt(e^2-1))*wn
%X(t)=C1*exp(S1*t)+C1*exp(S2*t)
%where C1=(x0*wn(e+sqrt(e^2-1))+dx0)/(2*wn*sqrt(e^2-1))
%      C2= (-x0*wn(e-sqrt(e^2-1))-dx0)/(2*wn*sqrt(e^2-1))
clc
clear all
%M=450;K=26519.2;C=1000; x0=0.539;dx0=1;t=0:0.001:10; %checking!!
M=4;K=4;C=20;x0=0.539;dx0=1;t=0:0.001:10;
[m,n]=size(C);
%if nargin<6 t=0:0.01:20;end
if m==n==1
wn=sqrt(K/M); e=C/(2*sqrt(K*M));             %damping radio

if e<1&e>0
    wd=sqrt(1-e^2)*wn;
    C1=x0;C2=(dx0+e*wn*x0)/wd;
Xt=exp(-e*wn.*t).*( C1*cos(wd.*t)+C2*sin(wd.*t));
elseif e==1
    C1=x0 ;C2=dx0+wn*x0;
    Xt=(C1*ones(size(t))+C2.*t).*exp(-wn.*t);
elseif e>1
    S1=(-e+sqrt(e^2-1))*wn ; S2=(-e-sqrt(e^2-1))*wn;
    C1=(x0*wn*(e+sqrt(e^2-1))+dx0)/(2*wn*sqrt(e^2-1)) ;
    C2= (-x0*wn*(e-sqrt(e^2-1))-dx0)/(2*wn*sqrt(e^2-1));
    Xt=C1.*exp(S1.*t)+C1.*exp(S2.*t);
end
figure(1)
plot(t,Xt,'r'); grid on; xlabel('time');ylabel('diaplacement');title('free vibration with vious damping');
end
%plot the Locus of S1 and S2 of this system
   wn=sqrt(K/M); Cc=2*sqrt(K*M);
    e1=[0:0.01:1];    e2=[1:0.1:1.1];
    S1=(-e1+i*sqrt(1-e1.^2)).*wn;  S2=(-e1 - i*sqrt(1-e1.^2)).*wn;
    S3=(-e2+sqrt(e2.^2-1)).*wn ;      S4=(-e2-sqrt(e2.^2-1)).*wn;
    S1r=real(S1);S1i=imag(S1);S2r=real(S2);S2i=imag(S2);
    S3r=real(S3);S3i=imag(S3);S4r=real(S4);S4i=imag(S4);
    figure(2)  
    plot(S1r,S1i,'r');grid on; hold on;
    plot(S2r,S2i,'b');title('locus of this systen');hold on
     plot(S3r,S3i,'*'); plot(S4r,S4i,'-');  xlabel('real axis');
     %plot Phane plane of damped system 使用R-K方法
     
     x1=x0;x2=dx0;T=0;h=0.001;n=1;x1old=x1;x2old=x2;
     while T<=10*pi/wn;         
         K11=x2old;
         K12=-(C*x2old+K*x1old)/M;
         T=T+h/2; x1=x1old+ K11*h/2;x2=x2old+0.5*h*K12;                 
         K21=x2;
         K22=-(C*x2+K*x1)/M;
         T=T+h/2;x1=x1old-h*K11+2*h*K21;x2=x2old-h*K12+2*h*K22;
         K31=x2;
         K32=-(C*x2+K*x1)/M;
         x1old=x1old+(h/6)*(K11+4*K21+K31);
         x2old=x2old+(h/6)*(K12+4*K22+K32);
         x11(1,n)=x1old;
         x21(1,n)=x2old;
         n=n+1;
        
     end
   
   
         figure(3)
   plot(x11,x21,'r');
grid on
xlabel('displacement')
ylabel('speed')
title('phase plane')
     
end
   
   





   
   
   

ChaChing

原创 or 转贴 !? 个人专业有限, 看似像分享!? 是吗?
若是, 若能稍作简易说明下相关输入/输出/算法, 甚至举例, 或许更适合学习! 个人浅见