gily 发表于 2006-12-8 09:38

运算双摆动力学时出现的问题,请指点迷津

我写了个计算一个普通双摆的动力学的代码,两杆的初始角度a1=30(degree),a2=45(degree),从静止开始摆动。通过龙格库塔法求解动力学算出两杆中心的位置,但发现杆的摆动越来越大,不知是什么原因。动力学方程应该没错误,并且已经进行了位置和速度的违约修正,误差已经相当小了。但输出的曲线说明摆的幅度越来越大,如图红线为摆1的中心点在X轴的运动轨迹,绿线为杆二中心点在X轴的运动轨迹。
请问这是怎么回事呢?

pengweicai 发表于 2006-12-8 12:34

铰链连接的双摆

程序来自 《理论力学计算机模拟(Matlab实现)》书中光盘内容。

l=9;
   =ode45('jjsbfun',,,[],l);
   y1=-l*cos(u(:,1));%\fs{计算球1的坐标}
   x1=l*sin(u(:,1));
   y2=y1-l*cos(u(:,3));%\fs{计算球2的坐标}
   x2=x1+l*sin(u(:,3));
   figure(1)
   set(gcf,'unit','normalized','position',);%\fs{设置窗口坐标}
   cla;
   plot(t,u(:,1),'g',t,u(:,3),'r')%\fs{画位移曲线}
   xlabel('时间');ylabel('摆角');
   legend('上面的摆','下面的摆');
   pause(0.5)
   
   figure(2)
   set(gcf,'unit','normalized','position',);%\fs{设置窗口坐标}
   cla;
   axis([-10 10 -20 20])
   axis equal
   hold on

a10=line([-9,9],,'color','k','linewidth',3.5);%\fs{画横梁}
%\fs{下面的循环语句是画横梁顶部的虚线}
a20=linspace(-9,9,36);   
for i=1:35
   a30=(a20(i)+a20(i+1))/2;
   plot(,,'color','b','linestyle','-','linewidth',1);
end


ball1a=line(x1(1),y1(1),'color',,'linestyle','-','linewidth',1,...
   'erasemode','none');
ball1=line(x1(1),y1(1),'color','r','marker','.','markersize',40,'erasemode','xor');
ball2a=line(x2(1),y2(1),'color',,'linestyle','-','linewidth',1,...
   'erasemode','none');
ball2=line(x2(1),y2(1),'color','r','marker','.','markersize',40,'erasemode','xor');
gan1=line(,,'color','g','linewidth',2,'erasemode','xor');
gan2=line(,,'color','g','linewidth',2,'erasemode','xor');

for i=1:length(u)
   set(ball1,'xdata',x1(i),'ydata',y1(i));
   set(ball2,'xdata',x2(i),'ydata',y2(i));
   set(ball1a,'xdata',x1(i),'ydata',y1(i));
   set(ball2a,'xdata',x2(i),'ydata',y2(i));
   set(gan1,'xdata',,'ydata',);
   set(gan2,'xdata',,'ydata',);
   drawnow
end


你可以参考一下。

gily 发表于 2006-12-8 13:53

非常感谢pengweicai的热心帮助,但上面这段代码的主要语句即 =ode45('jjsbfun',,,[],l); 而我用的动力学方法是带拉格朗日乘子的动力学方程,所以求解的是微分代数方程组。但这两种方法结果应该差不多呀。我现在奇怪的是我的结果怎么会发散呢?已经进行了违约修正啊?

xjzuo 发表于 2006-12-8 19:18

回复

建议先把代码传上来.

gily 发表于 2006-12-11 09:00

下面是双摆动力学在Maple下的代码,求解微分代数方程,先把它转化为线性代数方程得到加速度,再求解常微分方程,得到位置和速度。A、B为广义质量阵和广义力,以下按时间步长循环求解下一步的位置和速度。

循环求解每一时间步的位置和速度
> for compute_circle from 1 to circle do      
>phi1:=mid3; phi2:=mid6; phi11:=mid9; phi22:=mid12;e_x1:=mid1;e_y1:=mid2;e_x2:=mid4;e_y2:=mid5; e_x11:=mid7; e_y11:=mid8; e_x22:=mid10; e_y22:=mid11;

> A:=Matrix([,,,,,,,,,]);

> B:=Vector([0,-m*g,0,0,-m*g,0,
-line_L0*sin(phi1)*phi11^2-2*alpha*(e_x11-line_L0*cos(phi1)*phi11)-beta^2*(e_x1-line_L0*sin(phi1)),
line_L0*cos(phi1)*phi11^2-2*alpha*(e_y11-line_L0*sin(phi1)*phi11)-beta^2*(e_y1+line_L0*cos(phi1)),
-2*line_L0*sin(phi1)*phi11^2-line_L0*sin(phi2)*phi22^2-2*alpha*(e_x22-2*line_L0*cos(phi1)*phi11-line_L0*cos(phi2)*phi22)-beta^2*(e_x2-2*line_L0*sin(phi1)-line_L0*sin(phi2)),
2*line_L0*cos(phi1)*phi11^2+line_L0*cos(phi2)*phi22^2-2*alpha*(e_y22-2*line_L0*sin(phi1)*phi11-line_L0*sin(phi2)*phi22)-beta^2*(e_y2+2*line_L0*cos(phi1)+line_L0*cos(phi2)) ]);

> linearsolution:=LinearSolve(A,B);

解ODE
> eqn1:=diff(y1(t),t)=z7(t);
> eqn2:=diff(y2(t),t)=z8(t);
> eqn3:=diff(y3(t),t)=z9(t);
> eqn4:=diff(y4(t),t)=z10(t);
> eqn5:=diff(y5(t),t)=z11(t);
> eqn6:=diff(y6(t),t)=z12(t);
> eqn7:=diff(z7(t),t)=linearsolution;
> eqn8:=diff(z8(t),t)=linearsolution;
> eqn9:=diff(z9(t),t)=linearsolution;
> eqn10:=diff(z10(t),t)=linearsolution;
> eqn11:=diff(z11(t),t)=linearsolution;
> eqn12:=diff(z12(t),t)=linearsolution;

初值
> ICs:=y1(0)=mid1,y2(0)=mid2,y3(0)=mid3,y4(0)=mid4,y5(0)=mid5,y6(0)=mid6,z7(0)=mid7,z8(0)=mid8,z9(0)=mid9,z10(0)=mid10,z11(0)=mid11,z12(0)=mid12;

该微分方程可以比较简单,可以求得解析解。我也用4阶龙格库塔求数值解试了一下,结果差不多
> solution:=dsolve(,{y1(t),y2(t),y3(t),y4(t),y5(t),y6(t),z7(t),z8(t),z9(t),z10(t),z11(t),z12(t)});

subs(t=timestep,solution); 带入时间步长,求下一时间步的值

将解更新为下一步的初值mid1…mid12
od;

在得到位置和速度后还作了违约修正,使其误差几乎为零。之后把得到的位置连线画出来。发现摆的摆动幅度逐渐增大,不知是何原因,请指点迷津,万分感谢!

[ 本帖最后由 gily 于 2006-12-11 09:05 编辑 ]
页: [1]
查看完整版本: 运算双摆动力学时出现的问题,请指点迷津