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

gily

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

下面是双摆动力学在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;

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

suffer

能否贴一下你的动力学模型？