一个非线性粘弹性梁的Poincare截面图
这是刘延柱、陈立群书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);
=ode45(@ex6_1_5_01,,,options);
figure(1)
plot(X(1:100:end,1),X(1:100:end,2),'.b')
选取的周期为2*pi/w，这里w=1，所以每隔2*pi的时长截取一个点。
按照书中的讲述，Poincare截面应该是一个准周期环面，如下图所示：
但用我的程序画出的是这样的：
不知道错在哪里，请大家帮忙，非常感谢！
回复 1 # john152 的帖子
感觉计算时间有点不够，很难进入稳态响应。庞加莱截面需要去除瞬态解。
plot(X(1:100:end,1),X(1:100:end,2),'.b')包含了很多的顺态解。
建议积分时间取长，去除前面迭代的瞬态解，只取稳态画图！
meiyongyuandeze 发表于 2011-4-21 09:45 static/image/common/back.gif
回复 1 # john152 的帖子
感觉计算时间有点不够，很难进入稳态响应。庞加莱截面需要去除瞬态解。
感觉还是不行，我用了2000个周期画图：
beta=0.0000001;
tt=2*pi/w;
options=odeset('RelTol',1e-10);
=ode45(@ex6_1_5_01,,,options);
figure(1)
plot(X(100000:100:end,1),X(100000:100:end,2),'.b')
截取了周面一半的周期点，但画出的还是这个样子：
我定义的状态变量没有问题吧？
回复 3 # john152 的帖子
状态空间是正确的，我看了下你给的原图q的范围在[-6,3]之间，而楼主的在[-4,3]，如果对于的系统式相同的话，应该不会出现这么大的差别。是不是初值的问题。。。。。。。。。 回复 4 # meiyongyuandeze 的帖子
应该不是初值的问题，我试了其他初值，都是这个样子 那很难说是怎么回事了......，有时候也要坚信自己的结果........ meiyongyuandeze 发表于 2011-4-21 17:57 static/image/common/back.gif
那很难说是怎么回事了......，有时候也要坚信自己的结果........
嗯，谢谢，我用LET工具箱算了一下，beta这样取值时，它的Lyapunov指数也有大于0的。 呵呵，那些书本的漂亮的东西可能是个“传说”，呵呵！
页:
[1]