二自由度强迫振动 微分方程求解
本帖最后由 牛小贱 于 2014-4-4 20:34 编辑原方程和参数以及结果如图
写的代码如下：
M代码function dy=rigid(t,y)
%创建数组储存数据
dy=zeros(4,1); %原文，二元一次，现在分解应该四元一次了？
y(2)=dy(1);
y(4)=dy(3);
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;
dy(2)=-(52.50*y(2)+52.5*y(1)-52.50*y(4)-52.5*y(3))/1814;
%运行代码
= ode45('rigid',,);
plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+',t,y(:,4),'.')
应该四个曲线（机械上下部分的，位移和速度）都是接近正弦的，但计算的结果不是。请问，我这错在哪里了？
另外，我想求，x1，x2及其一次二次导数和欧米伽（接近w那个）的关系，应该怎想写代码？
本帖最后由 牛小贱 于 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;在命令窗口执行：
= ode45('rigid',,);
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里解这个微分方程组，做了一次变换而已。 米斯兰达 发表于 2013-6-25 18:56 static/image/common/back.gif
首先，指出楼主给的图片中第二个方程错了。根据楼主的程序，这个方程左边最后两项应该是-C1*x2'-K1*x2，而图 ...
我的式子中有，F0sinwt，这是一个持续的强制振动，这样还会衰减吗？ 我就是想问一下，楼主贴出的三个式子中的第二个式子后面出现了偏微分项，该如何求解？比如加一个Fx（x(t),y(t)）这一项，并不是恶意灌水，因为我做混沌里用到了偏微分，一直没有解出来。谢谢 猫头鹰先生 发表于 2014-4-4 20:46
我就是想问一下，楼主贴出的三个式子中的第二个式子后面出现了偏微分项，该如何求解？比如加一个Fx（x(t),y ...
这个还没解决。我自己的代码错了，米斯兰大的代码，结果也不是我预期的。
从理论上说，这个无所谓偏微分不偏微分吧。高次或者多元的，化成多元一次的。只是具体方法我没找到合适的。 长风 发表于 2014-4-11 10:30
这个还没解决。我自己的代码错了，米斯兰大的代码，结果也不是我预期的。
从理论上说，这个无所谓偏微分 ...
加油！ {:{10}:} {:{10}:}{:{10}:}{:{10}:}{:{10}:}
页:
[1]