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

gily

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

pengweicai

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

l=9;
   [t,u]=ode45('jjsbfun',[0:0.01:30],[0.5,0.2,0.1,2.8],[],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',[0.03 0.1 0.5 0.5]);%\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',[0.5 0.4 0.5 0.5]);%\fs{设置窗口坐标}
   cla;
   axis([-10 10 -20 20])
   axis equal
   hold on

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


ball1a=line(x1(1),y1(1),'color',[0.5 0.6 0.4],'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',[0.5 0.6 0.4],'linestyle','-','linewidth',1,...
   'erasemode','none');
ball2=line(x2(1),y2(1),'color','r','marker','.','markersize',40,'erasemode','xor');
gan1=line([0,x1(1)],[0,y1(1)],'color','g','linewidth',2,'erasemode','xor');
gan2=line([x1(1),x2(1)],[y1(1),y2(1)],'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',[0,x1(i)],'ydata',[0,y1(i)]);
   set(gan2,'xdata',[x1(i),x2(i)],'ydata',[y1(i),y2(i)]);
   drawnow
end


  你可以参考一下。

gily

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

xjzuo

建议先把代码传上来.

gily

下面是双摆动力学在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([[m,0,0,0,0,0,1,0,0,0],[0,m,0,0,0,0,0,1,0,0],[0,0,J,0,0,0,-line_L0*cos(phi1),-line_L0*sin(phi1),-2*line_L0*cos(phi1),-2*line_L0*sin(phi1)],[0,0,0,m,0,0,0,0,1,0],[0,0,0,0,m,0,0,0,0,1],[0,0,0,0,0,J,0,0,-line_L0*cos(phi2),-line_L0*sin(phi2)],[1,0,-line_L0*cos(phi1),0,0,0,0,0,0,0],[0,1,-line_L0*sin(phi1),0,0,0,0,0,0,0],[0,0,-2*line_L0*cos(phi1),1,0,-line_L0*cos(phi2),0,0,0,0],[0,0,-2*line_L0*sin(phi1),0,1,-line_L0*sin(phi2),0,0,0,0]]);

> 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[1];
> eqn8:=diff(z8(t),t)=linearsolution[2];
> eqn9:=diff(z9(t),t)=linearsolution[3];
> eqn10:=diff(z10(t),t)=linearsolution[4];
> eqn11:=diff(z11(t),t)=linearsolution[5];
> eqn12:=diff(z12(t),t)=linearsolution[6];

初值
> 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([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,eqn11,eqn12,ICs],{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 编辑 ]