thsb.m

关于弹簧连接耦合摆的计算模拟。求出系统的两个简正模和一般振动解。并做出摆的位移曲线

l=15;g=9.8;m=80;k=200;h=0.6
S=[g/l+k/m,-k/m;-k/m,g/l+k/m]
P=[1,0;0,1]
[M,L]=eig(S,P)
OL=sqrt(L)
A1=[0;0.02;0.02];A2=[0.03;0;-0.03];phi1=0;phi2=0
t=0:0.04:150;
str{1}='弹簧双摆的反相振动——简正模1';
str{2}='弹簧双摆的反相振动——简正模2';
str{3}='弹簧双摆一般振动';

for j=1:3   %%三种不同的振动情况
    theta1=M(1,1)*A1(j)*cos(OL(1,1)*t+phi1)...
        +M(1,2)*A2(j)*cos(OL(2,2)*t+phi2);
    theta2=M(2,1)*A1(j)*cos(OL(1,1)*t+phi1)...
        +M(2,2)*A2(j)*cos(OL(2,2)*t+phi2);

    figure(1)
    set(gcf,'unit','normalized','position',[0.03 0.1 0.5 0.5]);
    cla;
    subplot(2,1,1)
    axis([0 20 -0.3 0.3])
    plot(t,theta1)
    xlabel('时间');
    ylabel('摆角1');
    title(str{j});
    subplot(2,1,2)
    axis([0 20 -0.3 0.3])
    plot(t,theta2)
    xlabel('时间');
    ylabel('摆角2');
    pause(0.5)
    
    figure(2)
    set(gcf,'unit','normalized','position',[0.5 0.45 0.5 0.5]);
    cla;
    axis([-15 15 -8 10])
    title(str{j});
    hold on
    
    %%画横梁和小斜线
    a10=line([-9,9],[9,9],'color','k','linestyle','-','linewidth',3.5);
    a20=linspace(-9,9,36);
    for i=1:35
        a30=(a20(i)+a20(i+1))/2
        plot([a20(i),a30],[9,9.5],'color','b',...
            'linestyle','-','linewidth',1);
    end
    %%计算两小球的直角坐标
    x1=-5+1*sin(theta1);
    y1=9-1*cos(theta1);
    x2=5+1*sin(theta2);
    y2=9-1*cos(theta2);
    
    %%用正弦函数画弹簧
    a1=linspace(x1(1),x2(1),220);
    b1=0.7*sin((a1-x1(1))*40/(x2(1)-x1(1)))+y2(1);
    tan1=line(a1,b1,'color','m','linestyle','-',...
        'erasemode','xor','linewidth',1.5);
    yuan1=line([-5,x1(1)],[9,y1(1)],'color','b','erasemode','xor',...
        'linestyle','-','linewidth',2.5);
    yuan2=line([5,x2(1)],[9,y2(1)],'color','b','erasemode','xor',...
        'linestyle','-','linewidth',2.5);
    qiu1=line(x1(1),y1(1),'color','g','erasemode','xor',...
        'marker','.','markersize',60);
    qiu2=line(x2(1),y2(1),'color','g','erasemode','xor',...
        'marker','.','markersize',60);
    
    n=length(t);
    for i=1:n
        set(yuan1,'xdata',[-5,x1(i)],'ydata',[9,y1(i)]);
        set(yuan2,'xdata',[5,x2(i)],'ydata',[9,y2(i)]);
        set(qiu1,'xdata',x1(i),'ydata',y1(i));
        set(qiu2,'xdata',x2(i),'ydata',y2(i));
        %%用正弦函数画弹簧
        a1=linspace(x1(i),x2(i),220);
        b1=0.7*sin((a1-x1(i))*40/(x2(i)-x1(i)))+y1(i);
        set(tan1,'xdata',a1,'ydata',b1);
        drawnow;
    end
end