本帖最后由 牛小贱 于 2014-4-16 19:30 编辑
首先,指出楼主给的图片中第二个方程错了。根据楼主的程序,这个方程左边最后两项应该是-C1*x2'-K1*x2,而图片中的方程相应位置却成了-C1*x1'-K1*x1,这是错误的。
其次,检查楼主的程序,怀疑是
y(2)=dy(1);
y(4)=dy(3);
这两项写错了(顺序反了),正确的写法应该是:
dy(1)=y(2);
dy(3)=y(4);
这样,改正后完整的程序如下:
M程序:
- function dy=rigid(t,y)
- dy=zeros(4,1);
- dy(1)=y(2);
- dy(2)=-(52.50*y(2)+52.5*y(1)-52.50*y(4)-52.5*y(3))/1814;
- dy(3)=y(4);
- dy(4)=(60*sin(30*t)- (52.50+700.5)*y(4)-(52.5+140.1)*y(3)+52.50*y(2)+52.5*y(1))/2903;
复制代码 在命令窗口执行:
- [t,y] = ode45('rigid',[0 250],[0.2 0 0.3 0]);
- subplot(221)
- plot(t,y(:,1),'-')
- grid
- subplot(222)
- plot(t,y(:,2),'-')
- grid;
- subplot(223)
- plot(t,y(:,3),'-')
- grid
- subplot(224)
- plot(t,y(:,4),'-')
- grid
复制代码 这里为了能更好的观察图形,将时间延长成了250,画出来的图如附件1.jpg所示。
这里y(:,1)、y(:,2)分别代表m1的位移和速度,而y(:,3)、y(:,4)分别代表m2的位移和速度。从图形中可以看出,由于楼主考虑了阻尼因素,所以曲线都是振荡衰减的。
另外,关于楼主“分解应该是四元一次”的疑问,实际上并不是四元,仍然还是两元(x1,x2),只不过为了能够在matlab里解这个微分方程组,做了一次变换而已。 |