求助拉索参数振动，以及相图和庞加莱映射

petterchan

斜拉索的受轴向激励的参数振动分析，得到非线性振动微分方程，进行数值积分求解，用4解龙格-库塔法，绘出其跨中位移时程，相图，庞加莱映射图分别如下
但是不知道从该庞加莱图能得到什么结论。另外，因为该参数振动系统瞬时振动时段存在‘拍振’，想分析该非线性系统参数对拍振阶段的影响，不如拉索阻尼，垂度，倾斜角，但是不知如何下手。
请问大侠瞬态阶段分析和稳态阶段分析有什么好方法没有？请大家指教啊，毕业着急啊时间紧迫
另外，程序及非线性微分方程见附件
%下面方程为微分方程程序，命名为hnupr.m保存在matlab的工作路径内
function pdot = hnupr(t,p,flag,a1,a2,a3,a4,w1,ita,niu)
% 输入参数为各种初始条件的积分参数
% a1;   %平方非线性项系数
% a2;   %立方非线性项
% a3;   %与激励有关的刚度变化系数
% a4;   %激励系数，与拉索倾斜角、激励幅值、弦向分力有关
% ita;  %阻尼系数
% w0;   %标准弦的圆频率的平方
% w1;   %考虑拉索的垂度效应的拉索圆频率的平方
% Le;   %拉索总长度
% lamd; %Irvine常数
% 2*niu;     %激励的频率
% g = sin(2*niu*t); % forcing function
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
pdot = zeros(size(p));
pdot(1) = p(2);
pdot(2) = -ita*p(2)-(w1+a3*cos(2*niu*t))*p(1)-a1*p(1)^2-a2*p(1)^3-a4*cos(2*niu*t);
% eof - hnupr.m

%以下为主程序，直接复制粘贴到matlab的命令窗口就行了
clear
clc
theta=14.2; %拉索的倾斜角，单位为度
L=22.21;    %拉索的弦线长度，单位为米
E=190e9;    %弹性模量，单位为GPa
A=0.0001374;  %拉索横截面积，单位为m^-2
m=11.01;    %拉索的单位长度质量，单位为kg/m
g=9.81;     %重力加速度，单位为m/s^2
N=42800;    %初始张力，单位为N
H=34000;    %拉索静张力的弦向分量,单位为N
Le=L*(1+(m*g*L*cos(theta*pi/180)/H)^2/8);       %拉索的总长度
lamd=(m*g*L*cos(theta*pi/180)/H)^2*L/(H*Le/(E*A));  %Irvine常数的平方
w0=(pi^2*H)/(m*L^2);    %拉索的一阶振动圆频率的平方
w1=w0*(1+(2/pi)^4*lamd/2);      %考虑拉索的垂度和倾斜角的自振圆频率的平方
Xd=0.001;   %激励的幅值,单位为m
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
ita=2*0.002255*sqrt(w1);    %ita为衰减系数，一阶阻尼比为0.2255%,二阶的阻尼比为0.2088%,前面的数字为阻尼比，动力学上为ita=2*阻尼比*w1;
niu=sqrt(w1);  %niu为激励频率的二分之一，激励频率为圆频率
a1=3*pi*E*A*g*cos(theta*pi/180)/(H*L^2);    %平方非线性项系数
a2=pi^4*E*A/(4*m*L^4);  %立方非线性项
a3=pi^2*E*A*Xd/(m*L^3); %与激励有关的刚度变化系数
a4=4*E*A*g*Xd*cos(theta*pi/180)/(pi*H*L);   %激励系数，与拉索倾斜角、激励幅值、弦向分力有关
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
tspan=[0:0.01:500];   %积分的时间长度，可分为变步长和固定步长
initial_condition=[0 0];  %初始条件
[t,p] = ode45('hnupr',tspan,initial_condition,[],a1,a2,a3,a4,w1,ita,niu);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

%绘庞加莱图的程序如下，参照本论坛一个程序的代码，应用那个据说是频闪法的，也不是很了解，就是这个地方比较烦恼
m=zeros(500,1);
n=zeros(500,1);
for i=1:500
    m(i,1)=p(10*i,1);
    n(i,1)=p(10*i,2);
end
plot(m,n,'k.','markersize',10);

[ 本帖最后由 petterchan 于 2008-9-24 17:05 编辑 ]

petterchan

怎么没人回答啊？自己顶一下

无水1324

作分岔、Poincare等，都是对稳态进行分析，这时候必须消除瞬态！