请教动力学方程求解

zjplay · 发表于 2011-4-12 09:59

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

请教方程为何无法求解？

clear all;
clf;
fs=1000;
n=2048;
tspan=0:1/fs:(n-1)/fs;
x0 = [1; 1; 0; 0; 0; 0; 0; 0];
[t,x] = ode45('solution', tspan, x0);
plot(t,x(:,1));
xlabel('t');
ylabel('x(t)');
title('时域图');

function f = solution(t,x)
rb1=0.1;
rb2=0.2;
T1=10;
T2=10;
e=0.0001;
omega2=1;
kg=1e6;
cg=0.08;
J1=1.25e-3;
J2=2e-2;
m2=1;
k2=3.55e6;
c2=0.15;
f = zeros(8,1);
f(1) = x(2);
f(2) = (T1-(kg*(rb2*x(3)-rb1*x(1)+x(7))+cg*(rb2*x(4)-rb1*x(2)+x(8)))*rb1)/J1;
f(3) = x(4);
f(4) = (-T2-(kg*(rb2*x(3)-rb1*x(1)+x(7))+cg*(rb2*x(4)-rb1*x(2)+x(8)))*rb2)/J2;
f(5) = x(6);
f(6) = -m2*e*(omega2)^2*cos(omega2*t)-c2*x(6)-k2*x(5);
f(7) = x(8);
f(8) = -m2*e*(omega2)^2*sin(omega2*t)-c2*x(8)-k2*x(7)-(kg*(rb2*x(3)-rb1*x(1)+x(7))+cg*(rb2*x(4)-rb1*x(2)+x(8)));

zjplay · 发表于 2011-4-12 10:13

补充下问题，系统是线性的，本来想考察下外激励的影响，将T1写作T1+a*sin(omega2*t)，但仍然无法求解。问题出来哪里呢？

meiyongyuandeze · 发表于 2011-4-12 10:41

回复 2 # zjplay 的帖子

你的无解是什么意思？是循环不出来还是求出的解不收敛呢，能具体点吗？

zjplay · 发表于 2011-4-12 11:13

求出来的解是NaN

zjplay · 发表于 2011-4-12 11:15

在猜想是不是没有量纲化的原因呢？如果没有量纲化，会出现这样的问题吗？还是程序本身有问题呢？

meiyongyuandeze · 发表于 2011-4-12 11:30

回复 4 # zjplay 的帖子

我刚才大概运行了下，发现你的状态方程中f的值在时间0.1s后趋近与零，这就说明你的方程中出现了分母无穷大的情况，估计是你的方程有问题，请仔细检查

zjplay · 发表于 2011-4-12 14:31

能否留下联系方式，将方程给你看看。

meiyongyuandeze · 发表于 2011-4-12 14:57

回复 7 # zjplay 的帖子

把方程贴到论坛上吧，可能有人对你的问题比较熟悉，帮助更大！

zjplay · 发表于 2011-4-12 15:20

不知道能否看清？

zjplay · 发表于 2011-4-12 15:36

还有个问题，计算的poincare图怎么会与相图一样？
方程的w=1, 应该如何定间隔才合适呢？

meiyongyuandeze · 发表于 2011-4-12 15:47

回复 10 # zjplay 的帖子

感觉因该是拟周期吧，庞加莱映射取点太少了，建议夺取几个点，还要注意的是只取稳态点进行就算！

zjplay · 发表于 2011-4-12 16:10

但是我计算的谱图感觉应该不是拟周期呢

zjplay · 发表于 2011-4-12 16:14

都是截取的稳态点计算的，不知道问题在哪里？
感觉poincare图不应该是这样的。

zjplay · 发表于 2011-4-12 16:16

回复 11 # meiyongyuandeze 的帖子

刚才贴出的方程以及求解程序是否有问题？

zjplay · 发表于 2011-4-12 16:27

根据FFT的结果，周期该如何取呢？
这点儿有点不清楚，请教了。

账号		自动登录	找回密码
密码			我要加入

[自激振动] 请教动力学方程求解

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

点评

点评