二自由度强迫振动 微分方程求解

长风

原方程和参数以及结果如图
写的代码如下：

M代码

1. function dy=rigid(t,y)
   
2. %创建数组储存数据
   
3. dy=zeros(4,1); %原文，二元一次，现在分解应该四元一次了？
   
4. y(2)=dy(1);
   
5. y(4)=dy(3);
   
6. 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;
   
7. dy(2)=-(52.50*y(2)+52.5*y(1)-52.50*y(4)-52.5*y(3))/1814;
   
8. %运行代码
   
9. [t,y] = ode45('rigid',[0 1],[0.2 0 0.3 0]);
   
10. plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+',t,y(:,4),'.')

复制代码

应该四个曲线（机械上下部分的，位移和速度）都是接近正弦的，但计算的结果 不是。请问，我这错在哪里了？
另外，我想求，x1，x2及其一次二次导数和欧米伽（接近w那个）的关系，应该怎想写代码？

米斯兰达

首先，指出楼主给的图片中第二个方程错了。根据楼主的程序，这个方程左边最后两项应该是-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程序：

1. function dy=rigid(t,y)
   
2. dy=zeros(4,1);
   
3. dy(1)=y(2);
   
4. dy(2)=-(52.50*y(2)+52.5*y(1)-52.50*y(4)-52.5*y(3))/1814;
   
5. dy(3)=y(4);
   
6. 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;

复制代码

在命令窗口执行：

1. [t,y] = ode45('rigid',[0 250],[0.2 0 0.3 0]);
   
2. subplot(221)
   
3. plot(t,y(:,1),'-')
   
4. grid
   
5. subplot(222)
   
6. plot(t,y(:,2),'-')
   
7. grid;
   
8. subplot(223)
   
9. plot(t,y(:,3),'-')
   
10. grid
    
11. subplot(224)
    
12. plot(t,y(:,4),'-')
    
13. grid

复制代码

这里为了能更好的观察图形，将时间延长成了250，画出来的图如附件1.jpg所示。
这里y(:,1)、y(:,2)分别代表m1的位移和速度，而y(:,3)、y(:,4)分别代表m2的位移和速度。从图形中可以看出，由于楼主考虑了阻尼因素，所以曲线都是振荡衰减的。
另外，关于楼主“分解应该是四元一次”的疑问，实际上并不是四元，仍然还是两元（x1，x2），只不过为了能够在matlab里解这个微分方程组，做了一次变换而已。

长风

米斯兰达 发表于 2013-6-25 18:56

                               
登录/注册后可看大图

首先，指出楼主给的图片中第二个方程错了。根据楼主的程序，这个方程左边最后两项应该是-C1*x2'-K1*x2，而图 ...

我的式子中有，F0sinwt，这是一个持续的强制振动，这样还会衰减吗？

猫头鹰先生

我就是想问一下，楼主贴出的三个式子中的第二个式子后面出现了偏微分项，该如何求解？比如加一个Fx（x(t),y(t)）这一项，并不是恶意灌水，因为我做混沌里用到了偏微分，一直没有解出来。谢谢

长风

猫头鹰先生 发表于 2014-4-4 20:46
我就是想问一下，楼主贴出的三个式子中的第二个式子后面出现了偏微分项，该如何求解？比如加一个Fx（x(t),y ...

这个还没解决。我自己的代码错了，米斯兰大的代码，结果也不是我预期的。
从理论上说，这个无所谓偏微分不偏微分吧。高次或者多元的，化成多元一次的。只是具体方法我没找到合适的。

猫头鹰先生

长风 发表于 2014-4-11 10:30
这个还没解决。我自己的代码错了，米斯兰大的代码，结果也不是我预期的。
从理论上说，这个无所谓偏微分 ...

加油！                     

歪猴

luotao1014

{:{10}:}