运算双摆动力学时出现的问题，请指点迷津
我写了个计算一个普通双摆的动力学的代码，两杆的初始角度a1=30(degree),a2=45(degree)，从静止开始摆动。通过龙格库塔法求解动力学算出两杆中心的位置，但发现杆的摆动越来越大，不知是什么原因。动力学方程应该没错误，并且已经进行了位置和速度的违约修正，误差已经相当小了。但输出的曲线说明摆的幅度越来越大，如图红线为摆1的中心点在X轴的运动轨迹，绿线为杆二中心点在X轴的运动轨迹。请问这是怎么回事呢？
铰链连接的双摆
程序来自 《理论力学计算机模拟（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
你可以参考一下。 非常感谢pengweicai的热心帮助，但上面这段代码的主要语句即 =ode45('jjsbfun',,,[],l); 而我用的动力学方法是带拉格朗日乘子的动力学方程，所以求解的是微分代数方程组。但这两种方法结果应该差不多呀。我现在奇怪的是我的结果怎么会发散呢？已经进行了违约修正啊？
回复
建议先把代码传上来. 下面是双摆动力学在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]