彈簧擺matlab演示及源代碼

1、彈簧擺寫(xiě)出彈簧擺的運(yùn)動(dòng)微分方程，并用matlab編寫(xiě)程序，演示彈簧擺的運(yùn)動(dòng)彈簧擺的運(yùn)動(dòng)微分方程可以通過(guò)拉格朗日函數(shù)法求得，對(duì)于彈簧單擺系統(tǒng)，在極坐標(biāo)系中其拉格朗日函數(shù)為其中為自由度，分別為擺球質(zhì)量，彈簧剛度，重力加速度和擺初始長(zhǎng)。將拉格朗日函數(shù)帶入到拉格朗日方程得到系統(tǒng)的運(yùn)動(dòng)微分方程：對(duì)于這個(gè)非線性方程組，利用matlab里面自帶的ode45常微分方程組求解函數(shù)，對(duì)上面的方程組進(jìn)行數(shù)值求解，得到了擺球的運(yùn)動(dòng)軌跡，如下圖所示圖中給出了彈簧擺運(yùn)動(dòng)過(guò)程中自由度隨時(shí)間的變化關(guān)系，下面將給出這個(gè)程序的源代碼，演示程序時(shí)，首先停頓三秒，然后開(kāi)始擺動(dòng)，圖形上方的彈簧擺的示意圖開(kāi)始擺動(dòng)并描繪出軌跡，同時(shí)圖形下

2、方自由度隨時(shí)間變化關(guān)系也會(huì)動(dòng)態(tài)的給出，從這個(gè)過(guò)程中可以明顯看出彈簧擺的內(nèi)共振現(xiàn)象。程序源碼：function springmasstheta0=1.5*pi/10; %單擺的初始角度m = 20/9.8; k = 80; g = 9.8; %彈簧質(zhì)量的參數(shù)保證彈簧固有頻率是單擺固有頻率的2倍l0 = 1; l = l0 + m*g/k; %l0為彈簧原來(lái)長(zhǎng)度，l為彈簧靜止時(shí)長(zhǎng)度t,u1 = ode45(weifen,0 : 0.05: 15,l0 0 theta0 0,l,k,m,g);y1,x1 = pol2cart(u1(:,3),u1(:,1); y1 = -y1;%將極坐標(biāo)換為直角坐標(biāo)

3、axis(0 16 -1.5 4);% text(0,u1(1,1),擺長(zhǎng)l);% text(0,u1(1,3),擺角theta);%axis offtitle(彈簧擺內(nèi)共振動(dòng)畫(huà)模擬及自由度隨時(shí)間的變化曲線,fontsize,15)xlabel(時(shí)間t 單位秒);ylabel(對(duì)應(yīng)的自由度);hold on;r =0.055 ; %設(shè)置彈簧半徑y(tǒng)y = -l0 : 0.01 : 0;xx = r*sin(yy./l0*30*pi);%用正弦曲線表示彈簧a,r = cart2pol(xx,yy); %用坐標(biāo)變換來(lái)畫(huà)初始位置的彈簧a = a + theta0;xx,yy = pol2cart(a,

4、r);%彈簧的數(shù)據(jù)xx=2*xx+7;yy=yy+3.5;%xx(1)=x(1)+15;line(6 8,3.5 3.5,color,k,linewidth,4)%彈簧單擺系統(tǒng)初始狀態(tài)for i=1:9 line(6+0.2*i 6+0.2*i+0.2,3.5 3.7,color,k,linewidth,1);endline(7 7,3.5 2.3,color,k,linestyle,-.,linewidth,1)bigball = line(xx(1),yy(1),color,b,marker,.,.markersize,30,erasemode,xor);%球ball2 = line(xx

5、(1),yy(1),color,g,linestyle,-,.linewidth,1.3,erasemode,none);%軌線spring = line(xx,yy,color,k,linewidth,2,erasemode, xor);%彈簧linex = line(t(1),u1(1,1),color,r,marker,.,.markersize,10,erasemode,none);%擺長(zhǎng)曲線初始liney = line(t(1),u1(1,3),color,b,marker,.,.markersize,10,erasemode,none);%擺角曲線初始%以下為標(biāo)注text(6,-0

6、.75,uparrow);text(4,-1,擺角隨時(shí)間變化曲線);text(9.0,1.9,擺長(zhǎng)隨時(shí)間變化曲線);text(11.4,1.70,downarrow);%在初始狀態(tài)下停頓三秒pause(2)%以下程序開(kāi)始實(shí)現(xiàn)動(dòng)畫(huà)，在彈簧擺開(kāi)始振動(dòng)時(shí)畫(huà)出相應(yīng)的變化曲線for i = 2 : length(t)yy = -u1(i,1) : 0.01 : 0;%彈簧xx = r*sin(yy./u1(i,1)*30*pi);%xx(i)=x(i)+15;a,r = cart2pol(xx,yy);a = a + u1(i,3);% set(linex,xdata,t(i),ydata,u1(i,1

7、);%擺長(zhǎng)初始% set(ball2,xdata,2*x1(i)+7,ydata,y1(i)+3.5);%軌線% set(liney,xdata,t(i),ydata,u1(i,3);%擺角初始plot(t(i-1),t(i),u1(i-1),1),u1(i,1),xr);plot(t(i-1),t(i),u1(i-1),3),u1(i,3),*b);%legend(擺長(zhǎng)曲線,擺角曲線)xx,yy = pol2cart(a,r);xx=2*xx+7;yy=yy+3.5;set(bigball,xdata,2*x1(i)+7,ydata,y1(i)+3.5); %球%plot(2*x1(i)+7,y1(i)+3.5,b); %球%set(ball2,xdata,2*x1(i)+7,ydata,y1(i)+3.5);%軌線plot(2*x1(i-1)+7,2*x1(i)+7,y1(i-1)+3.5,y1(i)+3.5,-g); %軌線set(spring,xdata,xx,ydata,yy);%彈簧drawnow;end%以下為彈簧擺的運(yùn)動(dòng)微分方程function valu