一个非线性粘弹性梁的Poincare截面图

john152

这是刘延柱、陈立群书p257页的一个例子
一类几何非线性粘弹性梁运动的简化运动力学模型：

这是一个三阶微分方程，我定义方程状态变量为；
function dx=ex6_1_5_01(t,x)
global beta w e b a alpha F
dx=[x(2);
    x(3);
    -beta*x(3)-w^2*(1-e*cos(w*t))*x(2)-w^2*(e*w*sin(w*t)+beta*a*(1-e*cos(w*t)))*x(1)-...
    3/8*w^2*alpha^2*(1-3*e*cos(w*t))*x(1)^2*x(2)-1/8*w^2*alpha^2*(3*w*e*sin(w*t)+...
    beta*a-beta*(1-3*b)*e*cos(w*t))*x(1)^3+F*(w*sin(w*t)-beta*cos(w*t))];

绘制Poincare截面的程序为：
clear all;clc;
global beta w e b a alpha F
a=0.1;b=0.9;w=1.0;alpha=2.8284;e=0.01;F=34.4964;
beta=0.0000001;
tt=2*pi/w;
options=odeset('RelTol',1e-10);
[T,X]=ode45(@ex6_1_5_01,[0:tt/100:500*tt],[0.1,0,0],options);
figure(1)
plot(X(1:100:end,1),X(1:100:end,2),'.b')

选取的周期为2*pi/w，这里w=1，所以每隔2*pi的时长截取一个点。
按照书中的讲述，Poincare截面应该是一个准周期环面，如下图所示：

但用我的程序画出的是这样的：


不知道错在哪里，请大家帮忙，非常感谢！

meiyongyuandeze

回复 1 # john152 的帖子

感觉计算时间有点不够，很难进入稳态响应。庞加莱截面需要去除瞬态解。
plot(X(1:100:end,1),X(1:100:end,2),'.b')包含了很多的顺态解。
建议积分时间取长，去除前面迭代的瞬态解，只取稳态画图！

john152

meiyongyuandeze 发表于 2011-4-21 09:45

                               
登录/注册后可看大图

回复 1 # john152 的帖子

感觉计算时间有点不够，很难进入稳态响应。庞加莱截面需要去除瞬态解。

感觉还是不行，我用了2000个周期画图：
beta=0.0000001;
tt=2*pi/w;
options=odeset('RelTol',1e-10);
[T,X]=ode45(@ex6_1_5_01,[0:tt/100:2000*tt],[0.1,0,0],options);
figure(1)
plot(X(100000:100:end,1),X(100000:100:end,2),'.b')

截取了周面一半的周期点，但画出的还是这个样子：

我定义的状态变量没有问题吧？

meiyongyuandeze

回复 3 # john152 的帖子

状态空间是正确的，我看了下你给的原图q的范围在[-6,3]之间，而楼主的在[-4,3]，如果对于的系统式相同的话，应该不会出现这么大的差别。是不是初值的问题。。。。。。。。。

john152

回复 4 # meiyongyuandeze 的帖子

应该不是初值的问题，我试了其他初值，都是这个样子

meiyongyuandeze

那很难说是怎么回事了......，有时候也要坚信自己的结果........

john152

meiyongyuandeze 发表于 2011-4-21 17:57

                               
登录/注册后可看大图

那很难说是怎么回事了......，有时候也要坚信自己的结果........

嗯，谢谢，我用LET工具箱算了一下，beta这样取值时，它的Lyapunov指数也有大于0的。

meiyongyuandeze

呵呵，那些书本的漂亮的东西可能是个“传说”，呵呵！